Background:
I am working on a de novo transcriptome assembly for nutmeg (Myristica fragrans) using Trinity. There is no reference genome available for this species, so I am assembling from scratch using paired-end Illumina RNA-seq reads. My samples come from two tissue types: fruit/flower (MFF) and leaf (MFL), with three biological replicates each (n=6 total).
As part of my pre-assembly quality control pipeline, I ran Kraken2 against a standard database (bacteria, fungi, human, viral) to detect and remove contamination. Kraken2 outputs two sets of files per sample: classified reads (matched to a known organism) and unclassified reads (nutmeg)
Problem:
My contamination summary looks like this:
| Sample | Total reads | Unclassified % | Classified % | Fungal % | Bacterial % |
|---|---|---|---|---|---|
| MFF1 | 31,917,229 | 20.16 | 79.84 | 78.09 | 0.77 |
| MFF2 | 26,281,714 | 20.59 | 79.41 | 77.68 | 0.77 |
| MFF3 | 37,982,254 | 20.10 | 79.90 | 78.26 | 0.76 |
| MFL1 | 26,827,296 | 92.23 | 7.77 | 1.35 | 3.43 |
| MFL2 | 31,593,301 | 92.33 | 7.67 | 1.30 | 3.36 |
| MFL3 | 25,838,245 | 92.24 | 7.76 | 1.31 | 3.43 |
The dominant contaminant in all MFF samples is Aspergillus chevalieri (~74% of total reads), with minor contributions from other Aspergillus species. This is likely biological in origin, mold growing on the fruit/flower tissue at time of collection or during storage since the signal is highly consistent across all three MFF replicates.
After Kraken2 filtering, my unclassified ("clean") read counts per sample are approximately:
- MFF samples: ~6–8 million reads each
- MFL samples: ~24–29 million reads each
Question
Decide between two approaches before running Trinity:
Option A — Filtered reads only (recommended by literature)
Run Trinity using only the Kraken2-unclassified _filtered.fq.gz files from all 6 samples. This removes the Aspergillus contamination but leaves the MFF samples with relatively low read counts (~6–8M each). My concern is whether these MFF reads still contribute meaningful tissue-specific transcript coverage, or whether they are too few to be useful.
Option B — Unfiltered reads (all 6 samples as-is) Run Trinity on the original, unfiltered reads. This preserves read depth but would almost certainly result in Aspergillus transcripts being assembled alongside nutmeg transcripts in the MFF samples, producing a mixed transcriptome. I understand this is generally not recommended, but I want to understand how bad the consequences would actually be in practice.
I would greatly appreciate any help or resources!
•
written
by
soundslikedaisy •
1 answer
I don't recommend filtering reads using Kraken2, because it makes mistakes. Instead, I suggest assembling using all reads and then aligning the transcripts by DIAMOND to NCBI nr and filtering transcripts based on results. Since NCBI nr is to some extent spoiled by unremoved contamination, I recommend looking not only at the single best match, but at the top 5 matches. If more than half of them belong not to Embryophyta, remove the transcript.
Log in to answer this question.
It looks like there are a number of entries available in GenBank even if there is no reference genome available. See --> https://www.ncbi.nlm.nih.gov/bioproject/?term=txid51089[Organism:noexp] At least one appears to be for a transcriptome project (is that your own data?) --> https://www.ncbi.nlm.nih.gov/bioproject/1359417
Do you expect to have contamination in samples that were prepared directly from plant tissues with GLP no doubt. Samples of leaf etc must have been surface sterilized before isolation of RNA from them? It is difficult to believe that majority of the data you have is from organisms other than the one of interest.
If all you want to do is assemble a first-step, gene catalog transcriptome, what shelkmike suggested would be the way to go. Just pool all the reads together (from all 6 samples), assemble them de novo, and then exclude transcripts that have more than some
x%coverage andy%identity across the topmmatches (you could just run a fast translated search usingMMseqs2; or usingDIAMOND blastxas was suggested).And/or, you could assemble a MFF transcriptome, a MFL transcriptome, and keep only all those transcripts in the MFF transcriptome that have matches (at some sufficient coverage and identity) in the MFL transcriptome.
Or, given that you've already done the pre-filtering for the assembly with
Kraken2, you might as well go ahead and assemble the filtered reads. Unless you're able to regenerate data, I doubt your results will improve. Perhaps worth looking into theKraken2parameters to see if you can increase the stringency to limit false positive classifications.Technically, you could even do all three approaches (and throw in some publicly available data, if any, a la, what GenoMax suggested), compare the results, and take some kind of an intersection of the transcripts for your "final" transcriptome.