I am new to bioinformatics and have a question about potentially using the Sequence Reach Archive for taxonomic classification to improve the amount of classified reads obtained from a sample using WGS. I am not sure if this is a good idea, or even possible since it is more of an archive after all, it seems most people download single sequences from the SRA to use in experiments.
Here is a general explanation of this class project I working on for an intro class to bioinformatics: I have some metagenomic data from one seabird gut microbiome sample that was sequenced using Illumina MiSeq (paired-end mode, length approx. 400-500 bp (300+300 bp with ~100-200 bp overlap), it was already analyzed using this workflow: https://github.com/LangilleLab/microbiome_helper/wiki/Metagenomics-Standard-Operating-Procedure-v3. According to this workflow, they use Kraken2+Bracken and the NCBI Complete RefSeq database. The amount of classified reads that came back from the analysis was around 5%, extremely low!
So for the class project I want to try and increase the number of classified reads for this sample. At the advice of my instructor, I was going to try and use the Sequence Read Archive as an alternative database to classify the reads taxonomically. Since there is a larger amount of data in that archive and will probably give a higher level of classification, is this even possible to do?
Note: the Github linked above they recommended using their VM, VirtualBox (Microbiome Helper VBox), to run small-scale metagenomic analyses on your laptop, so I started my analysis by running through the code provided on their Github. Now that I've come to the step where I have to classify these reads to a new reference database, I'm stuck.
Seeing as SRA is such a large database, I was trying to learn how to download their databases through AWS, but I would have to do this through their virtual machine to run their code even if it is possible. Since I am so new to this I am having trouble understanding how I can proceed now, if anyone has any advice I would really appreciate it. Thank you!
2 answers
Can you confirm that you are working with shotgun metagenomic data and not 16S rRNA amplicons? It seems strange that shotgun metagenomic data would be generated on an Illumina MiSeq instrument and have a 100-200 bp overlap -- this sounds like an amplicon library, which might require a different reference database (e.g., SILVA, GreenGenes, or RDP).
To answer your original question, building a database from the SRA is probably not feasible in this context, because you wouldn't have any information about the taxonomy of those reads, which Kraken2 requires for building a custom database. If you are using a MiniKraken database, then that might explain the low classification rate; if you are using the standard Kraken database, then this would be somewhat surprising.
However, it might be that birds have a majority of microbes that simply are not present in current reference database, hence the low classification rate (still surprising though). Have you checked the proportion of host (e.g., bird) vs. non-host (e.g., microbial) reads in your dataset? Maybe your sample is overwhelmed with bird DNA and has very little microbial DNA -- it seems unlikely for a gut sample, but might be worth checking out.
I would recommend revisiting the quality control, trimming and classification parameters before moving forward. I hope this helps!
Chris's answer is very good. A few notes from my side:
- The SRA is a READ archive, not a genome archive. Don't use it.
- You could try a few blasts to see if you can get anything else other that what you have already. I wrote a cluster blast for this using nextflow, which is ok for ~10-100k reads mostly
- try centrifuge/metaphlan3 as an alternative
- check your reads aren't all duplicates / rubbish with fastqc/multiqc
- include the host genome (even a bunch of combined eg cat'ted contigs if available) in any custom reference sequence you put in.
Hi, thanks for your help. The fastqc did not yield terrible results, it didn't dip below a phred score of 20 at all.
When I stitch the reads of the sample using PEAR it gives me a warning to check the summary log because less than 90% of my reads were assembled. Do you know why this could be?
ID Assembled Discarded Unassembled
Min 77.490 0.258 22.253
Mean 77.490 0.258 22.253
Max 77.490 0.258 22.253
COMU02-MN_S5 77.490 0.258 22.253
When I stitch the reads of the sample using PEAR it gives me a warning to check the summary log because less than 90% of my reads were assembled. Do you know why this could be?
This means your inserts as short and the paired end reads are overlapping in the middle when they sample the library fragment.
Sounds alright, I wouldn't be worried. But I think you need to answer key questions first
- is this MiSeq data (2x300bp can only be MiSeq, as Nextseq can only do 2x150bp)
- is this 16S amplicon data from the 16S gene only ..... or is it whole genome metagenomics (WGS) data ?
You need completely different pipelines for the two approaches.
For 16S I would recommend ampliseq https://nf-co.re/ampliseq/2.1.1/output
Log in to answer this question.
-
Low percentage of classified reads (~18%) after Kraken2 analysis – is this expected?
written by sumah.1si21bt017 •
-
Blastn and Kraken2 short reads clssification
written by Giulia.cosenza •
-
Finding FASTA files to use with bowtie-2-build to create a custom database for Kneaddata
written by jobie1 •
-
Composition of a microbiome
written by vikasmh111 •
-
What Kraken2 reference database is being used in these SOPs
written by jobie1 •
-
Reads that used for re-constructing MAGs can not be classified by Kraken2
written by Yifeng •
-
Ribosomal reads in shotgun metagenomics data
written by grp2009 •
-
Shotgun rarefactions - metagenomics (microbiome), MetaPhlAn2
written by robert.kwapich •
-
Deep metagenomics - further analyses
written by robert.kwapich •
-
Classify 16S Reads
written by sebabiokr •
