This is a test version of Biostars. For the public version, visit https://www.biostars.org.
How do I extract specific read sequences from an SRA run file, based on taxonomy?

I am running into a challenge and looking for some advice.

I am looking at the following file: https://trace.ncbi.nlm.nih.gov/Traces/?run=ERR11476772

If you look in the taxonomy, you'll notice that a small number of reads are identified as belonging to SIV and HIV (that's simian immunodeficiency virus and human immunodeficiency virus). I know that these could be false hits, but I still want to know: what are the sequences that were identified by SRA as HIV and SIV? Is there a way to do this in BigQuery or on my own computer?

  • What I have already tried:

I downloaded the file (prefetch and fastq-dump), trimmed the reads and discarded low quality ones (bbduk.sh) and then ran kraken2 with a custom library consisting of the standard library + all genome assemblies of SIVs. I found no reads were identified as part of the lentiviral genome. I tried processing the SRA file multiple times in different ways to try to get my results to match the SRA taxonomic viewer. I have looked at the differences between Kraken2 and SRA Tools, and I think the difference between what I'm seeing and the SRA Tools result might have to do with the fact that SRA Tools has a different database, and that Kraken2 requires exact kmer matching. I expect that any viral sequence in this run will differ pretty significantly in sequence from what is available on NCBI and Kraken2 databases.

  • What I want to do:

Please let me know if there is any way for me to get these sequences out of SRA. Would I have to pay for BigQuery or is there a free way to do it?

virus taxonomy illumina ncbi sequencing
• 641 views • link
updated by GenoMax GenoMax • written by ggio • ggio

1 answer

I know that these could be false hits, but I still want to know:

You should be able to use bbduk.sh in filter mode to look for these sequences. Try something like this following (adjust references as needed)

$ bbduk.sh -Xmx4g in1=ERR11476772_1.fastq in2=ERR11476772_2.fastq ref=GCF_000863925.1_ViralProj15501_genomic.fna,GCA_003190765.1_ASM319076v1_genomic.fna k=24 outm1=match_R1.fastq.gz outm2=match_R2.fastq.gz

Add parameter hdist=N (replace a number) to allow for 1 or more errors (mismatches in reads).

Alternatively you could use the host genome (looks like it is yellow baboon?) and find reads that don't match the host genome (which should be all remaining reads). Something like following

$ bbduk.sh -Xmx4g in1=ERR11476772_1.fastq in2=ERR11476772_2.fastq ref=baboon_genome.fa k=24 outu1=ummatch_R1.fastq.gz outu2=unmatch_R2.fastq.gz
• 0 views • link
written by GenoMax GenoMax

Thank you so much! I just ran your code and got a nice set of sequences. I have a couple of questions if you don’t mind: (1) Is there a difference between what the SRA Toolkit and bbduk.sh would consider a “match”? For example, I see you used k=24 in the bbduk command; should I use k=32 if I want to match the SRA Toolkit (which uses a kmer length of 32)? (2) I got about 5,000 reads in total from your code. What would you recommend as a next step -- should I align these reads to the SIV/HIV reference genomes? Thanks so much again -- I've been scratching my head about this for a while.

• 0 views • link
written by ggio • ggio

The k value in bbduk is for initial matching of read fragments, I arbitrarily selected 24. You can use 32 if you want. Since the reads are 150 bp that would be well below half the value.

As for the next step, you can align the reads to SIV/HIV genomes (if you don't expect both viruses to be present at the same time, align them independently) to see what you get. bbmap.sh the aligner would be perfect/fast for these alignments. See https://bbmap.org/tools/bbmap for parameter explanations.

Something like following should work (make sure you have samtools or sambamba available, otherwise you will get SAM files)

$ bbmap.sh -Xmx20g in1=R1.fq.gz in2=R2.fq.gz ref=HIV.fa out=aligned.bam slow k=12
• 0 views • link
written by GenoMax GenoMax

Log in to answer this question.

More posts like this