This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Error in Bam file for mitoseek

I want to detect heteroplasmy in mitochondrial genomes . I have successfully installed mitoseek and mitoseek is working for default datasets .

This is the command and output for example datasets.

ammar@ammar-HP-15-Notebook-PC:~/MitoSeek-1.3$ perl mitoSeek.pl -i mito1.bam

[Tue Aug 16 22:20:44 2016] [WARN] folder mito1 already exists,will delete it first

==================================================

Steps will be run:

1,Extracting reads in mitochondria from '/home/ammar/MitoSeek-1.3/mito1.bam' (Output: mito1.bam)

2,Checking Reads number in 'mito1.bam'

3,Moving mitochondrial genome '/home/ammar/MitoSeek-1.3/Resources/genome/hg19.fasta' into mito1 (Output:mito.fasta)

4,Pileuping 'mito1.bam' (Output:mito1.pileup)

5,Parsing pileup file of 'mito1.pileup' (Output: mito1_basecall.txt)

6,Getting quality metrics on 'mito1.bam' 

7,Detecting heteroplasmy from 'mito1.bam' (Output: mito1_heteroplasmy.txt)

8,Detecting structure variants from 'mito1.bam' (Output: mito1_structure_discordant_mates.txt | mito1_structure_large_deletion.sam)

9,Generating report (Output: mitoSeek.html)

==================================================

==================================================

Start analyzing:

[Tue Aug 16 22:20:44 2016] [WARN] It seems like the mitochondrial reference used in the bam is rCRS, not hg19

[Tue Aug 16 22:20:44 2016] 1,Extracting reads in mitochondria from '/home/ammar/MitoSeek-1.3/mito1.bam' (Output: mito1.bam)

[Tue Aug 16 22:20:45 2016] 2,Checking Reads number in 'mito1.bam' n=46123

[Tue Aug 16 22:20:45 2016] 3,Moving mitochondrial genome '/home/ammar/MitoSeek-1.3/Resources/genome/hg19.fasta' into mito1 (Output:mito.fasta)

[Tue Aug 16 22:20:45 2016] 4,Pileuping 'mito1.bam' (Output:mito1.pileup)

[Tue Aug 16 22:20:49 2016] 5,Parsing pileup file of 'mito1.pileup' (Output: mito1_basecall.txt)

[Tue Aug 16 22:20:49 2016] 6,Getting quality metrics on 'mito1.bam'

[Tue Aug 16 22:22:48 2016] 7,Detecting heteroplasmy from 'mito1.bam' (Output: mito1_heteroplasmy.txt)

[Tue Aug 16 22:22:52 2016] 8,Detecting structure variants from 'mito1.bam' (Output: mito1_structure_discordant_mates.txt | mito1_structure_large_deletion.sam)

[Tue Aug 16 22:22:53 2016] 9,Generating report (Output: mitoSeek.html)

==================================================

Total Running Time: 0h:02m:09s

And this is the command for my dataset

ammar@ammar-HP-15-Notebook-PC:~/MitoSeek-1.3$ perl mitoSeek.pl -i data88.bam

[Tue Aug 16 23:09:38 2016] [WARN] folder data88 already exists,will delete it first

==================================================

Steps will be run: 1,Extracting reads in mitochondria from '/home/ammar/MitoSeek-1.3/data88.bam' (Output: mito1.bam)

2,Checking Reads number in 'mito1.bam'

3,Moving mitochondrial genome '/home/ammar/MitoSeek-1.3/Resources/genome/hg19.fasta' into data88 (Output:mito.fasta)

4,Pileuping 'mito1.bam' (Output:mito1.pileup)

5,Parsing pileup file of 'mito1.pileup' (Output: mito1_basecall.txt)

6,Getting quality metrics on 'mito1.bam' 

7,Detecting heteroplasmy from 'mito1.bam' (Output: mito1_heteroplasmy.txt)

8,Detecting structure variants from 'mito1.bam' (Output: mito1_structure_discordant_mates.txt | mito1_structure_large_deletion.sam)

9,Generating report (Output: mitoSeek.html)

==================================================

Start analyzing: [Tue Aug 16 23:09:38 2016] [ERROR] No mithochondrial chromosome detected from the header

snp genome next-gen
• 2,315 views • link
written by ammarsabir15 ammarsabir15

Yes, actually i tried to ask many things in the same post but it failed. Know I managed to install mitoseek but this issue is still there.

• 1 views • link
written by ammarsabir15 ammarsabir15

You checked the notation of chrMT/chrM/M/MT in fasta and bam file is the same?

• 1 views • link
written by WouterDeCoster WouterDeCoster

I have made bam file from fastq and in fastq there is no such notation, there is only information about reads. And would you please elaborate what is this notation ? Is it a tag that tells the sequence is from mitochondrial genome??

• 1 views • link
written by ammarsabir15 ammarsabir15

You created the bam file by aligning the fastq to a genome. Check the chromosome notations using grep '>' /home/ammar/MitoSeek-1.3/Resources/genome/hg19.fasta and samtools idxstats data88.bam

• 1 views • link
written by WouterDeCoster WouterDeCoster

I checked it and following are results:

ammar@ammar-HP-15-Notebook-PC:~/MitoSeek-1.3$ grep '>' /home/ammar/MitoSeek- 1.3/Resources/genome/hg19.fasta

gi|17981852|ref|NC_001807.4| Homo sapiens mitochondrion, complete genome

ammar@ammar-HP-15-Notebook-PC:~/MitoSeek-1.3$ samtools idxstats data88.bam

[bam_idxstats] fail to load the index

• 1 views • link
written by ammarsabir15 ammarsabir15

Then first sort and index you bam file using samtools.

• 1 views • link
written by WouterDeCoster WouterDeCoster

Thanks, I have indexed and sorted my bam file but the error remains to be the same as shown in very start of the post .

ammar@ammar-HP-15-Notebook-PC:~/My analysis_2$ samtools sort data88.bam > data.sorted.bam

ammar@ammar-HP-15-Notebook-PC:~/My analysis_2$ samtools index data.sorted.bam

ammar@ammar-HP-15-Notebook-PC:~/My analysis_2$ samtools idxstats data.sorted.bam

gi|251831106|ref|NC_012920.1| 16569 12095 0

  • 0 0 9025
• 1 views • link
written by ammarsabir15 ammarsabir15

So that's not the same chromosome.

• 1 views • link
written by WouterDeCoster WouterDeCoster

You haven't explained me the use of above command .

I haven't used /home/ammar/MitoSeek- 1.3/Resources/genome/hg19.fasta as reference genome but the used following reference genome.

This is the header of the genome.

gi|251831106|ref|NC_012920.1| Homo sapiens mitochondrion, complete genome

    ammar@ammar-HP-15-Notebook-PC:~/Myanalysis_2$ grep  '>'              

   /home/ammar/Myanalysis_2/reference.fa

  >gi|251831106|ref|NC_012920.1| Homo sapiens mitochondrion, complete genome

And as you can see from previous comment these two are same. Sorry for the confusion being created. Actually i didn't knew the purpose of command you have asked me to implement.

• 1 views • link
written by ammarsabir15 ammarsabir15

And earlier you were using /home/ammar/MitoSeek-1.3/Resources/genome/hg19.fasta...

In addition, I'm not required to explain you the use of any command, you can google it for example.

I can't code in perl but if I understand the script correctly it expects to find "chrM" as chromosome identifier (line1431 in https://github.com/riverlee/MitoSeek/blob/v1.3/mitoSeek.pl) so I think it would be the best if you repeated the mapping.

• 1 views • link
written by WouterDeCoster WouterDeCoster

0 answers

No answers yet.

Log in to answer this question.

More posts like this