search

dpuiu / MitoHPC

MIT License
11 stars 12 forks source link

No output after run "run.all.sh" #21

Open yuanhb1993 opened 5 days ago

yuanhb1993 commented 5 days ago

Hi, I am currently learning to use MitoHPC to analyze WGS data. After successfully setting up the environment and testing it according to the tutorial, the script did not perform the analysis or produce any output (just create some directory in output directory) when I ran it officially.

Could you please let me know if I missed anything or if there are any configuration errors?

`bash ./run.all.sh

Here is the content of my configuration file: init.sh: `cat init.sh

!/usr/bin/env bash

set -e

if [ -z $HP_SDIR ] ; then echo "Variable HP_SDIR not defined. Make sure you followed the SETUP ENVIRONMENT instructions" ; fi

##############################################################################################################

Program that setups the environmnet

Variable HP_SDIR must be pre-set !!!

##############################################################################################################

DIRECTORY PATHS

export HP_HDIR=readlink -f $HP_SDIR/.. #HP home directory export HP_BDIR=$HP_HDIR/bin/ #bin directory export HP_JDIR=$HP_HDIR/java/ #java directory

Human

export HP_RDIR=$HP_HDIR/RefSeq/ #reference directory

Mouse

export HP_RDIR=$HP_HDIR/RefSeqMouse/ #Mouse reference directory

###############################################################

SOFTWARE PATH

export PATH=$HP_SDIR:$HP_BDIR:$PATH

################################################################

ALIGNMNET REFERENCE

hs38DH(default)

export HP_RNAME=hs38DH export HP_RMT=chrM export HP_RNUMT="chr1:629084-634672 chr17:22521208-22521639" # 150bp reads chr2:32916199-32916630

export HP_RNUMT="chr1:629080-634925 chr1:76971223-76971280 chr2:148881723-148881858 chr5:80651184-80651597 chr11:10508892-10509738 chr13:109424123-109424381 chr17:22521208-22521639" # 100bp reads

export HP_RCOUNT=3366 # 195(hs38DH-no_alt); 194 (hs38DH-no_alt_EBV) export HP_RURL=ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa

hg19

export HP_RNAME=hg19

export HP_RMT=chrM

export HP_RNUMT="chr1:564465-569708 chr17:22020692-22020827"

export HP_RCOUNT=93

export HP_RURL=http://hgdownload.cse.ucsc.edu/goldenpath/hg19/bigZips/hg19.fa.gz

CHM13

export HP_RNAME=chm13

export HP_RMT=chrM

export HP_RNUMT="chr5:81136887-81137073 chr11:10594619-10594811 chr13:12167063-12168420 chr17:23209322-23209958"

export HP_RCOUNT=24

export HP_RURL=https://s3-us-west-2.amazonaws.com/human-pangenomics/T2T/CHM13/assemblies/chm13.draft_v1.1.fasta.gz

CHM13.v2; to compute NUMT regions

export HP_RURL=https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/009/914/755/GCF_009914755.1_T2T-CHM13v2.0/GCF_009914755.1_T2T-CHM13v2.0_genomic.fna.gz

Mouse

export HP_RNAME=mm39

export HP_RMT=chrM

export HP_RCOUNT=61

export HP_RNUMT="chr1:24650615-24655265 chr2:22477298-22480547 chr10:95913385-95913756 chr12:97028201-97028567 chr13:85274752-85275567"

export HP_RURL="https://hgdownload.soe.ucsc.edu/goldenPath/mm39/bigZips/mm39.fa.gz"

###############################################################

export HP_E=300 # extension(circularization)

################################################################

GENOME REFERENCES

export HP_O=Human # organism: Human, Mouse... export HP_MT=chrM # chrM, rCRS or RSRS, FASTA file available under $HP_RDIR export HP_MTC=chrMC export HP_MTR=chrMR export HP_MTLEN=16569 export HP_NUMT=NUMT # NUMT FASTA file under $HP_RDIR

Mouse

export HP_O=Mouse # organism: Human, Mouse...

export HP_MT=chrM # chrM, rCRS or RSRS, FASTA file available under $HP_RDIR

export HP_MTC=chrMC

export HP_MTR=chrMR

export HP_MTLEN=16299

export HP_NUMT=NUMT # NUMT FASTA file under $HP_RDIR

################################################################

OTHER

export HP_CN=1 # do compute mtDNA copy number export HP_L=222000 # number of MT reads to subsample; empty: no subsampling; 222000 150bp reads => ~2000x MT coverage export HP_FOPT="-q 15 -e 0" # FASTP options: Ex: " -q 20 -e 30 "; -q: min base quality; -e: avg quality thold export HP_DOPT="--removeDups" # samblaster option; leave empty if no deduplication should be done export HP_GOPT="--mitochondria-mode" # gatk mutect2 additional options : Ex "-max-reads-per-alignment-start 50" , "--mitochondria-mode"

export HP_M=mutect2 # SNV caller: mutect2,mutserve or freebayes export HP_I=2 # number of SNV iterations : 0,1,2

0: compute read counts,mtDNA-CN

                             #  1:1 iteration (mutect2,mutserve)
                             #  2:2 iterations (mutect2)

export HP_T1=03 # heteroplasmy tholds export HP_T2=05 export HP_T3=10

export HP_DP=50 # minimum coverage; lowered from 100 to 50 export HP_V=gridss # SV caller: gridssexport HP_DP=100 # minimum coverage: Ex 100

export HP_FRULE="perl -ane 'print unless(/strict_strand|strand_bias|base_qual|map_qual|weak_evidence|slippage|position|Homopolymer/ and /:0.[01234]\d+$/);' | bcftools filter -e 'DP<$HP_DP'" # filter rule (or just "tee")

export HP_P=10 # number of processors export HP_MM="30G" # maximum memory export HP_JOPT="-Xms$HP_MM -Xmx$HP_MM -XX:ParallelGCThreads=$HP_P" # JAVA options ################################################################

INPUT/OUTPUT

PWD=pwd -P

export HP_FDIR=$PWD/f # fastq input file directory ; .fq or .fq.gz file extension

export HP_ADIR=/mnt/data/WGS_sun/5.gatk/bqsr # bams or crams input file directory export HP_ODIR=$PWD/out/ # output dir export HP_IN=$PWD/in.txt # input file to be generated `

run.all.sh: `cat run.all.sh

!/usr/bin/env bash

set -eux

export HP_SDIR=/mnt/data/hbyuan/application/MitoHPC/MitoHPC/scripts export HP_BDIR=/mnt/data/hbyuan/application/MitoHPC/MitoHPC/bin/ export HP_JDIR=/mnt/data/hbyuan/application/MitoHPC/MitoHPC/java/ export HP_RDIR=/mnt/data/hbyuan/application/MitoHPC/MitoHPC/RefSeq/ export HP_ODIR=/mnt/data/hbyuan/project/WGS_TLS/out/ export HP_RNAME=hs38DH export HP_RMT=chrM export HP_RNUMT="chr1:629084-634672 chr17:22521208-22521639" export HP_RCOUNT=3366 export HP_RURL="ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa" export HP_O=Human export HP_MT=chrM export HP_MTR=chrMR export HP_MTC=chrMC export HP_MTLEN=16569 export HP_NUMT=NUMT export HP_CN=1 export HP_E=300 export HP_L=222000 export HP_IN=/mnt/data/hbyuan/project/WGS_TLS/in.txt export HP_ODIR=/mnt/data/hbyuan/project/WGS_TLS/out/ export HP_M=mutect2 export HP_I=2 export HP_T1=03 export HP_T2=05 export HP_T3=10 export HP_V=gridss export HP_DP=50 export HP_FRULE="perl -ane 'print unless(/strict_strand|strand_bias|base_qual|map_qual|weak_evidence|slippage|position|Homopolymer/ and /:0.[01234]\d+$/);' | bcftools filter -e 'DP<50'" export HP_FOPT="-q 15 -e 0" export HP_DOPT="--removeDups" export HP_GOPT="--mitochondria-mode" export HP_JOPT="-Xms30G -Xmx30G -XX:ParallelGCThreads=10" export HP_MM="30G" export HP_P=10 export PATH=/mnt/data/hbyuan/application/MitoHPC/MitoHPC/scripts:/mnt/data/hbyuan/application/MitoHPC/MitoHPC/bin/:$PATH

bash /mnt/data/hbyuan/application/MitoHPC/MitoHPC/scripts/filter.sh hanwei_bqsr /mnt/data/UKB/WGS_sun/5.gatk/bqsr/hanwei_bqsr.bam /mnt/data/hbyuan/project/WGS_TLS/out//hanwei_bqsr/hanwei_bqsr

bash /mnt/data/hbyuan/application/MitoHPC/MitoHPC/scripts/getSummary.sh /mnt/data/hbyuan/project/WGS_TLS/out/ `

dpuiu commented 2 days ago

Hi Yuan,

Sorry to hear about the difficulties in running MitoHPC.

Your run.all.sh script seems fine. You only have one sample (hanwei_bqsr) that you are trying to analyze, correct ?

Just make sure that the bam was generated by aligning the reads in paired-end mode using one of the standard aligners (bwa,bowtie2,minimap2...), sorted and indexed .

Then run "bash ./run.all.sh"

If you could send me the output of the "ls -l -t1 out/ " command, I might get a better idea what the issue is.