add cli option for 'jump' for hmmr_em .

Author: taoliu

MACS3/Commands/hmmratac_cmd.py

@@ -187,7 +187,7 @@ def run(args):
         options.info("#2 Use EM algorithm to estimate means and stddevs of fragment lengths")
         options.info("#  for mono-, di-, and tri-nucleosomal signals...")
         em_trainer = HMMR_EM(petrack, options.em_means[1:4],
-                             options.em_stddevs[1:4], seed=options.hmm_randomSeed)
+                             options.em_stddevs[1:4], seed=options.hmm_randomSeed, jump=options.em_jump)
         # the mean and stddev after EM training
         em_means = [options.em_means[0],]
         em_means.extend(em_trainer.fragMeans)

MACS3/Signal/HMMR_EM.py

@@ -103,7 +103,7 @@ class HMMR_EM:
     max_fraglen: cython.int
     epsilon: cython.float  # maximum difference to call the value is converged
     maxIter: cython.int  # maximum iteration
-    jump: cython.float
+    jump: cython.float # to amplify the difference for faster convergence
     seed: cython.int            # random seed for downsampling
     converged: cython.bint      # wheter the EM is converged
     sample_percentage: cython.float
@@ -118,7 +118,7 @@ def __init__(self, petrack: PETrackI, init_means: cython.list,
                  sample_percentage: cython.float = 10,
                  epsilon: cython.float = 0.05,
                  maxIter: cython.int = 20,
-                 jump: cython.float = 1.5, seed: cython.int = 12345):
+                 jump: cython.float = 0.5, seed: cython.int = 12345):
         """Initialize HMMR_EM object. The first three parameters are required.
 
         parameters:

bin/macs3

@@ -937,6 +937,10 @@ plus an extra option for the HMM model file like `macs3 hmmratac
     group_em.add_argument("--min-frag-p", dest="min_frag_p", type=float,
                           help="We will exclude the abnormal fragments that can't be assigned to any of the four signal tracks. After we use EM to find the means and stddevs of the four distributions, we will calculate the likelihood that a given fragment length fit any of the four using normal distribution. If a fragment length is not like ('not like' means the probability is lower than MIN_FRAG_P) any of short, mono, di, or tri-nuc fragment, we will exclude it while generating the four signal tracks for later HMM training and prediction. The value should be between 0 and 1. Smaller the value, more 'abnormal' fragments will be allowed. So if you want to include more fragments, make this value smaller. Default=0.001",
                           default=0.001)
+    # To amplify the difference for faster convergence
+    group_em.add_argument("--jump", dest="em_jump", type=float,
+                          help="This option will control the factor to amplify the difference between the observed and expected mean and variance values of the four distributions, then decide the expected mean and variance in the next EM step. Larger jump will lead to faster convergence but may be less stable. Default=0.5",
+                          default=0.5)
     # group for HMM
     group_hmm = argparser_hmmratac.add_argument_group("Hidden Markov Model arguments")
     #group_hmm.add_argument("-s", "--states", dest="hmm_states", type=int,

test/cmdlinetest

@@ -227,28 +227,28 @@ mkdir ${OUTPUTDIR_PREFIX}_run_hmmratac
 # macs3 hmmratac -i $ATACSEQBAM -n hmmratac_yeast500k --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac.log
 
 # check both gaussian and poisson hmm models, save everything
-macs3 hmmratac -i $ATACSEQBAM --hmm-type gaussian -n hmmratac_yeast500k --save-digested --save-states --save-likelihoods --save-training-data --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac.log
-macs3 hmmratac -i $ATACSEQBAM --hmm-type poisson -n hmmratac_yeast500k_poisson --save-digested --save-states --save-likelihoods --save-training-data --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac_poisson.log
+macs3 hmmratac -i $ATACSEQBAM --hmm-type gaussian -n hmmratac_yeast500k --jump 1.5 --save-digested --save-states --save-likelihoods --save-training-data --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac.log
+macs3 hmmratac -i $ATACSEQBAM --hmm-type poisson -n hmmratac_yeast500k_poisson --jump 1.5 --save-digested --save-states --save-likelihoods --save-training-data --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac_poisson.log
 
 # check bedpe for both gaussian, poisson
-macs3 hmmratac -i $ATACSEQBED -n hmmratac_yeast500k_bedpe -f BEDPE --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac_bedpe.log
-macs3 hmmratac -i $ATACSEQBED --hmm-type poisson -n hmmratac_yeast500k_bedpe_poisson -f BEDPE --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac_bedpe_poisson.log
+macs3 hmmratac -i $ATACSEQBED -n hmmratac_yeast500k_bedpe -f BEDPE --jump 1.5 --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac_bedpe.log
+macs3 hmmratac -i $ATACSEQBED --hmm-type poisson -n hmmratac_yeast500k_bedpe_poisson -f BEDPE --jump 1.5 --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac_bedpe_poisson.log
 
 TRAINING_REGIONS_BED=${OUTPUTDIR_PREFIX}_run_hmmratac/hmmratac_yeast500k_training_regions.bed
 HMM_MODEL=${OUTPUTDIR_PREFIX}_run_hmmratac/hmmratac_yeast500k_model.json
 TRAINING_REGIONS_BED_POISSON=${OUTPUTDIR_PREFIX}_run_hmmratac/hmmratac_yeast500k_poisson_training_regions.bed
 HMM_MODEL_POISSON=${OUTPUTDIR_PREFIX}_run_hmmratac/hmmratac_yeast500k_poisson_model.json
 
 echo "15.1 hmmratac load training regions from bedfile"
-macs3 hmmratac -i $ATACSEQBAM -n hmmratac_yeast500k_load_training_regions -t ${TRAINING_REGIONS_BED} --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac_load_training_regions.log
-macs3 hmmratac -i $ATACSEQBAM --hmm-type poisson -n hmmratac_yeast500k_poisson_load_training_regions -t ${TRAINING_REGIONS_BED_POISSON} --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac_load_training_regions_poisson.log
+macs3 hmmratac -i $ATACSEQBAM --jump 1.5 -n hmmratac_yeast500k_load_training_regions -t ${TRAINING_REGIONS_BED} --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac_load_training_regions.log
+macs3 hmmratac -i $ATACSEQBAM --jump 1.5 --hmm-type poisson -n hmmratac_yeast500k_poisson_load_training_regions -t ${TRAINING_REGIONS_BED_POISSON} --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac_load_training_regions_poisson.log
 
 echo "15.2 hmmratac load hmm model file"
-macs3 hmmratac -i $ATACSEQBAM -n hmmratac_yeast500k_load_hmm_model --model ${HMM_MODEL} --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac_load_hmm_model.log
-macs3 hmmratac -i $ATACSEQBAM --hmm-type poisson -n hmmratac_yeast500k_poisson_load_hmm_model --model ${HMM_MODEL_POISSON} --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac_load_hmm_model_poisson.log
+macs3 hmmratac -i $ATACSEQBAM --jump 1.5 -n hmmratac_yeast500k_load_hmm_model --model ${HMM_MODEL} --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac_load_hmm_model.log
+macs3 hmmratac -i $ATACSEQBAM --jump 1.5 --hmm-type poisson -n hmmratac_yeast500k_poisson_load_hmm_model --model ${HMM_MODEL_POISSON} --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac_load_hmm_model_poisson.log
 
 echo "15.3 hmmratac scATAC test"
-macs3 hmmratac -i $FRAGFILE -f FRAG --barcodes $BARCODESFILE -n hmmratac_scatac_test --hmm-type poisson -n hmmratac_scatac_test --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac_scatac.log
+macs3 hmmratac -i $FRAGFILE -f FRAG --barcodes $BARCODESFILE -n hmmratac_scatac_test --hmm-type poisson --jump 1.5 --outdir ${OUTPUTDIR_PREFIX}_run_hmmratac &> ${OUTPUTDIR_PREFIX}_run_hmmratac/run_hmmratac_scatac.log
 
 echo "16. search for errors in log files"
 flag=0