Question: (Closed) Tophat Error When "Reporting Output Tracks"
0
gravatar for Pierre Lindenbaum
6.3 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum118k wrote:

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

tophat • 5.2k views
ADD COMMENTlink modified 6.3 years ago • written 6.3 years ago by Pierre Lindenbaum118k

(I'm away from the lab) my fastqs are gzipped. Could it be a problem?

ADD REPLYlink written 6.3 years ago by Pierre Lindenbaum118k

It should work with TopHat even if they are. And you can see from the output that it has finished mapping - the failure comes at the reporting stage.

ADD REPLYlink modified 6.3 years ago • written 6.3 years ago by Mikael Huss4.6k

cross posted on seqanswers: http://seqanswers.com/forums/showthread.php?p=91519#post91519

ADD REPLYlink written 6.3 years ago by Pierre Lindenbaum118k

and I shouldn't have used the option '-c' in bowtie-index . I'll close that question.

ADD REPLYlink written 6.3 years ago by Pierre Lindenbaum118k
0
gravatar for Mikael Huss
6.3 years ago by
Mikael Huss4.6k
Stockholm
Mikael Huss4.6k wrote:

I don't know what is wrong. But I do know there can be some sort of stochastic behavior at this step. I have had runs fail in a similar way, restarted them, and suddenly they work. This typically happens at the Reporting output tracks stage. I suggest you try again.

Perhaps a not very helpful answer, but I wanted to give you some hope :-)

ADD COMMENTlink written 6.3 years ago by Mikael Huss4.6k
Please log in to add an answer.
The thread is closed. No new answers may be added.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1946 users visited in the last hour