6
$\begingroup$

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?

$\endgroup$
3
  • 3
    $\begingroup$ It's encouraged for Stack Exchange. There's even a FAQ 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). $\endgroup$
    – gringer
    Jul 13 '17 at 8:37
  • $\begingroup$ 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. $\endgroup$
    – bli
    Jul 18 '17 at 14:16
  • $\begingroup$ 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). $\endgroup$
    – gringer
    Jul 18 '17 at 14:30
6
$\begingroup$

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

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
$\endgroup$
2
  • $\begingroup$ I really like the apporach! How do you handle/treat replicates? $\endgroup$ Jul 12 '17 at 9:39
  • $\begingroup$ 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. $\endgroup$
    – gringer
    Jul 12 '17 at 11:12
3
$\begingroup$

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.

$\endgroup$
0
$\begingroup$

DETONATE offers a pipeline for filtering in this way as well. But the aligner used by default is Bowtie2...

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.