cannin/biostars_qa huggingface dataset:
Question: Example input: multi-sample VCF (adapted from [www.internationalgenome.org][1]):

**Note:** my actual file is bgzipped and tabixed with ~2Mln variants (rows) and ~1000 samples (columns).

##fileformat=VCFv4.0
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=1000GenomesPilot-NCBI36
##phasing=partial
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##INFO=
##FILTER=
##FILTER=
##FORMAT=
##FORMAT=
##FORMAT=
##FORMAT=
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT NA00001 NA00002 NA00003
20 14370 rs6054257 G A 29 PASS NS=3;DP=14;AF=0.5;DB;H2 GT:GQ:DP:HQ 0/0:48:1:51,51 1/0:48:8:51,51 1/1:43:5:.,.
20 17330 . T A 3 q10 NS=3;DP=11;AF=0.017 GT:GQ:DP:HQ 0/0:49:3:58,50 0/1:3:5:65,3 0/0:41:3
20 1110696 rs6040355 A G,T 67 PASS NS=2;DP=10;AF=0.333,0.667;AA=T;DB GT:GQ:DP:HQ 1/2:21:6:23,27 2/1:2:0:18,2 2/2:35:4
21 1230237 . T . 47 PASS NS=3;DP=13;AA=T GT:GQ:DP:HQ 0/0:54:7:56,60 0/0:48:4:51,51 0/0:61:2
21 1234567 microsat1 GTCT G,GTACT 50 PASS NS=3;DP=9;AA=G GT:GQ:DP 0/1:35:4 0/2:17:2 1/1:40:3

Expected output:

**20.txt**

NA00001	NA00002	NA00003
0/0	1/0	1/1
0/0	0/1	0/0
1/2	1/2	2/2

**21.txt**

NA00001	NA00002	NA00003
0/0	0/0	0/0
0/1	0/2	1/1

I was thinking of using `cut | sed` combo with some `regex`, but thought there must be already some tool out there, maybe *[bcftools][2]* (couldn't get the right flags to work) ?

Any other ideas?

[1]: http://www.internationalgenome.org/wiki/Analysis/Variant%20Call%20Format/vcf-variant-call-format-version-40/
[2]: http://www.htslib.org/doc/bcftools.html

Answer: bcftools annotate -x '^FORMAT/GT' input.vcf.gz | grep -v "^##" | cut -f 10-
Question: This dataset is for a specific disease-gene-test results. The dataset goes like this.

| Test | Gene | Relevance | Values |

Test and Gene are the two parameters on x and y axis. Values are the combined results for the parameters pair.
The problem for me here is that I have one more parameter called Relevance which represents the relevance of the test-gene pair and it is boolean (only two values YES/NO).

The dataset (Relevance) should be differentiated in the map with different colour (like red and green).

Gradient of that colour represents numerical values between that interaction.

The end result I was aiming at was, Test on x-axis and Gene on y-axis and the map for these interaction will be only with two colours (representing Relevance values) and the gradient of that colour representing Values.
Is this possible to achieve this kind of Heatmap, if Yes how can I achieve.
If not is there any other option to display such kind of data (similar to heatmap)

Help appreciated !!

Thanks,

RDS

something like this -

Answer: I am sure you can do it in R, a little bit fancy using the ggplot2 library and having your data in the form of dataframe.
Check these two posts on how to achieve it.
ggplot2-quick-heatmap-plotting and Constructing an Heatmap of "Distance of binding region relative to TSS".

Cheers
Question: I have a series of fastq files (with up to 4000 reads in each) that I want to parse based on the time of sequencing.

So in the fastq header, date/time is listed as "start_time=2017-10-09T18:54:24z"

If I wanted to extract all sequences between 18:00 hours and 20:00 hours, is there a tool I can use to find and extract them?

Answer: I wrote a python script now.

Execute as
`python timefilt.py part1.fastq.gz --time_from 2017-10-09T18:00:00Z --time_to 2017-10-09T20:00:00Z`

https://gist.github.com/wdecoster/1ab9adac7c8095498ff91ee22468eaac
Question: Heys,

I imagine this is an odd question for most of you, but for me would be really useful to know your opinion on this. I already have some experience with bash, I know the basic commands and I even run for loops and make some functions, but now I achieved a point were I have the impression I can't go forward without help. However, I've been trying to look for tutorials or whatever option that teaches me a bit further without success. I'm in datacamp, but they don't have a lot of courses on bash (although the ones they have helped me quite a lot) and in internet I did not find any practical help.

So my question is, where could I find some tutorials/information about more advanced help in bash? I still need to learn a lot and I do have a lot of (i imagine) simple doubts (for example, is it better to use a for loop or a function, how to optimize bash scripts, ...), but without an external help I don't know how to answer them. What did you do where you achieved my situation? I'm a biologist without background on informatics, so I imagine that this not help in my situation. Whatever option you think will be more than welcome!

Thanks!!

Answer: The Linux Foundation has some advanced tutorials https://tldp.org/LDP/abs/html/index.html
Question: Hi there,
I'm using bowtie2 to align some fastq files downloaded from [ENA][1]

I've first run a quality check to see if there is something annoying, but in general terms it is fine.
Then, I've run bowtie-2 as follow:

bowtie2 -x human_index_bowtie/human_index_bwt -U SRR11557616.fastq.gz -S 5h.sam

But I got this error message:

Error: Read SRR11557616.1 ugc_599_6_F3_0002_0/1 has more quality values than read characters.
terminate called after throwing an instance of 'int'
Aborted (core dumped)
(ERR): bowtie2-align exited with value 134

the head of the file is a bit weird:

head -n 8

@SRR11557614.1 ugc_599_3_F3_0001_0 length=74
T3.0321.1000.133...11...011..20...11...21...10...11...31....0....2....1....
+SRR11557614.1 ugc_599_3_F3_0001_0 length=74
!aB_aaaBa^aQBWaYBBBaYBBBOa_BBWSBBBKaBBBNWBBBWOBBBS^BBBN`BBBBSBBBBPBBBBaBBBB
@SRR11557614.2 ugc_599_3_F3_0001_1 length=74
T2.3313.1031.030...30...231..10...22...00...02...22...13....0....2....2....
+SRR11557614.2 ugc_599_3_F3_0001_1 length=74
!aBaaaaBaa_]BPa^BBBbaBBBaaaBBaaBBB\`BBBa\BBBaaBBBaaBBBWPBBBBaBBBBaBBBBaBBBB

Any ideas how to solve the problem? or I should select a different file?

Thanks!

[1]: https://www.ebi.ac.uk/ena/browser/home

Answer: The problem is that this is a colorspace file (like ABI SOLiD system, not Illumina), and `bowtie2` does not support colorspace alignment. `bowtie` (aka `bowtie1`) used to support it, but support got dropped in recent versions afaik.

See: http://bowtie-bio.sourceforge.net/news.shtml

Colorspace support was dropped in versions 1.3.0 and later so you would need to download a version prior to that, then build an index `bowtie-build` with the `-C` option to enable colorspace, and then align your data. See manual for details, I never used colorspace myself.
Question: Hallo,

I would like to create such an expression matrix:

GENE ID/SAMPLE sample1 sample2 sample3 sample4 sample-n
gene1 logFC ... ... ... ...
gene2 ... ... ... ... ...
gene3 ... ... ... ... ...
gene4 ... ... ... ... ...
gene-n ... ... ... ... ...

using data from GSE at GEO. I will be grateful for any suggestions or just codes.

Second thing:

How to get such an expression matrix after limma DE analysis? How to modify limma codes? I will appreciate any suggestions or codes.

Thanks in advance, one more time.

Regards!

Answer: Something like this should about do it.

library(GEOquery)
# assumes only one platform in the GSE
gse = getGEO('GSEXXXX')[[1]]

gse is an ExpressionSet, so all the usual ExpressionSet methods work as expected. In particular, we can use `fData()` to get the feature information (gene information) and the `exprs()` method to get the actual values.

write.table(data.frame(fData(gse),exprs(gse)),sep="\t",row.names=FALSE,file='abc.txt')

You may need to use a subset of the columns in `fData` to match your needs.
Question: Is there a way to sample sequences from the `fastq.gz` file? I have files that contain hundreds of millions of reads, I want to sample about 100 million reads from each file. Is there a way to do this?

Answer: [BBTools's][1] reformat tool does this, keeping pairs together if they are interleaved:

reformat.sh in=reads.fq.gz out=sampled.fq.gz samplereadstarget=100000000

or for reads in 2 files:

reformat.sh in1=reads1.fq.gz in2=reads2.fq.gz out=sampled1.fq.gz out2=sample2.fq.gz samplereadstarget=100000000

You can also just sample a fraction of the reads without respect to the final number, which is faster since it only needs to read the file once:

reformat.sh in=reads.fq.gz out=sampled.fq.gz samplerate=0.1

or just sample the first X reads, which is faster still:

reformat.sh in=reads.fq.gz out=sampled.fq.gz reads=100000000

Reformat is multithreaded and extremely fast.

[1]: https://sourceforge.net/projects/bbmap/
Question: I am using IGV to look at some RNA-seq alignment. In general I find this tool extremely helpful. However, when looking at some longer genes I am unable to visualize the alignment and I receive the message: "Zoom in to see alignment".
Is there a way to change this setting to visualize the entire gene (maybe by increasing memory usage of the tool?). Does anybody know how to do this?

Answer: Go to `View > Preferences`. Then in `Alignments` tab, change the `Visibility range threshold` field.
Question: Hello:

I would like to convert .gff3 file to 12-column .bed file, as in this link under "BED Format" (http://genome.ucsc.edu/FAQ/FAQformat.html#format1).

I have thus far used Galaxy from Penn State, but it outputs a 6-column .bed file.

Any advice is greatly appreciated!

Thank you

Answer: The UCSC Genome Browser [hosts conversion utilities][1] that you can run from your command line to accomplish the gff3 to BED12 conversion. Note utilities are OS specific and need to be given permission to execute with "chmod +x utilityName".

Here's an example of how I did a conversion using the following steps:

- wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_31/gencode.v31.basic.annotation.gff3.gz
- gunzip gencode.v31.basic.annotation.gff3.gz
- ./gff3ToGenePred
gencode.v31.basic.annotation.gff3 gencode.v31.basic.genePred
- ./genePredToBed gencode.v31.basic.genePred gencode.v31.basic.bed

Disclaimer that I work for the UCSC Genome Browser. :)

[1]: http://hgdownload.soe.ucsc.edu/downloads.html#utilities_downloads
Question: I have a Mac with OS 10.11.4, which has both Python 2.7.3 (OS Build-in) and Python 3.5.1 (just installed). With Python2, I can import Biopython (version 1.61). The Python3, however, cannot see the Biopython package. I want to update Biopython to the latest version 1.66 and link it to Python3. Do I need to uninstall version 1.61 first? Is yes, how? How to make sure the installation from tarball will link it to Python3? Thanks.

Answer: This is somewhat off-topic for the site, but if you `pip3 install --user biopython` it'll install the most recent version of biopython. Packages are always specific to a version of python, so your python2 and python3 installed versions can coexist.
Question: Hi,

I've a GRanges representing a intervals for all genes in the genome. A lot of these intervals are overlapping. I would like to use the `reduce()` from `GenomicRanges` package in order to make a non-overlapping set of interval. However I would like to do it for each gene separately. Thus for one specific gene, intervals for this gene should not overlap ; but intervals for different genes may overlap. One solution would like to split the GRanges by gene and apply reduce() on each subset but I'm wondering if there is a more efficient way ?

Thanks

Actual data

chrom	start	end	hgnc
1	100	200	MYC
1	150	300	MYC
1	400	500	MYC
1	150	230	TP53
1	200	350	TP53
1	420	550	TP53

expected result

chrom	start	end	hgnc
1	100	300	MYC
1	400	500	MYC
1	150	350	TP53
1	420	550	TP53

My actual solution :

# gene is the dataframe used to create the initial GRanges
do.call(rbind,lapply(
split(gene,gene$hgnc),
function(x){
as.data.frame(
reduce(
GRanges(x$chrom,IRanges(x$start,x$end))))}))

Answer: Looks like we can't avoid *reduce()*, but we can avoid `do.call(rbind, lapply(...` using *data.table*:

library(data.table)
library(GenomicRanges) # reduce function

# example data
x <- fread("
chrom start end hgnc
1 100 200 MYC
1 150 300 MYC
1 400 500 MYC
1 150 230 TP53
1 200 350 TP53
1 420 550 TP53")

x[, as.data.table(reduce(IRanges(start, end))), by = .(chrom, hgnc)]
# chrom hgnc start end width
# 1: 1 MYC 100 300 201
# 2: 1 MYC 400 500 101
# 3: 1 TP53 150 350 201
# 4: 1 TP53 420 550 131

Code stolen and adapted from StackOverflow post:

- [Finding all overlaps in one iteration of foverlap in R's data.table](https://stackoverflow.com/a/45574386/680068)
Question: How can separate one column in two column in R?

id rs143 rs148 rs149 rs1490
1 02003s NA 11 22 11
2 02003s NA 10 11 22
3 02003s NA 11 11 12
4 02003s NA 10 11 11
5 02003s NA 10 11 11

in result i want this format:

id rs143 rs143.1 rs148 rs148.1 rs149 rs149.1 rs1490 rs1490.1
1 02003s NA NA 1 1 2 2 1 1
2 02003s NA NA 1 0 1 1 2 2
3 02003s NA NA 1 1 1 1 1 2
4 02003s NA NA 1 0 1 1 1 1
5 02003s NA NA 1 0 1 1 1 1

Answer: Requires a nice mixture of diverse functions here. This could be done in a single line, but would be complex.

> df
id rs143 rs148 rs149 rs1490
1 02003s NA 11 22 11
2 02003s NA 10 11 22
3 02003s NA 11 11 12
4 02003s NA 10 11 11
5 02003s NA 10 11 11

#Ensure that NAs are encoded as characters
> df[is.na(df)] <- "NA"

#split each value by an empty string delimiter, then re-merge all columns back together
> df2 <- data.frame(df$id, do.call(cbind, lapply(df[,2:ncol(df)], function(x) t(as.data.frame(strsplit(as.character(x), split=""))))), row.names=c(1:nrow(df)))

> df2
df.id X1 X2 X3 X4 X5 X6 X7 X8
1 02003s N A 1 1 2 2 1 1
2 02003s N A 1 0 1 1 2 2
3 02003s N A 1 1 1 1 1 2
4 02003s N A 1 0 1 1 1 1
5 02003s N A 1 0 1 1 1 1

#Now fix the colnames
> index1 <- seq(from=2, to=ncol(df2), by=2)
> index2 <- seq(from=3, to=ncol(df2), by=2)
> colnames(df2)[index1] <- colnames(df[2:ncol(df)])
> colnames(df2)[index2] <- paste(colnames(df[2:ncol(df)]), ".1", sep="")

> df2
df.id rs143 rs143.1 rs148 rs148.1 rs149 rs149.1 rs1490 rs1490.1
1 02003s N A 1 1 2 2 1 1
2 02003s N A 1 0 1 1 2 2
3 02003s N A 1 1 1 1 1 2
4 02003s N A 1 0 1 1 1 1
5 02003s N A 1 0 1 1 1 1

#Restore he NAs
> df2[df2=="N"] <- NA
> df2[df2=="A"] <- NA

df
id rs143 rs148 rs149 rs1490
1 02003s NA 11 22 11
2 02003s NA 10 11 22
3 02003s NA 11 11 12
4 02003s NA 10 11 11
5 02003s NA 10 11 11

df2
df.id rs143 rs143.1 rs148 rs148.1 rs149 rs149.1 rs1490 rs1490.1
1 02003s 1 1 2 2 1 1
2 02003s 1 0 1 1 2 2
3 02003s 1 1 1 1 1 2
4 02003s 1 0 1 1 1 1
5 02003s 1 0 1 1 1 1
Question: Hi,

I was wondering if anyone knows if there is a default protein database linked to the MAF files of the TCGA project?

Basically what I want to do is to get the mutated protein sequences that correspond to the missense mutations listed in the MAF files. For that, I need the correct protein sequence for each gene with a missense mutation, since the MAF file has the position of the mutation in the gene and in the protein, I can write a simple script that changes the wild type amino acid with the mutation. However, it is crucial to get the correct protein sequence, so that the position and wild type amino acid stored in the MAF file corresponds to the same amino acid in that position of the protein sequence in the database.

Thanks

Answer: The latest TCGA MAF standard is to choose the "[worst affected][1]" isoform per variant, from among the Gencode Basic v19 isoforms. These standards have changed over time, so not all MAFs will use the same isoform database. You'll find GAF files listed [here][2], originally based on UCSC KnownGenes, but now based on Gencode Basic v19.

If you want to map each variant to [Uniprot's canonical isoform][3] per gene, then pull the Ensembl ENST IDs of all Uniprot's canonical isoforms, dump them into a text file one ID per line, and pass it to [maf2maf][4] under argument `--custom-enst`, to re-annotate all TCGA MAFs.

[1]: http://useast.ensembl.org/info/genome/variation/predicted_data.html#consequences
[2]: https://tcga-data.nci.nih.gov/docs/GAF/
[3]: http://www.uniprot.org/help/canonical_and_isoforms
[4]: https://github.com/ckandoth/vcf2maf
Question: I am reading a paper about miRNA and confused by a fig.
![enter image description here][1]

[1]: http://i.imgur.com/yKXw8Ja.png

In "Novel miRNA for which miRNA* were found", what does miRNA* mean in the context?
Does it mean reverse-complement sequence of miRNA?

I hope someone could give me some suggestions.

Answer: miRNA* refers to the antisense sequence to the mature miRNA. miRDeep classifies miRNA reads based on the location of reads in the pre-miRNA: miRNA (mature miRNA sequence), miRNA* (as sequence of miRNA), and loop (matching loop region of pre-miRNA hairpin).
Question: The Homo sapiens Genome (Build 37.2) has just been published by the NCBI : http://www.ncbi.nlm.nih.gov/mapview/stats/BuildStats.cgi?taxid=9606&build=37&ver=2

Do the Fasta sequences of the human genome have changed since the version 37.1 (e.g: some large contigs of NNNNN would have been solved ? ) or is it just a matter of annotations ?

Answer: Short answer

Afaik, the sequences are not changed for the primary assembly, so if you don't go looking for the new sequences, you won't find them. However, some new corrections and new sequence are available if you want them. The annotations have been updated.

Long answer

NCBI 37.2, like Ensembl 60, is based on GRCh37.p2.

GRCh37.p2 is, as the name implies, just a patch release for GRCh37. It contains two kinds of patches, which should be seen as temporary updates until they are fully incorporated into the next major release of the genome.

A fix patch represents corrections to an existing sequence.

A novel patch represents novel sequence (perhaps filling some of the N's of runs).

As I understand it, these are not intended to replace the original sequences while the main release, GRCh37, is still in effect. They can be seen as a sneak preview of the next major release for those who are interested. This means that the GRCh37 primary assembly is the same between GRCh37, GRCh37.p1 and GRCh37.p2, and the patch sequences exist separately.

If you download the primary assembly sequence for Chromosome 5, say, from the GRC, you will get the original sequence, which doesn't include the updates. To get the updated sequences, you would need to download a separate patch file. Small annotation files are also provided that explain where the patch sequence "fits" in the original sequence.

Ensembl follows the pattern of keeping the original sequences separate from the patches - you would download the original sequences separately from the patch sequences.

I'm less sure of this, but I think NCBI also keeps the sequences separate. NCBI 37.2 contains a 'GRCh37.p2-Primary Assembly', and the patch sequences seem to be represented as a separate assembly - "GRCh37.p2-PATCHES".

Note also that GRCh37.p2 contains the mitochondrial sequence - previous versions of the assembly did not contain this.
Question: I have the sequenced data of an organism. but it has three 16srRNA which belong to 3 different organisms. I guess it could be contaminated. How could I extract the contigs belonging to each organism present in the sequence data?

Answer: If you truly feel that there are three organisms then you can use `bbsplit.sh` (from [BBMap suite][1]) to bin your reads into respective organismal pools. This will generally work well as long as the bacterial are distinct enough. You are able to decide what you want to do with reads that multi-map (map to all three reference genomes). e.g. keep in all bins, toss etc.

Use the answer here and ask if you have any questions: https://www.biostars.org/p/143019/#143040

Since you have bacterial data you could turn off `maxindel=0`.

[1]: https://sourceforge.net/projects/bbmap/
Question: Hi,

I want to extract read-pairs that aligned concordantly exactly 1 time to the reference genome. Is there any way to parse the output SAM file? I would really appreciate your suggestions.

Best regards,

Rahul

Answer: The simplest way is likely something along the lines of:

samtools view -hf 0x2 alignments.bam | grep -v "XS:i:" > filtered.alignments.sam

As I mentioned in the comment above, the -f 0x2 part will get only "properly paired" alignments, which will effectively be the concordant alignments. For the option #2 definition of "map only once", you can take advantage of the fact that bowtie2 will add an XS auxiliary tag to reads that have another "valid" mapping. So, a quick inverse grep (grep -v) can get rid of those. One possible problem with this is if only one read of a pair can map to multiple places (e.g., the original fragment partly overlapped a simple tandem repeat) then you'd end up with orphans. The easiest fix to that (if it's a problem) would be to just check if the read names of pairs of alignments are the same. I'm sure someone will come up with a way to do that in awk, but the following python script is likely simpler:

#!/usr/bin/env python
import csv
import sys

f = csv.reader(sys.stdin, dialect="excel-tab")
of = csv.writer(sys.stdout, dialect="excel-tab")
last_read = None
for line in f :
#take care of the header
if(line[0][0] == "@") :
of.writerow(line)
continue

if(last_read == None) :
last_read = line
else :
if(last_read[0] == line[0]) :
of.writerow(last_read)
of.writerow(line)
last_read = None
else :
last_read = line

If you saved that as foo.py and made it executable and in your PATH, then the following would solve the aforementioned issue:

samtools view -hf 0x2 alignments.bam | grep -v "XS:i:" | foo.py > filtered.alignments.sam

Note that this won't work with coordinate-sorted alignments, but bowtie2 doesn't produce those.
Question: Hi,

I'm starting out new in bioinformatics. I have couple of questions on searching dbSNP.

1. I'm searching (a list of rs# in [batch mode][1] in browser) dbSNP for SNP coordinates. dbSNP returns the coordinates in hg38 assembly build (I'm requesting a bed file for output format). I'd like to retrieve the coordinates in hg19 version. Is there a way to achieve this? dbSNP FAQ section doesn't mention if this could be done.

2. I would also like to know if it's possible to search genotype information for a given SNP (rs#). I would also like this in batch mode.

Any help is appreciated!

[1]: http://www.ncbi.nlm.nih.gov/projects/SNP/dbSNP.cgi?list=rslist

Answer: One way to do this is via the command line. You could download SNP annotations via `wget`. For example:

$ wget -qO- ftp://ftp.ncbi.nih.gov/snp/organisms/human_9606_b151_GRCh37p13/VCF/common_all_20180423.vcf.gz | gunzip -c | convert2bed --input=vcf --output=bed --sort-tmpdir=${PWD} - > hg19.snp151.bed

Filter via `grep` for the SNP of interest. For example, to search on a single SNP ID:

$ grep -F rs554008981 hg19.snp151.bed
1 13549 13550 rs554008981 . G A . RS=554008981;RSPOS=13550;dbSNPBuildID=142;SSR=0;SAO=0;VP=0x050000000005000026000100;GENEINFO=DDX11L1:100287102;WGT=1;VC=SNV;ASP;KGPhase3;CAF=0.9966,0.003395,.;COMMON=1;TOPMED=0.99221139143730886,0.00778064475025484,0.00000796381243628

To search on a file of IDs, e.g. a list of SNP IDs in `rsIDs.txt`:

$ grep -fF rsIDs.txt hg19.snp151.bed > matches.bed
Question: Heys,

I'm working with two ssh clusters and I want to transfer some data from one to another. I could download it to my pc and then to the other cluster, but is more than 100 GB and I imagined there would be an easy way to transfer the data between both clusters.

I've been trying to look for info in internet, but I didn't fully got it. If anybody could help me would be great!

Thanks in advance!

Answer: Anything wrong with just `rsync` it from A to B directly?

Something like this, assuming only the files to be transferred were in the current directory on the source cluster:

rsync --progress ./* username@servername:/path/to/destination