Problem using the Trinity script align_and_estimate_abundance.pl with RSEM

[gwilymh]

Hi,
I have assembled RNA-Seq data do novo in Trinity, and am now trying to estimate the abundance of the transcripts that were identified. I used TopHat to align the RNA-Seq fastq reads against the Trinity-generated fasta file of transcripts.

All attempts to use the align_and_estimate_abundance.pl script that is packaged with Trinity have failed, however.

My input code is:

/usr/local/trinity/util/align_and_estimate_abundance.pl \
--transcripts "sampleName.Trinity.fasta" \
--seqType fq \
--left "/filepath/reads.1.fastq" \
--right "/filepath/reads.2.fastq" \
--est_method RSEM \
--aln_method accepted_hits.bam \
--thread_count 1 \
--trinity_mode \
--fragment_length 180 \
--output_dir "./" \
--output_prefix "S1"

sampleName.Trinity.fasta was generated in Trinity, and accepted_hits.bam generated in TopHat

I always get the following error message:

CMD: rsem-calculate-expression --paired-end -p 4 --bam accepted_hits.bam /filepath/S1_RSEM/CB1_S1_L001.Trinity.fasta.RSEM S1
/usr/local/bin/rsem-parse-alignments /filepath/S1_RSEM/CB1_S1_L001.Trinity.fasta.RSEM S1.temp/S1 S1.stat/S1 b accepted_hits.bam -t 3 -tag XM
RSEM does not support gapped alignments, sorry!
"/usr/local/bin/rsem-parse-alignments /filepath/S1_RSEM/CB1_S1_L001.Trinity.fasta.RSEM S1.temp/S1 S1.stat/S1 b accepted_hits.bam -t 3 -tag XM" failed! Plase check if you provide correct parameters/options for the pipeline!
Error, cmd: rsem-calculate-expression --paired-end -p 4 --bam accepted_hits.bam /filepath/S1_RSEM/CB1_S1_L001.Trinity.fasta.RSEM S1 died with ret: 65280 at /usr/local/trinity/util/align_and_estimate_abundance.pl line 620

I have tried using --alignment_method bowtie2 instead of accepted_hits.bam but got the same error message.

Replace the --seqType flag from fq to fa changed the result -t 3 to -t 2 in the output, but otherwise produced the same error.

Any suggestions as to why align_and_estimate_abundance.pl does not seem to be working correctly?

[GenoMax]

Looking at the help page (http://trinityrnaseq.sourceforge.net...ml#RSEM_output) you need to use bowtie as the aligner (since RSEM needs ungapped alignments). So tophat output won't work unless you choose bowtie v.1 as the aligner for tophat.

[jordi]

Hi gwilymh / all!
I got the same error message. However I mapped reads back to transcripts using Bowtie2 allowing just the best alignment, as my reference transcriptome was a subset of filtered sequences by their annotation.
In fact I don't know what exactly the error refers to. Is there a way to see the ungapped alignment from the sam/bam file??
Moreover, testing the resulting bam file from the alignment with rsem-sam-validator, the bam file looks fine since I get "The input file is valid!".
Could anyone shed light over this?
Thanks!

[jordi]

Update: using samtools depth I was unable to find a position with 0 coverage. In fact it seems that I get the error message with rsem-parse-alignments since:

RSEM does not support gapped alignments, sorry!
"rsem-parse-alignments BLAST_5-nucleotidase.orf.RSEM BLAST_5-nucleotidase_DE.temp/BLAST_5-nucleotidase_DE BLAST_5-nucleotidase_DE.stat/BLAST_5-nucleotidase_DE b BLAST_5-nucleotidase_cobalt.sam_RSEM.bam -t 1 -tag XM" failed! Plase check if you provide correct parameters/options for the pipeline!

Advices would be appreciate