Off topic:Tophat Error When "Reporting Output Tracks"

Question

[Pierre Lindenbaum]

I'm learning tophat . here is the Makefile I'm using:

include tools.mk
REF=./chr22
REF.bowtie=$(foreach S,1.bt2 2.bt2 3.bt2 4.bt2 rev.1.bt2 rev.2.bt2, $(addsuffix .$S,$(REF)) )
GTF=Homo_sapiens.GRCh37.69.gtf
export PATH:=$(PATH):${BOWTIE2.dir}:${samtools.dir}

JETER: s1_r1.fastq.gz s1_r2.fastq.gz ${REF.bowtie} $(REF).fa $(GTF)
    echo ${PATH}
    $(TOPHAT.dir)/tophat2 -G $(GTF) -o $@ $(REF) s1_r1.fastq.gz s1_r2.fastq.gz
s1_r1.fastq.gz:
    gunzip -c data/1_7_1.fastq.gz | head -n 400000 | gzip --best > $@
s1_r2.fastq.gz:
    gunzip -c data/1_7_2.fastq.gz | head -n 400000 | gzip --best > $@

$(GTF):
    gunzip -c /commun/data/pubdb/ensembl/release-69/gtf/homo_sapiens/Homo_sapiens.GRCh37.69.gtf.gz| egrep '^22' > $@

$(REF).fa: 
    curl -o $(REF).fa.gz "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/chr22.fa.gz"
    gunzip  -c $(REF).fa.gz | sed 's/^>chr/>/' > $@
    rm $(REF).fa.gz


${REF.bowtie}:$(REF).fa
    ${BOWTIE2.dir}/bowtie2-build  -f -c $(REF) $(REF)
    ${BOWTIE2.dir}/bowtie2-inspect -s $(REF)

When 'make' is invoked, I get the following messages:

[lindenb@master tophat]$ make -I /commun/data/packages 

/commun/data/packages/tophat-2.0.6.Linux_x86_64/tophat2 -G Homo_sapiens.GRCh37.69.gtf -o JETER ./chr22 s1_r1.fastq.gz s1_r2.fastq.gz

[2012-12-12 15:33:51] Beginning TopHat run (v2.0.6)
-----------------------------------------------
[2012-12-12 15:33:51] Checking for Bowtie
          Bowtie version:     2.0.2.0
[2012-12-12 15:33:51] Checking for Samtools
        Samtools version:     0.1.18.0
[2012-12-12 15:33:51] Checking for Bowtie index files
[2012-12-12 15:33:51] Checking for reference FASTA file
[2012-12-12 15:33:51] Generating SAM header for ./chr22
    format:         fastq
    quality scale:     phred33 (default)
[2012-12-12 15:33:51] Reading known junctions from GTF file
[2012-12-12 15:33:51] Preparing reads
     left reads: min. length=76, max. length=76, 100000 kept reads (0 discarded)
    right reads: min. length=76, max. length=76, 99889 kept reads (111 discarded)
[2012-12-12 15:33:55] Creating transcriptome data files..
[2012-12-12 15:34:21] Building Bowtie index from Homo_sapiens.GRCh37.69.fa
[2012-12-12 15:35:28] Mapping left_kept_reads to transcriptome Homo_sapiens.GRCh37.69 with Bowtie2 
[2012-12-12 15:35:58] Mapping right_kept_reads to transcriptome Homo_sapiens.GRCh37.69 with Bowtie2 
[2012-12-12 15:36:27] Resuming TopHat pipeline with unmapped reads
[2012-12-12 15:36:27] Mapping left_kept_reads.m2g_um to genome chr22 with Bowtie2 
[2012-12-12 15:36:34] Mapping left_kept_reads.m2g_um_seg1 to genome chr22 with Bowtie2 (1/3)
[2012-12-12 15:36:35] Mapping left_kept_reads.m2g_um_seg2 to genome chr22 with Bowtie2 (2/3)
[2012-12-12 15:36:37] Mapping left_kept_reads.m2g_um_seg3 to genome chr22 with Bowtie2 (3/3)
[2012-12-12 15:36:39] Mapping right_kept_reads.m2g_um to genome chr22 with Bowtie2 
[2012-12-12 15:36:45] Mapping right_kept_reads.m2g_um_seg1 to genome chr22 with Bowtie2 (1/3)
[2012-12-12 15:36:47] Mapping right_kept_reads.m2g_um_seg2 to genome chr22 with Bowtie2 (2/3)
[2012-12-12 15:36:48] Mapping right_kept_reads.m2g_um_seg3 to genome chr22 with Bowtie2 (3/3)
[2012-12-12 15:36:50] Searching for junctions via segment mapping
[2012-12-12 15:36:53] Retrieving sequences for splices
[2012-12-12 15:37:00] Indexing splices
[2012-12-12 15:37:52] Mapping left_kept_reads.m2g_um_seg1 to genome segment_juncs with Bowtie2 (1/3)
[2012-12-12 15:38:44] Mapping left_kept_reads.m2g_um_seg2 to genome segment_juncs with Bowtie2 (2/3)
[2012-12-12 15:39:33] Mapping left_kept_reads.m2g_um_seg3 to genome segment_juncs with Bowtie2 (3/3)
[2012-12-12 15:40:23] Joining segment hits
[2012-12-12 15:40:27] Mapping right_kept_reads.m2g_um_seg1 to genome segment_juncs with Bowtie2 (1/3)
[2012-12-12 15:41:17] Mapping right_kept_reads.m2g_um_seg2 to genome segment_juncs with Bowtie2 (2/3)
[2012-12-12 15:42:07] Mapping right_kept_reads.m2g_um_seg3 to genome segment_juncs with Bowtie2 (3/3)
[2012-12-12 15:42:56] Joining segment hits
[2012-12-12 15:42:59] Reporting output tracks
    [FAILED]
Error running /commun/data/packages/tophat-2.0.6.Linux_x86_64/tophat_reports --min-anchor 8 --splice-mismatches 0 --min-report-intron 50 --max-report-intron 500000 --min-isoform-fraction 0.15 --output-dir JETER/ --max-multihits 20 --max-seg-multihits 40 --segment-length 25 --segment-mismatches 2 --min-closure-exon 100 --min-closure-intron 50 --max-closure-intron 5000 --min-coverage-intron 50 --max-coverage-intron 20000 --min-segment-intron 50 --max-segment-intron 500000 --read-mismatches 2 --read-gap-length 2 --read-edit-dist 2 --read-realign-edit-dist 3 --max-insertion-length 3 --max-deletion-length 3 -z gzip --inner-dist-mean 50 --inner-dist-std-dev 20 --gtf-annotations Homo_sapiens.GRCh37.69.gtf --gtf-juncs JETER/tmp/Homo_sapiens.juncs --no-closure-search --no-coverage-search --no-microexon-search --sam-header JETER/tmp/chr22_genome.bwt.samheader.sam --report-discordant-pair-alignments --report-mixed-alignments --samtools=/commun/data/packages/samtools-0.1.18/samtools --bowtie2-max-penalty 6 --bowtie2-min-penalty 2 --bowtie2-penalty-for-N 1 --bowtie2-read-gap-open 5 --bowtie2-read-gap-cont 3 --bowtie2-ref-gap-open 5 --bowtie2-ref-gap-cont 3 ./chr22.fa JETER/junctions.bed JETER/insertions.bed JETER/deletions.bed JETER/fusions.out JETER/tmp/accepted_hits JETER/tmp/left_kept_reads.m2g.bam JETER/tmp/left_kept_reads.bam JETER/tmp/right_kept_reads.m2g.bam JETER/tmp/right_kept_reads.bam
Loaded 8026 junctions

make: *** [JETER] Error 1

what's wrong ?

Thank you.

EDIT: I've narrowed the error:

The exception is raised from

 (int)length(*ref_str)

in

src/tophat_reports.cpp

in the line

   if (new_left >= 0 && new_bh.right() <= (int)length(*ref_str))

it seems that ref_str is NULL when length is invoked

Pierre