Error in Bam file for mitoseek

Question

[ammarsabir15]

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

---

---

[WouterDeCoster]

This looks like a more informative question, but it's quite similar to C: Mitoseek not installing correctly isn't it?

[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.

[WouterDeCoster]

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

[ammarsabir15]

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??

[WouterDeCoster]

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

[ammarsabir15]

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

---

---

[WouterDeCoster]

Then first sort and index you bam file using samtools.

[ammarsabir15]

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

---

---

[WouterDeCoster]

So that's not the same chromosome.

[ammarsabir15]

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.

[WouterDeCoster]

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.