Why Does Tophat Fail On Me?

Question

[learnwell]

Hi BioStar,

I am fairly new to RNA-seq. I am trying to use Tophat for RNA-seq read alignment, and am trying to save the Bowtie index built as part of the alignment process. According to the Tophat Manual the way to do this is to use the -G flag, followed by the --transcriptome-index flag. To that end I have written a shell script that makes use of this feature (relevant portion below):

ebwt=/home/Homo_sapiens/UCSC/hg19/Sequence/Bowtie2Index/genome
gtf=/home/Homo_sapiens/UCSC/hg19/Annotation/Genes/genes.gtf
refFastq=/home/Homo_sapiens/UCSC/hg19/Sequence/Bowtie2Index/genome.fa
known=/home/Homo_sapiens/UCSC/hg19/Annotation/Genes/known

output="tophat_20120917"

dir="/home/"$lane"/"

cd $dir
tophat -G $gtf --transcriptome-index=/home/known  --solexa1.3-quals -o $output --library-type fr-unstranded $ebwt $left 1>stdout.txt 2>stderr.txt

However, when I run the script, I get the following output to stderrr.txt:

[2012-09-26 11:16:48] Checking for Bowtie index files
Error: Could not find Bowtie 2 index files (/home/known.*.bt2)

What am I doing wrong? It seems like the --transcriptome-index flag is overriding the -G flag, but this shouldn't be happening.

Any advice is appreciated!

Answer 1

[Istvan Albert]

With this tool the index base is specified as the first non-flag parameter.

Since you use the = sign with the options that throws off the command line parser. It now thinks that what you are after is setting --transcriptome-index to = then wants to use the next value home/known as index base.

The solution is to remove the = sign from the command line:

tophat -G $gtf --transcriptome-index /home/known  --solexa1.3-quals -o $output --library-type fr-unstranded $ebwt $left 1>stdout.txt 2>stderr.txt

Disclaimer: I did not actually try or verify the above - it but appears to be a simple logical explanation to what you see