I'm trying to assemble a bivalve transcriptome. Since bivalves are filter feeders, their transcriptomes tend to be highly contaminated by bacteria, algae and whatnot. Since I pooled several transcriptomes, I have a high amount of reads (>2B reads). I also have a not-so-high quality reference genome.

Now, I'm not sure when to perform a BLAST-based filter on these reads:

  • Do I filter reads with BLAST before assembling them?
  • Do I assemble the transcriptome and then filter out contaminant contigs?
  • Since I have a reference genome, do I simply not care about contaminant reads because they will not align to it?

Having this amount of reads, I'd be inclined to first filter out those reads that don't have a match against Metazoa, and then to perform a second filter after assembly in order to retain only those contigs that match against Lophotrochozoa.

If you know resources on how to manage read contamination, feel free to share them.


First, to answer your question about mapping to a low-quality reference:

1. Mapping

For mapping, low genome contiguity (low N50) doesn't really matter. You will be using a spliced aligner and short reads, so even small portions of your reads will map if they match the reference genome. What does matter though is that you are confident that the draft genome assembly does not contain any scaffolds from contaminating species like bacteria, algae, et cetera.

2. Reasons to not BLAST your raw reads

There are a few reasons that I do not suggest starting with BLAST, but other options may better suit your needs.

  • BLASTing 2 billion reads will take several weeks to months even on a machine with many cores (90+).
  • Since you're not working with a model organism, there will be many reads that do not have a BLAST hit but still are probably derived from your bivalve. I do not suggest using BLAST to filter your reads for that reason.
  • Similarly, many contaminating reads will not hit to anything in NCBI databases, so you will miss many of the contaminants as well. Leaving them in for now and removing with the strategy I describe in #4, then reassembling, may work better.

3. Using multiple individuals of highly heterozygous species in assembly.

Most importantly for non-model organisms and highly heterozygous species - do not assemble a transcriptome using data from more than one individual. The number of allelic differences in highly heterozygous species will make the transcript graphs a jumbled mess and you'll end up with a bunch of partial and erroneous transcripts. You can merge transcriptomes from multiple individuals later with tools like orthoFinder.

4. A strategy to filter transcripts based on mapping and BLAST

Given #1, #2, and #3 above, I propose the following strategy to filter your reads that may prevent removing true bivalve reads and will guarantee removal of reads from known contaminating species.

I would try this:

  • Assemble the transcriptome of each individual separately.
    • This avoids misassemblies due to heterozygous sites and over-filtering by Trinity or other assemblers.
  • Concatenate all of the transcriptomes into one file. From this point on we will use the information from all individuals to help filter.
    • Use somthing like cat txome_individual1.fasta txome_individual2.fasta > all_individuals.fasta.
  • Map all the transcripts to the reference genome s using a spliced aligner like STAR or HISAT2.
  • Get a bam file of all the transcripts that didn't map to the reference genome. These may be from contaminants.
    • Use samtools view -f 4 -b mapped_transcripts.bam | samtools sort - > output.sorted.bam to only keep the transcripts that did not mapped to the reference. -f 4 is the magic flag here that keeps things that didn't map.
  • Make a list of all the transcripts that did not map to the reference...
    • samtools view output.sorted.bam | cut -f1 | sort | uniq > unmapped_transcripts.txt
  • ...and get only the transcripts from your fasta file that are potential contaminants.
    • seqtk subseq all_transcripts.fa unmapped_transcripts.txt > unmapped.fasta
  • BLAST all of the unmapped.fasta and identify which ones are from contaminating species. blobtools works well for this if there are too many to keep track of.
  • Use seqtk again to only get the known contaminating seqs from BLAST or blobtools. Let's say you call this known_contaminants.fasta
  • Map all reads from all individuals to the known_contaminants.fasta.
  • Get all the reads that map to known contaminating transcripts, known_contaminants.fasta using:
    • samtools view -F 4 all_reads_to_contams.bam | cut -f1 | sort | uniq > contaminating_reads.txt
  • Now remove all of the contaminating reads from your individual fastq files using bbmap filterbyname.sh.
  • Reassemble the filtered fastq files.

This strategy has the following benefits:

  • Does not require weeks/months of compute time to BLAST 2 billion reads.
  • Ensures that all reads from the reference genome are left untouched.
  • Ensures that all reads from known contaminants are removed.
  • Leaves in reads with no hits to BLAST, which could be from unassembled parts of the reference genome.
  • Avoids potentially filtering out low-abundance contaminants by not inputting all 2 billion reads to the transcriptome assembler, but assembling on a per-individual basis.
  • Leverages contaminants from different individuals to filter reads from all individuals.
  • $\begingroup$ Thanks for your answer, there are some interesting ideas. I want to use BLAST mainly to filter out those reads/transcript that do have a bacterial/algal hit - on my first assembly attempt, I ended up with >20kb long transcripts that were just bacterial genomes. Do allelic differences distort the graphs even with a genome-based assembly? If I assemble single individuals first, I'll have 35 individuals - not counting the pooled larvae transcriptomes, which in themselves contain 50k individuals. $\endgroup$ – LinuxBlanket Aug 22 '18 at 7:14
  • $\begingroup$ Allelic differences also distort the graph with genome-based assembly, causing many misassemblies and fragmenting the output. May I ask what sort of analysis that you're doing downstream with the transcriptomes? It will only take a small amount of time to make a shell script to assemble 35 transcriptomes compared to the amount of time it will take to discover all the incomplete and incorrect transcripts if all the reads are lumped into one assembly, IMHO. I have updated the answer to give more specific instructions and the reasoning behind each step. Hope it helps! $\endgroup$ – conchoecia Aug 23 '18 at 15:57
  • $\begingroup$ Wow, that's a very thorough answer! Thank you very much, it seems a very clever strategy! I want to use these transcriptomes for a variety of analyses, GSEA, 3'-UTR motifs and other. $\endgroup$ – LinuxBlanket Aug 23 '18 at 22:37
  • $\begingroup$ @conchoecia sorry for necro'ing this post but does point 3 also hold for de novo transcriptome assembly? (Because it appears your comments hint at this being relevant for genome-based assembly.) $\endgroup$ – Dunois Jun 29 '20 at 9:31
  • $\begingroup$ Hi @Dunois - I think I was being a bit dramatic above, but yes, point #3 is specifically relevant for de novo transcriptome assembly. You'll get more accurate results if a transcriptome is only assembled from reads from 1 individual. If you need different life stages or tissue types, my suggestion is to assemble the transcriptomes from different individuals independently, then remove duplicate transcripts using orthofinder. With that being said, if you have multiple libraries from different tissues from the same individual, then assemble all those together. $\endgroup$ – conchoecia Jun 30 '20 at 22:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.