Answers.Body: <p>I would recommend <a href="http://hibberdlab.com/transrate/index.html" rel="nofollow noreferrer">Transrate</a>. It generates a score for each transcript based on many metrics (read their paper <a href="http://dx.doi.org/10.1101/gr.196469.115" rel="nofollow noreferrer">here</a>) and, using that score, gives a predicted cut-off to produce the optimal transcript set.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-13T08:35:41.727Z
Answers.LastActivityDate: 2017-07-13T08:44:04.673Z
Answers.LastEditDate: 2017-07-13T08:44:04.673Z
Answers.UserDisplayName: ithinkiam
Answers.Body: <p>I used a combination of BUSCO and Salmon to filter transcripts based on their abundance in the RNASeq read dataset. The approach was roughly as follows:</p>
  
  <ol>
  <li>Run BUSCO in short/transcript mode on the Trinity-generated sequences</li>
  <li>Run Salmon on the Trinity-generated sequences, using the RNASeq reads that were used to generate the transcriptome</li>
  <li>Merge the BUSCO full results table with the Salmon results table to determine distribution of the number of reads mapped for BUSCO genes</li>
  <li>Choose an appropriate "real" read count threshold based on the read count distribution for BUSCO genes</li>
  <li>Sub-select the transcripts based on this threshold, using the Salmon results table</li>
  </ol>
  
  <h2>Step 1 -- BUSCO</h2>
  
  <pre><code>python ~/install/busco/BUSCO.py -i ../Trinity-GG.fasta -o BUSCO_Trinity-GG_nematodes -l ~/install/busco/nematoda_odb9 -m tran -c 10
  </code></pre>
  
  <h2>Step 2 -- Salmon</h2>
  
  <pre><code>~/install/salmon/Salmon-0.8.2_linux_x86_64/bin/salmon index -t Trinity-GG.fasta -i Trinity-GG.fasta.sai
  ~/install/salmon/Salmon-0.8.2_linux_x86_64/bin/salmon quant -i Trinity-GG.fasta.sai -1 reads-all_R1_trimmed.fastq.gz -2 reads-all_R2_trimmed.fastq.gz -p 10 -o quant/reads-all_quant -l A
  </code></pre>
  
  <h2>Step 3 -- Results merge (using R)</h2>
  
  <pre><code>#!/usr/bin/Rscript
  busco.table &lt;-
      read.delim("BUSCO/run_BUSCO_Trinity-GG_nematodes/full_table_BUSCO_T-BNOCFED_nematodes.tsv", comment.char="#",
                 header=FALSE, stringsAsFactors = FALSE,
                 col.names=c("Busco id","Status","Sequence","Score","Length"));
  busco.table &lt;- subset(busco.table, Status != "Missing");
  
  count.table &lt;- read.delim("quant/reads-all_quant_ISR/quant.sf",
                            stringsAsFactors=FALSE);
  
  busco.expr.df &lt;- merge(busco.table, count.table,
                         by.x="Sequence", by.y="Name", all.x=TRUE);
  </code></pre>
  
  <h2>Step 4 -- Count thresholding (visual exploration)</h2>
  
  <pre><code>pdf("busco_NumReads.pdf", paper="a4r", width=11, height=8);
  options(scipen=10);
  busco.numReads &lt;- sort(tapply(busco.expr.df$NumReads, busco.expr.df$Busco.id,
                                function(x){
                                    exp(mean(log(x[x&gt;0])))}), decreasing=TRUE);
  par(mfrow=c(1,2));
  plot(x=seq(0,1,length.out=length(busco.numReads)),
       busco.numReads, log="y", ylab="Number of Mapped Reads",
       yaxt="n",
       xlab="Proportion of Complete BUSCO sequences");
  axis(2, at=rep(c(1,2,5),each=6) * 10^(0:5), las=2, cex.axis=0.71);
  abline(col="#00000020", h=rep(1:9,each=6) * 10^(0:5), lwd=0.1);
  abline(col="#80808080", h=rep(1,each=6) * 10^(0:5), lwd=2);
  rect(ybottom=10^(par("usr")[3]), ytop=10^(par("usr")[4]),
       xleft=0.9, xright=par("usr")[2],
       col="#80808020", border=NA);
  plot(x=seq(0,1,length.out=length(busco.numReads)), xlim=c(0.9,1),
       busco.numReads, log="y", ylab="Number of Mapped Reads", yaxt="n",
       xlab="Proportion of Complete BUSCO sequences");
  axis(2, at=rep(c(1,2,5),each=6) * 10^(0:5), las=2, cex.axis=0.71);
  abline(col="#00000020", h=rep(1:9,each=6) * 10^(0:5), lwd=0.1);
  abline(col="#80808080", h=rep(1,each=6) * 10^(0:5), lwd=2);
  invisible(dev.off());
  </code></pre>
  
  <p><a href="https://i.sstatic.net/7iw47.png" rel="noreferrer"><img src="https://i.sstatic.net/7iw47.png" alt="choosing read count threshold from BUSCO genes"></a></p>
  
  <p>Looks like a value of 50 will be suitable.</p>
  
  <h2>Step 5 -- Sub-select the transcripts</h2>
  
  <pre><code>cat(file="HighCount_Transcripts_Num50.txt",
      subset(count.table, NumReads &gt;= 50)$Name, sep="\n");
  </code></pre>
  
  <p>Then filter using a suitable <a href="https://github.com/gringer/bioinfscripts/blob/master/fastx-fetch.pl" rel="noreferrer">FASTA selection utility</a>:</p>
  
  <pre><code>pv Trinity-GG.fasta | ~/scripts/fastx-fetch.pl -i HighCount_Transcripts_Num50.txt &gt; HighCount50_Trinity-GG.fasta
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-12T09:39:48.067Z
Answers.Text: I really like the apporach! How do you handle/treat replicates?
Answers.UserDisplayName: Kristoffer Vitting-Seerup
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-12T11:12:42.807Z
Answers.Text: Do you mean identical transcripts (or sub-transcripts)? I've used a very strict cd-hit-est clustering to clean up a few (98% identity), but the BUSCO results still indicate that there's a lot of duplication going on. I'm hesitant to use a looser identity threshold because I don't want to include any of the real gene duplications.
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-07-12T08:41:00.213Z
Answers.LastActivityDate: 2017-07-12T08:54:27.503Z
Answers.LastEditDate: 2017-07-12T08:54:27.503Z
Answers.UserDisplayName: gringer
Answers.Body: <p><a href="http://deweylab.biostat.wisc.edu/detonate/" rel="nofollow noreferrer">DETONATE</a> offers a pipeline for filtering in this way as well. But the aligner used by default is Bowtie2...</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-12-20T16:18:01.330Z
Answers.LastActivityDate: 2018-12-20T16:18:01.330Z
Answers.UserDisplayName: winni2k
Body: <p>I have recently generated a genome-guided transcriptome with <a href="https://github.com/trinityrnaseq/trinityrnaseq/wiki" rel="noreferrer">Trinity</a>, and would like to apply an additional filter to exclude transcripts that don't have good support from the RNASeq reads. This is with the goal of trying to reduce the initial dataset down to something a bit more manageable (I have about 300k transcripts covering 250Mb in total, but would prefer about 1/10 of that number).</p>
  
  <p>I see that Trinity has a <a href="https://github.com/trinityrnaseq/trinityrnaseq/wiki/RNA-Seq-Read-Representation-by-Trinity-Assembly" rel="noreferrer">workflow using Bowtie2</a> for evaluating the quality of assembled transcripts, but Bowtie2 is getting a bit old in the mapping world now. Are any of the super-fast transcript mappers (e.g. <a href="https://combine-lab.github.io/salmon/" rel="noreferrer">Salmon</a> or <a href="https://pachterlab.github.io/kallisto/about" rel="noreferrer">Kallisto</a>) appropriate to use for transcript filtering?</p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-13T08:37:09.387Z
Comments.Text: It's encouraged for Stack Exchange. There's even a [FAQ Answer](https://bioinformatics.stackexchange.com/help/self-answer) about it. It serves as a documentation of a tricky problem/solution I have, and also allows others to contribute and offer other solutions that I wasn't aware of (like your own).
Comments.UserDisplayName: gringer
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-18T14:16:38.783Z
Comments.Text: You may try using hisat2 instead of bowtie2 in the workflow you mention: hisat2 almost has the same user interface as bowtie2, and could be faster.
Comments.UserDisplayName: bli
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-07-18T14:30:27.753Z
Comments.Text: I've tried both HISAT2 and Bowtie2 for mapping. HISAT2 is splice-aware, so *should* be better for mapping RNASeq reads to a perfect genome, but Bowtie2 gives slightly better BUSCO scores for sequencing error correction based on mapping (in my particular situation).
Comments.UserDisplayName: gringer
Tags: <rna-seq><read-mapping>
Title: Filter Trinity transcriptome based on RNASeq reads
ThreadText: # Filter Trinity transcriptome based on RNASeq reads
  
  (Asked by: username_0 on 2017-07-12)
  
  I have recently generated a genome-guided transcriptome with [Trinity](https://github.com/trinityrnaseq/trinityrnaseq/wiki), and would like to apply an additional filter to exclude transcripts that don't have good support from the RNASeq reads. This is with the goal of trying to reduce the initial dataset down to something a bit more manageable (I have about 300k transcripts covering 250Mb in total, but would prefer about 1/10 of that number).
  
  I see that Trinity has a [workflow using Bowtie2](https://github.com/trinityrnaseq/trinityrnaseq/wiki/RNA-Seq-Read-Representation-by-Trinity-Assembly) for evaluating the quality of assembled transcripts, but Bowtie2 is getting a bit old in the mapping world now. Are any of the super-fast transcript mappers (e.g. [Salmon](https://combine-lab.github.io/salmon/) or [Kallisto](https://pachterlab.github.io/kallisto/about)) appropriate to use for transcript filtering?
  
  **Tags:** rna-seq, read-mapping
  
  > **Comments:**
  > **username_0**: It's encouraged for Stack Exchange. There's even a [FAQ Answer](https://bioinformatics.stackexchange.com/help/self-answer) about it. It serves as a documentation of a tricky problem/solution I have, and also allows others to contribute and offer other solutions that I wasn't aware of (like your own).  
  > **username_1**: You may try using hisat2 instead of bowtie2 in the workflow you mention: hisat2 almost has the same user interface as bowtie2, and could be faster.  
  > **username_0**: I've tried both HISAT2 and Bowtie2 for mapping. HISAT2 is splice-aware, so *should* be better for mapping RNASeq reads to a perfect genome, but Bowtie2 gives slightly better BUSCO scores for sequencing error correction based on mapping (in my particular situation).  
  
  ---
  
  ## 3 answers
  
  ### Answer by username_0 on 2017-07-12. Score: 6
  
  I used a combination of BUSCO and Salmon to filter transcripts based on their abundance in the RNASeq read dataset. The approach was roughly as follows:
  
  1. Run BUSCO in short/transcript mode on the Trinity-generated sequences
  2. Run Salmon on the Trinity-generated sequences, using the RNASeq reads that were used to generate the transcriptome
  3. Merge the BUSCO full results table with the Salmon results table to determine distribution of the number of reads mapped for BUSCO genes
  4. Choose an appropriate "real" read count threshold based on the read count distribution for BUSCO genes
  5. Sub-select the transcripts based on this threshold, using the Salmon results table
  
  ## Step 1 -- BUSCO
  
  ```
  python ~/install/busco/BUSCO.py -i ../Trinity-GG.fasta -o BUSCO_Trinity-GG_nematodes -l ~/install/busco/nematoda_odb9 -m tran -c 10
  
  ```
  
  ## Step 2 -- Salmon
  
  ```
  ~/install/salmon/Salmon-0.8.2_linux_x86_64/bin/salmon index -t Trinity-GG.fasta -i Trinity-GG.fasta.sai
  ~/install/salmon/Salmon-0.8.2_linux_x86_64/bin/salmon quant -i Trinity-GG.fasta.sai -1 reads-all_R1_trimmed.fastq.gz -2 reads-all_R2_trimmed.fastq.gz -p 10 -o quant/reads-all_quant -l A
  
  ```
  
  ## Step 3 -- Results merge (using R)
  
  ```
  #!/usr/bin/Rscript
  busco.table <-
      read.delim("BUSCO/run_BUSCO_Trinity-GG_nematodes/full_table_BUSCO_T-BNOCFED_nematodes.tsv", comment.char="#",
                 header=FALSE, stringsAsFactors = FALSE,
                 col.names=c("Busco id","Status","Sequence","Score","Length"));
  busco.table <- subset(busco.table, Status != "Missing");
  
  count.table <- read.delim("quant/reads-all_quant_ISR/quant.sf",
                            stringsAsFactors=FALSE);
  
  busco.expr.df <- merge(busco.table, count.table,
                         by.x="Sequence", by.y="Name", all.x=TRUE);
  
  ```
  
  ## Step 4 -- Count thresholding (visual exploration)
  
  ```
  pdf("busco_NumReads.pdf", paper="a4r", width=11, height=8);
  options(scipen=10);
  busco.numReads <- sort(tapply(busco.expr.df$NumReads, busco.expr.df$Busco.id,
                                function(x){
                                    exp(mean(log(x[x>0])))}), decreasing=TRUE);
  par(mfrow=c(1,2));
  plot(x=seq(0,1,length.out=length(busco.numReads)),
       busco.numReads, log="y", ylab="Number of Mapped Reads",
       yaxt="n",
       xlab="Proportion of Complete BUSCO sequences");
  axis(2, at=rep(c(1,2,5),each=6) * 10^(0:5), las=2, cex.axis=0.71);
  abline(col="#00000020", h=rep(1:9,each=6) * 10^(0:5), lwd=0.1);
  abline(col="#80808080", h=rep(1,each=6) * 10^(0:5), lwd=2);
  rect(ybottom=10^(par("usr")[3]), ytop=10^(par("usr")[4]),
       xleft=0.9, xright=par("usr")[2],
       col="#80808020", border=NA);
  plot(x=seq(0,1,length.out=length(busco.numReads)), xlim=c(0.9,1),
       busco.numReads, log="y", ylab="Number of Mapped Reads", yaxt="n",
       xlab="Proportion of Complete BUSCO sequences");
  axis(2, at=rep(c(1,2,5),each=6) * 10^(0:5), las=2, cex.axis=0.71);
  abline(col="#00000020", h=rep(1:9,each=6) * 10^(0:5), lwd=0.1);
  abline(col="#80808080", h=rep(1,each=6) * 10^(0:5), lwd=2);
  invisible(dev.off());
  
  ```
  
  [![choosing read count threshold from BUSCO genes](https://i.sstatic.net/7iw47.png)](https://i.sstatic.net/7iw47.png)
  
  Looks like a value of 50 will be suitable.
  
  ## Step 5 -- Sub-select the transcripts
  
  ```
  cat(file="HighCount_Transcripts_Num50.txt",
      subset(count.table, NumReads >= 50)$Name, sep="\n");
  
  ```
  
  Then filter using a suitable [FASTA selection utility](https://github.com/username_0/bioinfscripts/blob/master/fastx-fetch.pl):
  
  ```
  pv Trinity-GG.fasta | ~/scripts/fastx-fetch.pl -i HighCount_Transcripts_Num50.txt > HighCount50_Trinity-GG.fasta
  
  ```
  
  > **Comments:**
  > **username_3**: I really like the apporach! How do you handle/treat replicates?  
  > **username_0**: Do you mean identical transcripts (or sub-transcripts)? I've used a very strict cd-hit-est clustering to clean up a few (98% identity), but the BUSCO results still indicate that there's a lot of duplication going on. I'm hesitant to use a looser identity threshold because I don't want to include any of the real gene duplications.  
  
  ---
  
  ### Answer by username_2 on 2017-07-13. Score: 3
  
  I would recommend [Transrate](http://hibberdlab.com/transrate/index.html). It generates a score for each transcript based on many metrics (read their paper [here](http://dx.doi.org/10.1101/gr.196469.115)) and, using that score, gives a predicted cut-off to produce the optimal transcript set.
  
  
  ---
  
  ### Answer by username_4 on 2018-12-20. Score: 0
  
  [DETONATE](http://deweylab.biostat.wisc.edu/detonate/) offers a pipeline for filtering in this way as well. But the aligner used by default is Bowtie2...
  
  
  ---
  
  