merge iterations into extra params

Author: pdimens

config.yaml

@@ -1,14 +1,14 @@
 # Configuration file for LepWrap
 
-###################
-#### Lep-Map 3 ####
-###################
+#########################################################
+#                       Lep-Map 3                       #
+#########################################################
 
 #---------------#
 #  ParentCall2  #
 #---------------#
 # the filtered VCF file with your genotype likelihoods:
-vcf: "out.15.miss95recode.vcf"
+vcf: "out.15.missrecode90.vcf"
 
 # Instructions to create pedigree file: https://sourceforge.net/p/lep-map3/wiki/software/LepMap3 Home/#parentcall2
 # the pedigree file associated with your data
@@ -35,16 +35,16 @@ extra_params_Filtering: ""
 # LOD score in the range of lod_min to lod_max
 
 # The minimum LOD for SeperateChromosomes2
-lod_min: 20
+lod_min: 10
 
 # The maximum LOD for SeperateChromosomes2
-lod_max: 40
+lod_max: 50
 
 # Use only markers with informative father (1), mother(2), both parents(3) or neither parent(0)
-informative: "informativeMask=123"
+informative: "informativeMask=3"
 
 # If there are any additional parameters you would like to use for SeparateChromosomes2 (e.g. distrotionLOD=1), add them here
-extra_params_SeparateChromosomes: "sizeLimit=5 distortionLOD=1"
+extra_params_SeparateChromosomes: "sizeLimit=5 distortionLod=1"
 
 
 #-------------------#
@@ -61,7 +61,7 @@ lod_limit: "lodLimit=2"
 lod_difference: "lodDifference=2"
 
 # If there are any additional parameters you would like to use for JoinSingles2All (e.g. iterate=1), add them here
-extra_params_JoinSingles: "iterate=1 distortionLOD=1"
+extra_params_JoinSingles: "iterate=1 distortionLod=1"
 
 
 #-----------------#
@@ -72,11 +72,9 @@ extra_params_JoinSingles: "iterate=1 distortionLOD=1"
 # Set exp_lg to your expected number of chromosomes
 exp_lg: 24
 
-# Set iterations to the number of iterations you want per chromosome (more is better). Recommend 30 or more
-iterations: 100
-
 # If there are any additional parameters you would like to use for OrderMarkers2 (e.g. hyperPhaser=1), add them here
-extra_params_OrderMarkers: "hyperPhaser=1 useKosambi=1 phasingIterations=2"
+# I recommend setting numMergeIterations to ~100 (Lep-Map3 default is 6)
+extra_params_OrderMarkers: "hyperPhaser=1 useKosambi=1 phasingIterations=2 numMergeIterations=100"
 
 
 #-----------------#
@@ -99,25 +97,25 @@ trim_cutoff: 100
 #--------------------#
 # The second round of OrderMarkers will use the same basic parameters as the first round (but not the extra params)
 # If there are additional parameters you would like to use, add them here:
-extra_params_reOrderMarkers: "improveOrder=1 useKosambi=1"
-
+extra_params_reOrderMarkers: "improveOrder=1 useKosambi=1 numMergeIterations=75"
 
 #-----------------------#
 #  Calculate Distances  #
 #-----------------------#
 # If you used useKosambi=1 or useMorgan=1 for Ordering/reOrdering, add that same
 # parameter to distance_method: 
-distance_method: ""
+distance_method: "useKosambi=1"
+
 
+#########################################################
+#                       Lep-Anchor                      #
+#########################################################
 
-####################
-#### Lep-Anchor ####
-####################
 # change this to true if you also want to run Lep-Anchor
 run_lepanchor: true
 
 # the path to the genome assembly you are trying to anchor
-# ideally it *is not* gzipped and ends in .fa, but there are built-in workarounds if it is. 
+# ideally it is *not* gzipped and ends in .fa, but there are built-in workarounds if it is. 
 assembly: "YFT.genome.latest.fasta"
 
 # the number of linkage groups you have
@@ -182,4 +180,4 @@ extra_params_PlaceOrient: "keepEmptyIntervals=1 numRuns=10"
 LA_edge_length: 20
 
 # Set trim_cuttoff to the centiMorgan distance cutoff (5 is reasonable)
-LA_trim_cutoff: 3
+LA_trim_cutoff: 5

rules/LepMap3/LepMap3.smk

@@ -13,7 +13,7 @@ data_tol=config["data_tol"]
 filtering_extra = config["extra_params_Filtering"]
 # separate chromosomes #
 lod_max = str(config["lod_max"])
-lod_range = list(range(config["lod_min"], config["lod_max"]+1))
+lod_range = list(range(config["lod_min"], lod_max+1))
 informative = config["informative"]
 sepchrom_extra = config["extra_params_SeparateChromosomes"]
 # join singles #
@@ -22,16 +22,14 @@ lod_lim = config["lod_limit"]
 lod_diff = config["lod_difference"]
 js2a_extra = config["extra_params_JoinSingles"]
 # ordering #
-lg_range = list(range(1, config["exp_lg"]+1))
 lg_count = config["exp_lg"]
-ITER = config["iterations"]
+lg_range = list(range(1, lg_count+1))
 order_extra = config["extra_params_OrderMarkers"]
 # trimming #
 edge_len = str(config["edge_length"])
 trim_thresh = str(config["trim_cutoff"])
 # ordering II #
 reorder_extra = config["extra_params_reOrderMarkers"]
-ITER2 = round(ITER/2)
 # distances #
 dist_method = config["distance_method"]
 

rules/LepMap3/order.smk

@@ -11,12 +11,11 @@ rule order_markers:
     message: "Ordering linkage group {params.chrom} with {params.iterations} iterations"
     params:
         chrom = "{lg_range}",
-        iterations = ITER,
         extra = order_extra
     threads: 2
     shell:
         """
-        zcat {input.datacall} | java -cp software/LepMap3 OrderMarkers2 map={input.filt_map} {params.extra} data=- numThreads={threads} numMergeIterations={params.iterations} chromosome={params.chrom} &> {output.runlog}
+        zcat {input.datacall} | java -cp software/LepMap3 OrderMarkers2 map={input.filt_map} {params.extra} data=- numThreads={threads} chromosome={params.chrom} &> {output.runlog}
         sed -n '/\*\*\* LG \=/,$p' {output.runlog} > {output.lg}
         grep "recombin" {output.runlog} > {log.recomb}
         awk '/#java/{{flag=1}} flag; /logL/{{flag=0}}' {output.runlog} > {log.run}

rules/LepMap3/reorder.smk

@@ -12,12 +12,11 @@ rule reorder_markers:
     message: "Reordering linkage group {params.lg} with {params.iterations} iterations"
     params:
         lg = "{lg_range}",
-        iterations = ITER2,
         extra = reorder_extra
     threads: 2
     shell:
         """
-        zcat {input.datacall} | java -cp software/LepMap3 OrderMarkers2 {params.extra} map={input.filt_map} data=- numThreads={threads} evaluateOrder={input.lg_order} numMergeIterations={params.iterations} &> {output.runlog}
+        zcat {input.datacall} | java -cp software/LepMap3 OrderMarkers2 {params.extra} map={input.filt_map} data=- numThreads={threads} evaluateOrder={input.lg_order} &> {output.runlog}
         sed -n '/\*\*\* LG \=/,$p' {output.runlog} > {output.lg}
         grep "recombin" {output.runlog} > {log.recomb}
         awk '/#java/{{flag=1}} flag; /logL/{{flag=0}}' {output.runlog} > {log.run}