Answers.Body: <p>First, to answer your question about mapping to a low-quality reference:</p>
  
  <h3>1. Mapping</h3>
  
  <p>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.</p>
  
  <h3>2. Reasons to not BLAST your raw reads</h3>
  
  <p>There are a few reasons that I do not suggest starting with BLAST, but other options may better suit your needs.</p>
  
  <ul>
  <li>BLASTing 2 billion reads will take several weeks to months even on a machine with many cores (90+).</li>
  <li>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.</li>
  <li>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 <strong>#4</strong>, then reassembling, may work better.</li>
  </ul>
  
  <h3>3. Using multiple individuals of highly heterozygous species in assembly.</h3>
  
  <p>Most importantly for non-model organisms and highly heterozygous species - <strong>do not assemble a transcriptome using data from more than one individual</strong>. 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 <a href="https://github.com/davidemms/OrthoFinder" rel="nofollow noreferrer">orthoFinder</a>.</p>
  
  <h3>4. A strategy to filter transcripts based on mapping and BLAST</h3>
  
  <p>Given <strong>#1</strong>, <strong>#2</strong>, and <strong>#3</strong> 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.</p>
  
  <p>I would try this:</p>
  
  <ul>
  <li>Assemble the transcriptome of each individual separately.
  
  <ul>
  <li>This avoids misassemblies due to heterozygous sites and over-filtering by Trinity or other assemblers.</li>
  </ul></li>
  <li>Concatenate all of the transcriptomes into one file. From this point on we will use the information from all individuals to help filter. 
  
  <ul>
  <li>Use somthing like <code>cat txome_individual1.fasta txome_individual2.fasta &gt; all_individuals.fasta</code>.</li>
  </ul></li>
  <li>Map all the transcripts to the reference genome s using a spliced aligner like <a href="https://github.com/alexdobin/STAR" rel="nofollow noreferrer">STAR</a> or <a href="https://ccb.jhu.edu/software/hisat2/index.shtml" rel="nofollow noreferrer">HISAT2</a>. </li>
  <li>Get a bam file of all the transcripts that didn't map to the reference genome. These may be from contaminants. 
  
  <ul>
  <li>Use <code>samtools view -f 4 -b mapped_transcripts.bam | samtools sort - &gt; output.sorted.bam</code> to only keep the transcripts that <strong>did not</strong> mapped to the reference. <code>-f 4</code> is the magic flag here that keeps things that didn't map.</li>
  </ul></li>
  <li>Make a list of all the transcripts that <strong>did not</strong> map to the reference...
  
  <ul>
  <li><code>samtools view output.sorted.bam | cut -f1 | sort | uniq &gt; unmapped_transcripts.txt</code></li>
  </ul></li>
  <li>...and get only the transcripts from your fasta file that are potential contaminants.
  
  <ul>
  <li><code>seqtk subseq all_transcripts.fa unmapped_transcripts.txt &gt; unmapped.fasta</code></li>
  </ul></li>
  <li>BLAST all of the <code>unmapped.fasta</code> and identify which ones are from contaminating species. <a href="https://blobtools.readme.io/docs" rel="nofollow noreferrer">blobtools</a> works well for this if there are too many to keep track of.</li>
  <li>Use seqtk again to only get the known contaminating seqs from BLAST or blobtools. Let's say you call this <code>known_contaminants.fasta</code></li>
  <li>Map <strong>all reads</strong> from <strong>all individuals</strong> to the <code>known_contaminants.fasta</code>.</li>
  <li>Get all the reads that map to known contaminating transcripts, <code>known_contaminants.fasta</code> using:
  
  <ul>
  <li><code>samtools view -F 4 all_reads_to_contams.bam | cut -f1 | sort | uniq &gt; contaminating_reads.txt</code></li>
  </ul></li>
  <li>Now remove all of the contaminating reads from your individual fastq files using <a href="https://www.biostars.org/p/127141/#127143" rel="nofollow noreferrer">bbmap filterbyname.sh</a>.</li>
  <li>Reassemble the filtered fastq files.</li>
  </ul>
  
  <p>This strategy has the following benefits:</p>
  
  <ul>
  <li>Does not require weeks/months of compute time to BLAST 2 billion reads.</li>
  <li>Ensures that all reads from the reference genome are left untouched.</li>
  <li>Ensures that all reads from known contaminants are removed.</li>
  <li>Leaves in reads with no hits to BLAST, which could be from unassembled parts of the reference genome.</li>
  <li>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.</li>
  <li>Leverages contaminants from different individuals to filter reads from all individuals.</li>
  </ul>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-08-22T07:14:54.890Z
Answers.Text: 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.
Answers.UserDisplayName: LinuxBlanket
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-08-23T15:57:01.927Z
Answers.Text: 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!
Answers.UserDisplayName: conchoecia
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-08-23T22:37:55.653Z
Answers.Text: 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.
Answers.UserDisplayName: LinuxBlanket
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-06-30T22:52:51.843Z
Answers.Text: 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.
Answers.UserDisplayName: conchoecia
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-08-22T00:22:19.983Z
Answers.LastActivityDate: 2018-08-23T15:54:56.220Z
Answers.LastEditDate: 2018-08-23T15:54:56.220Z
Answers.UserDisplayName: conchoecia
Body: <p>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.</p>
  
  <p>Now, I'm not sure when to perform a BLAST-based filter on these reads:</p>
  
  <ul>
  <li>Do I filter reads with BLAST before assembling them?</li>
  <li>Do I assemble the transcriptome and then filter out contaminant contigs?</li>
  <li>Since I have a reference genome, do I simply not care about contaminant reads because they will not align to it? </li>
  </ul>
  
  <p>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.</p>
  
  <p>If you know resources on how to manage read contamination, feel free to share them.</p>
  
Tags: <transcriptome><reference-genome><read-correction>
Title: At what stage of a transcriptome assembly is it better to perform read contaminant filter?
ThreadText: # At what stage of a transcriptome assembly is it better to perform read contaminant filter?
  
  (Asked by: username_0 on 2018-08-21)
  
  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.
  
  **Tags:** transcriptome, reference-genome, read-correction
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2018-08-22. Score: 3 | Accepted answer
  
  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](https://github.com/davidemms/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](https://github.com/alexdobin/STAR) or [HISAT2](https://ccb.jhu.edu/software/hisat2/index.shtml).
  * 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](https://blobtools.readme.io/docs) 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](https://www.biostars.org/p/127141/#127143).
  * 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.
  
  > **Comments:**
  > **username_0**: 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.  
  > **username_1**: 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!  
  > **username_0**: 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.  
  > **username_1**: 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.  
  
  ---
  
  