Filter Trinity transcriptome based on RNASeq reads
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 &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);


Step 4 -- Count thresholding (visual exploration)

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());




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 &gt;= 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 &gt; HighCount50_Trinity-GG.fasta