Question: Mitoseek not installing correctly
0
gravatar for ammarsabir15
2.7 years ago by
ammarsabir1550
ammarsabir1550 wrote:

I have installed mitoseek1.3 on ubuntu 15.10 (64 bit) but when i run it from terminal it gives following error. The command i gave and error is shown below:



I made following changes in mitoSeek_new.pl

my $samtools = "/usr/local/bin/samtools"; #Where is the samtools file

my $mitomap = "mitomap.pl"; #Where is the samtools file



And following changes in mitomap.pl

my $samtools = "/usr/local/bin/samtools"; #Where is the samtools file

my $bwa = "/usr/bin/bwa";



ammar@ammar-HP-15-Notebook-PC:~/Downloads/mitoseek 1.3/MitoSeek-1.3$ sudo perl mitoSeek_new.pl -i data88.bam

[sudo] password for ammar:

[Sun Aug 14 01:30:47 2016] [ERROR] input bam file (-i)

'/workspace/guoy1/GeneTorrent/lij17/download/HNSC/DNA_NB/e448ac62-9d33-4583-b8ea-2da157cfe0c2/C495.TCGA-C'

does not exists13-08.1.bam

Usage: perl mitoSeek.pl -i inbam -i [bam] Input bam file

-j [bam] Input bam file2, if this file is provided, it will conduct somatic mutation mining, and it will be taken as normal tissue.

-t [input type] Type of the bam files, the possible choices are 1=exome, 2=whole genome, 3= RNAseq, 4 = mitochondria only,default = 1

-d [int] The minimum recommended depth requirement for detecting heteroplasmy. Lower depth will severely damage the confidence of heteroplasmy calling, default=50

-ch Produce circos plot input files and circos plot figure for heteroplasmic mutation, (-noch to turn off and -ch to turn on), default = on

-hp [int] Heteroplasmy threshold using [int] percent alternative allele observed, default = 5

-ha [int] Heteroplasmy threshold using [int] allele observed, default = 0

-alpha [float] Shape1 parameter of Beta prior distribution, default is 3.87 which is estimated from 600 BRCA samples

-beta [float] Shape1 parameter of Beta prior distribution, default is 174.28, which is estimated from 600 BRCA samples

-A If - A is used, the total read count is the total allele count of all allele observed. Otherwise, the total read count is the sum of major and minor allele counts. Default = off

-mmq [int] Minimum map quality, default =20

-mbq [int] Minimum base quality, default =20

-sb [int] Remove all sites with strand bias score in the top [int] %, default = 10

-cn Estimate relative copy number of input bam(s), does not work with mitochondria targeted sequencing bam files, (-noch to turn off and -ch to turn on) default = off.

-sp [int] Somatic mutation detection threshold,int = percent of alternative allele observed in tumor, default int=5

-sa [int] Somatic mutation detection threshold,int = number of alternative allele observed in tumor, default int=3

-cs Produce circos plot input files and circos plot figure for somatic mutation, (-nocs to turn off and -cs to turn on), default = off

-r [ref] The reference used in the bam file, the possible choices are hg19 and rCRS, default=hg19

-R [ref] The reference used in the output files, the possible choices are hg19 and rCRS, default=hg19

-str [int] Structure variants cutoff for those discordant mapping mates, int = number of spanning reads supporting this structure variants, default = 2

-strf [int] Structure variants cutoff for those large deletions, int = template size in bp, default=500

-QC Produce QC result, (--noQC to turn off and -QC to turn on), default=on

-samtools[samtools] Tell where is the samtools program, default is your mitoseek directory/Resources/samtools/samtools

-bwa [bwa] Tell where is the bwa program, default value is 'bwa' which is your $PATH

-bwaindex [bwaindex] Tell where is the bwa index of non-mitochondrial human genome, no default value

-advance Will get mitochondrial reads in an advanced way, generally followed by 1) Initially extract mitochrodrial reads from a bam file, then 2) remove those could be remapped to non-mitochondrial human genome by bwa. Advanced extraction needs -bwaindex option. Default extraction without removing step.

@------------------------------------------------------------@ | Statistical framework for heteroplasmy detection | |------------------------------------------------------------| | Fisher test for heterplasmy | |------------------------------------------------------------| | major minor | | observed n11 n12 | n1p | | expected n21 n22 | n2p | | ----------------- | | np1 np2 npp | | n21 = (n11+n12)(1-hp/100) in which hp is defined by -hp | | n22 = (n11+n12)hp/100 in which hp is defined by -hp | |------------------------------------------------------------| | Empirical Bayesian method | |------------------------------------------------------------| | _Inf | | / | | probability = | f(x)dx | | _/p | | probablity is the calculus of f(x) from p to Inf | | x=hp/100 in which hp is defined by -hp | | f(x) = 1/beta(b+A,a+B)x^(A+b-1)(1-x)^(B+a-1) | | in which a/b is the number of major/minor allele, | | A and B are estimated from 600 BRCA samples. | |------------------------------------------------------------| | For documentation, citation & bug-report instructions: | | https://github.com/riverlee/MitoSeek | @------------------------------------------------------------@

Program exist with Error


Anyone who knows?

ADD COMMENTlink modified 2.7 years ago by haansi50 • written 2.7 years ago by ammarsabir1550
1
gravatar for haansi
2.7 years ago by
haansi50
Austria
haansi50 wrote:

@WouterDeCoster is right - the name of the reference is required to be called "M" or "chrM" - if (/SN:M|SN:chrM/i) { see https://github.com/riverlee/MitoSeek/blob/v1.3/mitoSeek_new.pl

Another suggestion: you could also try to use lofreq see previous discussion in biostars

ADD COMMENTlink written 2.7 years ago by haansi50
0
gravatar for WouterDeCoster
2.7 years ago by
Belgium
WouterDeCoster38k wrote:

Can you post the output of readlink -f data88.bam? Is it a (broken?) symbolic link?

Unrelated, why do you use sudo with that command? It's always better not to use sudo if you don't have to, certainly with code you're not certain of what it does.

ADD COMMENTlink written 2.7 years ago by WouterDeCoster38k
0
gravatar for ammarsabir15
2.7 years ago by
ammarsabir1550
ammarsabir1550 wrote:

Thanks, I am a beginner in ngs as well as linux and actually I was using sudo because sometimes errors were coming like permission denied etc.

This is the output of readlink -f data88.bam.



ammar@ammar-HP-15-Notebook-PC:~/Downloads/mitoseek 1.3/MitoSeek-1.3$ readlink -f data88.bam

/home/ammar/Downloads/mitoseek 1.3/MitoSeek-1.3/data88.bam



ADD COMMENTlink written 2.7 years ago by ammarsabir1550

Please use ADD COMMENT to reply to earlier questions, that way we keep this thread logically ordered and structured.

Do the provided test datasets work? Is your samtools and bwa the recommended version?

ADD REPLYlink modified 2.7 years ago • written 2.7 years ago by WouterDeCoster38k

Note: mitoSeek_new.pl is an updated version of mitoseek in mitoseek1.3.

When I use the mitoSeek_new.pl with example datasets then output is same as appeared in above screen mentioned in very start of my post (i.e the program exist with error). And the samtools version which is available in resources (of mitoseek repository) is version is 0.1.18 which is very old while I installed samtools version 1.3.1. In the same way circos is 0.56 while i installed 0.67-7.



But when I used the example datasets (mito1.bam) given in https://github.com/riverlee/MitoSeek with mitoSeek.pl then the output is generarated :



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

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:

[Mon Aug 15 08:38:46 2016] [WARN] It seems like the mitochondrial reference used in the bam is rCRS, not hg19

[Mon Aug 15 08:38:46 2016] 1,Extracting reads in mitochondria from '/home/ammar/MitoSeek-1.3/mito1.bam' (Output:  
mito1.bam)

[Mon Aug 15 08:38:47 2016] 2,Checking Reads number in 'mito1.bam' n=46123

[Mon Aug 15 08:38:48 2016] 3,Moving mitochondrial genome '/home/ammar/
MitoSeek-      1.3/Resources/genome/hg19.fasta' into mito1 (Output:mito.fasta)

[Mon Aug 15 08:38:48 2016] 4,Pileuping 'mito1.bam' (Output:mito1.pileup)

[Mon Aug 15 08:38:51 2016] 5,Parsing pileup file of 'mito1.pileup' (Output: mito1_basecall.txt)

[Mon Aug 15 08:38:51 2016] 6,Getting quality metrics on 'mito1.bam'

[Mon Aug 15 08:40:54 2016] 7,Detecting heteroplasmy from 'mito1.bam' (Output: mito1_heteroplasmy.txt)

[Mon Aug 15 08:40:55 2016] 8,Detecting structure variants from 'mito1.bam' ( Output: mito1_structure_discordant_mates.txt | mito1_structure_large_deletion.sam)

[Mon Aug 15 08:40:55 2016] 9,Generating report (Output: mitoSeek.html)

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


As i am a beginner so I even don't know that whether the output was produced correctly or not but it did't produced circos plots, I configured circos and it was working rightly but it was not working with mitoseek. And for giving the right [ath of circod to mitoseek i made only following single change in Circoswrap.pm



_circosbin=>"/usr/bin/circos",     #where is the circos program


When I used mitoseek.pl with my dataset then it produced the following output:



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

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

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:

[Mon Aug 15 08:41:10 2016] [ERROR] No mithochondrial chromosome detected from the header

Program exist with Error



I also don't know what this error is refering to.!!!!!!!!!!!

ADD REPLYlink modified 2.7 years ago • written 2.7 years ago by ammarsabir1550
1

Could you check that in your bam file and fasta file the mitochondrial genome is named the same? Maybe one has chrMT and the other just M or something.

ADD REPLYlink written 2.7 years ago by WouterDeCoster38k
Please log in to add an answer.

Help
Access

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