At what stage of a transcriptome assembly is it better to perform read contaminant filter?

Question

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.

Answer 1

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.