text: user73: I have recently generated a genome-guided transcriptome with Trinity, 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 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 or Kallisto) appropriate to use for transcript filtering?
  
  user1092: I would recommend Transrate. It generates a score for each transcript based on many metrics (read their paper here) and, using that score, gives a predicted cut-off to produce the optimal transcript set.
  
  user73: 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:
  
  Run BUSCO in short/transcript mode on the Trinity-generated sequences
  Run Salmon on the Trinity-generated sequences, using the RNASeq reads that were used to generate the transcriptome
  Merge the BUSCO full results table with the Salmon results table to determine distribution of the number of reads mapped for BUSCO genes
  Choose an appropriate "real" read count threshold based on the read count distribution for BUSCO genes
  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());
  
  
  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:
  pv Trinity-GG.fasta | ~/scripts/fastx-fetch.pl -i HighCount_Transcripts_Num50.txt > HighCount50_Trinity-GG.fasta
  
  
original_text: 73: <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>
   1092: <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>
   73: <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.stack.imgur.com/7iw47.png" rel="noreferrer"><img src="https://i.stack.imgur.com/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>
  