cannin/biostars_qa huggingface dataset:
Question: I want to extract trinity gene id only from below text file: like

TRINITY_DN33489_c0_g1_i1
TRINITY_DN33489_c0_g2_i1
TRINITY_DN33447_c0_g1_i1

# Query: TRINITY_DN33489_c0_g1_i1 len=657 path=[818:0-102 898:103-168 964:169-182 978:183-192 988:193-217 1013:218-224 1968:225-248 1044:249-249 1045:250-265 1061:266-273 1069:274-289 1085:290-329 1965:330-353 1149:354-439 1969:440-463 1259:464-656] [-1, 818, 898, 964, 978, 988, 1013, 1968, 1044, 1045, 1061, 1069, 1085, 1965, 1149, 1969, 1259, -2]

# Query: TRINITY_DN33489_c0_g2_i1 len=816 path=[261:0-148 387:149-278 1963:279-302 541:303-433 672:434-434 1964:435-458 697:459-567 806:568-579 25:580-591 37:592-598 44:599-601 47:602-612 58:613-619 65:620-622 68:623-636 82:637-643 89:644-644 1966:645-668 114:669-749 1967:750-773 219:774-815] [-1, 261, 387, 1963, 541, 672, 1964, 697, 806, 25, 37, 44, 47, 58, 65, 68, 82, 89, 1966, 114, 1967, 219, -2]

# Query: TRINITY_DN33447_c0_g1_i1 len=566 path=[1:0-68 47:69-90 69:91-92 71:93-114 93:115-174 807:175-198 177:199-207 186:208-225 204:226-231 210:232-249 228:250-266 808:267-290 269:291-339 806:340-363 342:364-479 458:480-480 459:481-483 462:484-502 481:503-504 483:505-507 486:508-526 505:527-565] [-1, 1, 47, 69, 71, 93, 807, 177, 186, 204, 210, 228, 808, 269, 806, 342, 458, 459, 462, 481, 483, 486, 505, -2]

Answer: That was not clear from your initial question. You also didn't show how the data looks like for ids without hit.

What about:

cat yourfile.txt | tr ' ' '\n' | grep '^TRINITY' > output.txt
Question: I am trying to use a pipeline that uses formatdb command of old BLAST. This option is been replaced by makeblastdb in newer versions of blast+ . I could not find archive for blast to download the old version. I am using CentOS 6.5.

Here is the program, if any body is interested to look at the source code:

https://github.com/chiulab/surpi/blob/master/plot_reads_to_gi.sh

Answer: The term you're looking for is "Legacy BLAST" :-)

Here we go: ftp://ftp.ncbi.nlm.nih.gov/blast/executables/release/LATEST/

Also, your source code can be upgraded by modifying

formatdb -p F -i $3.$2.fasta

to

makeblastdb -dbtype nucl -in $3.$2.fasta

and

blastall -p blastn -m 8 -a 1 -b 1 -K 1 -d $3.$2.fasta -i $f -o $f.$2.$3.$4.blastn -e $4 >& $f.$2.$3.$4.error

to

legacy_blast.pl blastall -p blastn -m 8 -a 1 -b 1 -K 1 -d $3.$2.fasta -i $f -o $f.$2.$3.$4.blastn -e $4 >& $f.$2.$3.$4.error --path /path/to/blast/bin

EDIT: I ***highly*** recommend updating your copy of the code to use blast+. In fact, you could do the modifications and submit a PR, and the repo owner will appreciate it.
Question: I am new to NGS data analysis and I'm working in a multiple-sample variant calling workflow. I have Illumina-Miseq fastq files (paired end, raw reads) for a father, mother and child trio, one pair for each individual, totalling 6 files. I could trim, align, do the pre-processing and variant calling for each individual pair separately (I'm skipping indel-realignment and quality recalibration, for the sake of simplicity, as this workflow is intended for learning only), but I wish to merge the samples into a single file. I wish that the alignment step (with BWA-MEN), the pre-processing steps (with Picard) and the variant calling step (with FreeBayes), are done at once for all samples, if possible and correct, while taking in consideration the correct paired end mates and the respective read groups (when applicable).

My final goal is to obtain a single vcf file from which I'll compute the total number of different kinds of variants.

At which step, in which file format and with which Galaxy tools can I merge the samples in a manner that I can get correct, meaninful results at the variant calling step?

Answer: avoid to merge the VCFS : call them at the same time:

samtools mpileup [options...] father.bam mother.bam child.bam | bcftools call [options...] > result.vcf

or something like:

java -jar GATK.jar -T HaplotypCaller [options...] -o result.vcf -i father.bam -i mother.bam child.bam

or you can use GATK + gVCF calling : http://gatkforums.broadinstitute.org/gatk/discussion/4017/what-is-a-gvcf-and-how-is-it-different-from-a-regular-vcf

or you can merge the 3 VCF later with https://software.broadinstitute.org/gatk/documentation/tooldocs/3.5-0/org_broadinstitute_gatk_tools_walkers_variantutils_CombineVariants.php , but you'll get some missing genotypes (./.) where you can't tell if the sample was HOM_REF or if there was not enough information to do the calling.
Question: I'm reviewing the UniProt data model and I'm confused by what is a protein. Originally, I thought that a protein was defined by its sequence of amino acids. However, the fact that isoforms of a gene are stored in one entry for UniProt/Swiss-Prot leads me to believe a protein is defined by the gene from which it originates. Otherwise, these alternative splicing would receive distinct UniProt/Swiss-Prot entries. Or perhaps it is more complicated than that and if isoforms are distinctly different enough they receive different entries? Just confused a bit by the definition of protein in this light.

Looking at the word itself, iso- means equal and isoform would seem to mean equal form. Understanding that structure is more conserved than sequences I would understand that form is better to define a protein than sequence. But I don't think this is what isoform really means from looking at the data.

Any help appreciated.

Answer: https://ghr.nlm.nih.gov/primer/howgeneswork/protein

Protein is a molecule made of amino acids that have a specific function. Sometimes changes in protein sequence due to alternative splicing are very little and have almost no effect on function in other cases differences in sequence and function are dramatic. So one gene can encode multiple isoforms of the same protein and/or multiple different proteins.
Question: Dear All,

I’m making a command line for HISAT following the manual instruction.
What is not clear to me is how to use the genome index files.
I downloaded from HISAT website the reference genome (genome.1.ht2… until genome.8.ht2)

My question is how to use this? the manual says “any of the index files up to but not including the final .1.ht2”
can I just pick one, insert into the line command and remove thee others from the folder?

My command would be like

hisat2 -q -x genome -1 mymate1.fq -2 mymate2.fq -S -o resultpath

I know it's not correct.. but can you kindly help me to fix it?

Thank you

Answer: what the manual is trying to explain to you with

> “any of the index files up to but not including the final .1.ht2”

is that you need to provide the file name you used for the index but omit the '`.1.ht2`' part from it. So using simply 'genome' as you did in your example is correct (given the index names you provided).

and, no, you can not remove the others from your directory, the software needs all of them to work.

----------

EDIT: to resolve all issues with the posted cmdline, you will need to add a correct value (== a filename to write the sam output to) for the `-S` parameter.
Question: Do freeBayes and platypus support gVCF?

https://sites.google.com/site/gvcftools/home/about-gvcf

With GATK you used to have the 1000 + 1 genome variant calling problem, where you had to redo the whole expensive multi-sample variant calling (gatk haplotype caller) if you wanted to add a single genome.

Now GATK Haplotype caller can create intermediate gVCF files for each sample and do a much less expensive merging of the gVCF files to a multi-sample vcf.

https://www.broadinstitute.org/gatk/gatkdocs/org_broadinstitute_gatk_tools_walkers_variantutils_GenotypeGVCFs.php

Do freeBayes and and platypus also support this or are they planning to support this?

Answer: Basic support for gVCF has been added to Freebayes just this weekend. More support for gVCF will be added in the future.....

There are some inherent problems with gVCF and complex multi-allelic variant sites. For just snps it should work:

See also this issue at the Freebayes github site, the last comment is the update from Erik stating the (upcoming) gVCF support.

https://github.com/ekg/freebayes/issues/76
Question: Hello,

In choosing a Next-gen seq platform, given the choices are HiSeq and NextSeq, which one would be more benefical? I know that HiSeq is very high throughput, compared to Nextseq but Nextseq produces sequences faster. Considering our sample size would be anywhere between 1500-3000 and mostly our experiments would be WGS, WES and a bit of functional work as well (methylation, RNA-seq), does anyone have any views from a data-analysis perspective as to pros and cons of the 2 systems ?

Many thanks

Answer: There have been a couple of reports about data quality for PCR sample preps vs PCRFree preps as well as reports of the NextSeq platform on seqanswers comparing it to HiSeq V3. There is also another report (Slideshare) comparing HiSeq V4 to HiSeq V3.

Overall, the trend in those reports is that the data quality is, in order of best to worst:

Illumina PCRFree > Illumina Nano PCR >> Other PCR sample preps

Sample prep is the biggest factor for data quality at the moment, then followed by instrument differences:

HiseqV4 PCRFree > HiseqV3 PCRFree > NextSeqV1 PCRFree

There is one sweet-spot for the NextSeq which was reported by Dale Yuzuki, who works for Thermo (Ion Torrent), which is the low price per 100M reads for 1x75bp runs in the NextSeq, which are commonly used in RNA tag counting experiments:

http://www.yuzuki.org/clarifications-ion-torrent-pii-nextseq-500/

In terms of data management and analysis, the NextSeq is highly integrated with Illumina BaseSpace. From the design of the experiment, barcoding of the samples, uploading of the resulting read sets up to the cloud, and then running of analysis in BaseSpace apps and sharing the results with collaborators.

The HiSeq machines are gaining data integration, but it is not as seamless as the NextSeq at the moment (2014-11-28).

If you don't want to upload your datasets on the BaseSpace private cloud, you can use BaseSpace OnSite, which is reportedly as good for the NextSeq as it is for the HiSeq.

In terms of throughput, find the relevant information below:

A run on the HiSeq with V4 chemistry 2x125bp takes 6 days and produces between 0.9-1Tbp of data:

http://systems.illumina.com/systems/hiseq_2500_1500/performance_specifications.html

At 15-20x per lane, and 16 lanes per run, that's 8 WGS 30-40x 2-lane genomes every 6 days, and 8*50 = 400 WGS per year not running it on Sundays and 2 weeks holiday period.

For exomes, multiply the 400 WGS by 12, and it gives you 4800 WES per year.

If you do half and half WGS and WES, you can do 200 WGS and 2400 WES per year with a HiSeq using V4 chemistry.

For RNA profiling, multiply 400 by 24, and it gives you 9600 RNAseq per year.

If you can do 200 WGS, 1200 WES and 2400 RNAseq per year with a HiSeq using V4 chemistry.

http://systems.illumina.com/systems/sequencing-platform-comparison.html
Question: Hi,

I have been using the package [fgsea][1], and I would like to use a for loop in R to output multiple Enrichment plots.

So far I can generate a plot of the desired pathway using:

plotEnrichment(pathway =
pathways.hallmark[["HALLMARK_EPITHELIAL_MESENCHYMAL_TRANSITION"]],
gseaParam = 1, ticksSize = 0.5, stats= ranks)
+ labs(title="Epithelial to Mesenchymal Transition")
+ theme(plot.title = element_text(hjust = 0.5, face="bold"))

I decided to capture the statistically relevant pathways, and attempted to use a for loop to output them all at once:

filtered_pathway <- subset(fgseaResTidy, pval < 0.05)

filt_p <- as.vector(filtered_pathway$pathway)

filt_p
[1] "HALLMARK_APICAL_JUNCTION" "HALLMARK_MYOGENESIS"
[3] "HALLMARK_HEDGEHOG_SIGNALING" "HALLMARK_PI3K_AKT_MTOR_SIGNALING"
...

for (i in filt_p){
plotEnrichment(pathway = pathways.hallmark[[i]],
gseaParam = 1, ticksSize = 0.5, stats= ranks) +
labs(title=i) + theme(plot.title = element_text(hjust = 0.5, face="bold"))
}

This code does not generate anything. I am aware this is not a bioinformatics question, however I have no doubt many would be able to point out the problem in 2 seconds. Hopefully this post will be helpful to people using fgsea in the future.

Kind regards,

Barry

[1]: http://bioconductor.org/packages/release/bioc/html/fgsea.html

Answer: The loop you are running will generate plot but will not save anything .
I have modified the code, I hope it helps.

for (i in filt_p){
pdf(paste0(i,".pdf"),height=5,width=5) # change height and width parameter
plt <- plotEnrichment(pathway = pathways.hallmark[[i]],
gseaParam = 1, ticksSize = 0.5, stats= ranks) +
labs(title=i) + theme(plot.title = element_text(hjust = 0.5, face="bold"))
print(plt)
dev.off()
}
Question: Hello everyone,

I have previously worked with gene tree/species tree evolution tools but not tumor evolution so if it looks like a silly question, I beg your pardon.

I have variant calling data (using Mutect and VarScan) for whole-exome sequencing data. In the review paper, such as http://sysbio.oxfordjournals.org/content/64/1/e1, the authors say that, for instance, for PyClone, Phylosub and SciClone, the input data is single nucleotide variant (SNV). When I check, for instance, in test data in PyClone, it has a nicely formatted input file which is used to generate the output.

I am wondering if someone has worked with any pipeline for phylogenetic tree reconstruction using "any" of tumor evolution software? I would be grateful if you can share your experience (and, if possible, the steps).

Thanks in advance

Answer: Once you've clustered your data to identify subclonal populations (with the algorithms you list above), there are a few different tools that will attempt to construct evolutionary trees. The ones I'm aware of are below:

ClonEvol: https://github.com/hdng/clonevol

SubcloneSeeker: http://www.genomebiology.com/2014/15/8/443

Schism: http://karchinlab.org/apps/appSchism.html
Question: Hi there,

I have a rather large VCF file annotated using VEP, and I am trying to pick the most severe consequence of the variant for each variant in the file. Now I know that VEP automatically arranges the output to be listed as most-severe to least according to this (http://uswest.ensembl.org/info/genome/variation/prediction/predicted_data.html#consequences). However, my VCF doesn't have it in that format (no idea why, but I spot checked, and it doesn't) .

I cannot re-make the VCF using --severe because I didn't create it in the first place. My solution is to write a script to do this manually. But before I do that, I'd be interested if anyone out there has encountered a similar situation and have better suggestions for this.

Thanks for the help!

Answer: Hi shwethacm,

I'm afraid that I'm not able to help with any particular suggestions or custom scripts that will help to filter or rearrange your VCF output. However, I wanted to correct a small error in your assumption- the VEP does not necessarily arrange the output from most to least severe for each variant.

Some considerations for your filtering and reordering:

- You may be able to run your VCF through the VEP again using the different --pick options:
http://www.ensembl.org/info/docs/tools/vep/script/vep_options.html#filt
- You may be able to use the filter_vep script to filter your results:
http://www.ensembl.org/info/docs/tools/vep/script/vep_filter.html

Best wishes

Ben
Ensembl Helpdesk
Question: Hi,

Could you please suggest to me how I can get all these 82697 sequences from this website using the linux command:

http://fungene.cme.msu.edu/hmm_detail.spr?hmm_id=721&page=1

Thanks

Answer: (edited) To download the nucleotide sequences into 4 parts simultaneously:

```
API_KEY="REPLACE_WITH_YOUR_NCBI_EUTILS_API_KEY"

for p in {1..83}; do
curl -s "http://fungene.cme.msu.edu/hmm_detail.spr?hmm_id=721&page=${p}" \
| grep "gbnucdata" \
| sed -r "s|.+seqAccno=(.+)\&format=Genbank.+|http://www.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nucleotide\&retmode=text\&rettype=fasta\&id=\1\&api_key=${API_KEY}|" \
>> nucl_urls.txt
done

split -d -l 21000 nucl_urls.txt part
ls part0* | parallel -j 4 wget -O {}.fasta -i {}
```

---

(original answer to download protein sequences)
You can download the sequences in this way (here for example all the protein sequences):

```
# First get all the accession numbers of protein sequences in all the pages
for p in {1..83}; do
curl -s "http://fungene.cme.msu.edu/hmm_detail.spr?hmm_id=721&page=${p}" \
| grep "gpprotdata.jsp?seqAccno" \
| sed -r 's|.+>(.+)<\/a><\/td>|\1|' \
>> prot_acc.txt
done

for acc in $(cat prot_acc.txt); do
wget -O ${acc}.fasta "http://www.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=protein&retmode=text&rettype=fasta&id=${acc}"
done
```
Question: Hi,

Can anybody please tell me how to install samtools v.1.2 on Ubuntu 14.0 32 bit?

I am very new to Linux environment and NGS tools. I have tried the following but the version is not what I want.

sudo apt-get install samtools

Please help me with the appropriate link to download samtools and the commands to install it.

Many thanks

Answer: wget "https://github.com/samtools/samtools/releases/download/1.2/samtools-1.2.tar.bz2" && tar xfj samtools-1.2.tar.bz2 && (cd samtools-1.2 && make)

then put samtools in your path: https://linuxconfig.org/permanently-add-a-directory-to-shell-path

export PATH=/full/path/to/samtools-1.2:${PATH}
Question: Tree #1 from https://www.ncbi.nlm.nih.gov/pubmed/26773003 :

![enter image description here][1]

Tree #2 from http://www.fisht1k.org/subproject/phylogeny :

![enter image description here][2]

I would like to group sequences on the tree and have a table/graphics on the side. Also, it would be great if at least some part of this task can be done programmatically. I was wondering what could be a good tool (or an R library, for example) to do that?

Thank you in advance,
Ivan

[1]: http://ivanya.com/tmp/2016/Tree1.png
[2]: http://ivanya.com/tmp/2016/Tree2.png

Answer: I'm not sure if it will fit your needs exactly, but you can do many many things using ggtree in R http://bioconductor.org/packages/release/bioc/html/ggtree.html
Question: Is it enough to find a transgene sequence in the genetically engineered bacteria (via BLAST server) to confirm that the engineering was successful? In general, how to prove that genetic engineering was successful by using bioinformatics tools? Is the phylogenetic tree would be recommended? What else could you suggest?

Answer: If you are engineering a bacteria, you are probably working with a single-isolate culture. Then, the complete workflow would be: sequence the original strain, sequence the putative transformed strain - best to sequence both with PacBio or Nanopore, assemble both genomes. Then there is a number of analyses you should perform: align the genomes, search the transgenes in both genomes (with Blast or other program), check if the transgene is complete, if the marker gene (if any, e.g., some antibiotic resistance gene) is present, and so on.

If you are engineering an eukaryote, you can follow the same workflow, but the more complex the genome (e.g., some polyploid plants), the more difficult will be to fully characterize the transgene (where it inserted, how many copies, if it integral or not, etc). In this case, wet lab techniques can be faster and cheaper than sequencing plus bioinformatics.
Question: I use picard software to mark duplicates, and here is my command:

```
java -d64 -server -XX:+UseParallelGC -XX:ParallelGCThreads=2 -Xms8g -Xmx16g -Djava.io.tmpdir=tmp -jar ./picard.jar MarkDuplicates I=input.bam O=out_markdup.bam METRICS_FILE=out.metrics ASO=coordinate VALIDATION_STRINGENCY=LENIENT
```

It works well but when the `input.bam` file gets bigger, the speed is very slow! I found that the picard MarkDuplicates doesn't support multiple threads. So, is there anyway to speedup picard? Another way, is there any better software to do the same as picard MarkDuplicates but with less time? I know [elprep](https://github.com/ExaScience/elprep) is another choice, but it needs very large memory!

Besides, I found that samtools can also remove duplicates, but according to my search, samtools can not remove the duplicates cross different chromosomes, so picard is more universe.

Any reply will be much appreciated!

Answer: I have this toolkit which I am currently developing and is therefore not to be "browsed too much", where the *deduplicate* module is completed for the paired-end reads. I would be interested in a test from another user, if you want to give it a try!
https://github.com/MatteoSchiavinato/SPIQR

The principle is that it removes PCR duplicates by sequence identity and not by mapping position, so it will speed us the mapping step as well and not only the deduplication. We tried with our data and the intersection of reads deduplicated by this tool and by picard was ~99%.
Question: This might be easy but I have not been able to solve it - I have a normal fasta file as below:

```
>seq1
ATCTAC
>seq2
AATCGCATCG
>seq3
ATATACAGC
>seq1
ATCGCGGGC
>seq4
ATTAATTTTAT
```

what I want o do is to reverse the order of lines in the file. So the desired output should be:

```
>seq4
ATTAATTTTAT
>seq1
ATCGCGGGC
>seq3
ATATACAGC
>seq2
AATCGCATCG
>seq1
ATCTAC
```

I tried using ```tail -r file.fa``` but then it does not work because ids get placed after the sequences.

How can I do this?

Thank you

Answer: Another way this can be done is to read the FASTA file from its end to its start.

Reading the file backwards avoids having to use hash tables or sorting, which can be expensive on very large inputs. This approach will be fast and use very little memory:

#!/usr/bin/env python

import sys
import os

ifn = sys.argv[1] if len(sys.argv) == 2 else None

def reversed_lines(file):
"Generate the lines of file in reverse order."
part = ''
for block in reversed_blocks(file):
for c in reversed(block):
if c == '\n' and part:
yield part[::-1]
part = ''
part += c
if part: yield part[::-1]

def reversed_blocks(file, blocksize=4096):
"Generate blocks of file's contents in reverse order."
file.seek(0, os.SEEK_END)
here = file.tell()
while 0 < here:
delta = min(blocksize, here)
here -= delta
file.seek(here, os.SEEK_SET)
yield file.read(delta)

def main():
ifh = sys.stdin if not ifn else open(ifn, "r")
lines = reversed_lines(ifh)
sequence = ""
header = ""
for line in lines:
if line.startswith('>'):
header = line.rstrip()
sys.stdout.write('{}\n{}\n'.format(header, sequence))
sequence = ""
else:
sequence = line.rstrip() + sequence
ifh.close()

if __name__ == "__main__":
main()

Usage:

$ ./reverse_fasta_records.py < in.fa > out.fa

Or:

$ ./reverse_fasta_records.py in.fa > out.fa

Edit: Modified this script to support the input file either coming from standard input or a regular file.
Question: Where can I do an online project to get some hands on experience with NGS data analysis?

Answer: For gaining experience on NGS data analysis, you don't get any online projects. You need to join any lab doing the kind of work you would like to do or you need to learn on your own. If the former is not possible go through review papers on NGS (decide which one WES/WGS/RNA-seq...etc by reading wiki of each) then download the raw data from NCBI-SRA and try to reproduce the results. It will give you some idea how things work in the process. Previous posts that are helpful

- Register for Genomics related courses on Coursera like this
- https://www.biostars.org/p/80739/
- https://www.biostars.org/p/142494/
- https://www.biostars.org/p/76171/
- Introduction to AWS cloud computing - This entire wikipage on RNA-seq (Informatics of RNA-seq) would be *very \* Inf* helpful.

And many other, if you spend time searching.