<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>

CC BY-SA 3.0
2017-07-13T08:35:41.727Z
2017-07-13T08:44:04.673Z
2017-07-13T08:44:04.673Z
ithinkiam
<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>

CC BY-SA 3.0
2017-07-12T09:39:48.067Z
I really like the apporach! How do you handle/treat replicates?
Kristoffer Vitting-Seerup
CC BY-SA 3.0
2017-07-12T11:12:42.807Z
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.
gringer
CC BY-SA 3.0
2017-07-12T08:41:00.213Z
2017-07-12T08:54:27.503Z
2017-07-12T08:54:27.503Z
gringer
<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>

CC BY-SA 4.0
2018-12-20T16:18:01.330Z
2018-12-20T16:18:01.330Z
winni2k
<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>

CC BY-SA 3.0
2017-07-13T08:37:09.387Z
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).
gringer
CC BY-SA 3.0
2017-07-18T14:16:38.783Z
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.
bli
CC BY-SA 3.0
2017-07-18T14:30:27.753Z
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).
gringer
<rna-seq><read-mapping>
Filter Trinity transcriptome based on RNASeq reads
# 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:

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

> **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...


---

