[miso-users] Error: tagBam call failed when running pe_utils.py --compute-insert-len
Yarden Katz
yarden at MIT.EDU
Sun Jul 7 15:05:42 EDT 2013
Hi,
It does look like a headers mismatch then. Your BAM file contains "chr" style chromosomes (e.g. "chr10" and not "10"). I believe your GFF, /directories/exons/Homo_sapiens.GRCh37.65.min_1000.const_exons.gff, is from Ensembl which would not contain chr-prefixes. Just look in that gff file and see what the chromosome entries are like. If they don't have chr, the operation will fail. All you need to do is generate a constitutive exons file from a UCSC gff which contains chromosome headers that match your .bam file.
See:
http://genes.mit.edu/burgelab/miso/docs/#human-mouse-gene-models-for-isoform-centric-analyses
If you're using hg19, you can use this GFF:
http://genes.mit.edu/burgelab/miso/annotations/ucsc_tables/hg19/ensGene.gff3
Use our exon_utils program to generate constitutive exons from this file and then rerun pe_utils with that, instead of the GRCh37 gff file.
--Yarden
On Jul 7, 2013, at 2:27 PM, Gu Mi wrote:
> Dear Yarden:
>
> Thank you for your prompt reply! Here is the output from the first command:
>
> @HD VN:1.0 SO:coordinate
> @SQ SN:chr1 LN:249250621
> @SQ SN:chr10 LN:135534747
> @SQ SN:chr11 LN:135006516
> @SQ SN:chr12 LN:133851895
> @SQ SN:chr13 LN:115169878
> @SQ SN:chr14 LN:107349540
> @SQ SN:chr15 LN:102531392
> @SQ SN:chr16 LN:90354753
> @SQ SN:chr17 LN:81195210
> @SQ SN:chr18 LN:78077248
> @SQ SN:chr19 LN:59128983
> @SQ SN:chr2 LN:243199373
> @SQ SN:chr20 LN:63025520
> @SQ SN:chr21 LN:48129895
> @SQ SN:chr22 LN:51304566
> @SQ SN:chr3 LN:198022430
> @SQ SN:chr4 LN:191154276
> @SQ SN:chr5 LN:180915260
> @SQ SN:chr6 LN:171115067
> @SQ SN:chr7 LN:159138663
> @SQ SN:chr8 LN:146364022
> @SQ SN:chr9 LN:141213431
> @SQ SN:chrM LN:16569
> @SQ SN:chrX LN:155270560
> @SQ SN:chrY LN:59373566
> @PG ID:TopHat CL:/informatics/tools/Linux-AS5/bin/tophat -o dir/Lane3 -g 1 --coverage-search --microexon -r 100 --phred64-quals --library-type fr-unstranded -p 4 -G gene_models/Homo_sapiens.GRCh37.72_norm.gtf --transcriptome-index=gene_models/transcripts dir/Human/bowtie/human_ref_genome dir/L1_1.fq.gz dir/L1_2.fq.gz VN:1.4.1
>
> The output for the second command is exactly the same as above. As they're not empty, I guess there might be other reasons for the error?
>
> Thanks again!
>
> Best,
> Gu
>
>
> On Sun, Jul 7, 2013 at 2:15 PM, Yarden Katz <yarden at mit.edu> wrote:
> Hi,
>
> See replies below:
>
> On Jul 7, 2013, at 2:04 PM, Gu Mi wrote:
>
>> Dear All:
>>
>> I am using MISO to test for differential exon usage between a control and a treatment group. I got an error when computing the insert length distribution using pe_utils.py --compute-insert-len. I list the steps I used below:
>>
>> 1. sort the BAM file from TopHat (by coordinate):
>> samtools sort control.bam control_sorted
>>
>> 2. index the BAM file:
>> samtools index control_sorted.bam control_sorted.bai
>>
>> 3. run pe_utils.py:
>> python pe_utils.py --compute-insert-len controlam /directories/exons/Homo_sapiens.GRCh37.65.min_1000.const_exons.gff --output-dir /directories/insert-dist/
>>
>> After the command above, I got the error message:
>>
>> Preparing to call bedtools 'tagBam'
>> tagBam -i control.bam -files /directories/exons/Homo_sapiens.GRCh37.65.min_1000.const_exons.gff -labels gff -intervals -f 1 | samtools view - -h | egrep '^@|:gff:' | samtools view - -Shb -o /directories/insert-dist/bam2gff_Homo_sapiens.GRCh37.65.min_1000.const_exons.gff/control.bam
>> [samopen] SAM header is present: 25 sequences.
>> [sam_read1] reference 'ID:TopHat CL:/informatics/tools/Linux-AS5/bin/tophat -o Lane3 -g 1 --coverage-search --microexon -r 100 --phred64-quals --library-type fr-unstranded -p 4 -G gene_models/Homo_sapiens.GRCh37.72_norm.gtf --transcriptome-index=gene_models/transcripts /directories/Genomes/NCBI_Jul-09-2012/Human/bowtie/human_ref_genome Lane3_1.fq.gz Lane3_2.fq.gz VN:1.4.1
>> ' is recognized as '*'.
>> [main_samview] truncated file.
>> Traceback (most recent call last):
>> File "/pe_utils.py", line 520, in <module>
>> main()
>> File "pe_utils.py", line 517, in main
>> sd_max=sd_max)
>> File "pe_utils.py", line 271, in compute_insert_len
>> output_dir)
>> File "exon_utils.py", line 185, in map_bam2gff
>> raise Exception, "Error: tagBam call failed."
>> Exception: Error: tagBam call failed.
>
> It sounds like tagBam could not map any of your reads to the GFF, which is usually a headers mismatch issue.
>
> Could you confirm that your BAM file contains Ensembl style headers, and not UCSC ones? What is the output of this command command for you:
>
> samtools view -H /directories/insert-dist/bam2gff_Homo_sapiens.GRCh37.65.min_1000.const_exons.gff/control.bam
>
> Also, what is the output of the following for you?
>
> tagBam -i control.bam -files /directories/exons/Homo_sapiens.GRCh37.65.min_1000.const_exons.gff -labels gff -intervals -f 1 | samtools view - -h | egrep '^@|:gff:'
>
> If this yields something empty, it means that tagBam could not match your BAM to the GFF, probably because of the chromosome headers mismatch issue that I mentioned.
>
> Best, --Yarden
>
>>
>> I used Homo_sapiens.GRCh37.72_norm.gtf from Ensembl as the annotation file when preparing my data, but downloaded
>> • Human genome (hg19) alternative events v2.0
>> from the MISO website and unzipped. I saw it is based on Homo_sapiens.GRCh37.65. Is this the version problem? If so, could anyone provide the latest GFF3 file for use? Thank you for your suggestions!
>>
>> Best,
>> Gu
>> _______________________________________________
>> miso-users mailing list
>> miso-users at mit.edu
>> http://mailman.mit.edu/mailman/listinfo/miso-users
>
>
>
>
> --
> Ph.D. student
> Department of Statistics
> Oregon State University
> ~~~~~~~~~~~~~~~~~~~~~~~~~
> "Some people dream of success, while others wake up and work hard at it." ~ Author Unknown.
>
More information about the miso-users
mailing list