Merge pull request #103 from Cristianetaniguti/develop

Develop

.dockerfiles/reads2map/Dockerfile

@@ -16,7 +16,7 @@ RUN apt update \
    libxml2-dev \
    libnlopt-dev=2.4.2+dfsg-4 \
    devscripts=2.17.12ubuntu1.1 \
-   software-properties-common=0.96.24.32.18 \
+   software-properties-common=0.96.24.32.22 \
    libgit2-dev=0.26.0+dfsg.1-1.1ubuntu0.2 \
    && pip install matplotlib==2.0.2 numpy 
 
@@ -53,7 +53,7 @@ RUN Rscript -e 'remotes::install_version("gsalib",upgrade="never", version = "2.
 RUN Rscript -e 'remotes::install_github("tpbilton/GUSbase", ref = "92119b9c57faa7abeede8236d24a4a8e85fb3df7")'
 
 RUN Rscript -e 'remotes::install_github("tpbilton/GUSMap", ref = "4d7d4057049819d045750d760a45976c8f30dac6")'
-RUN Rscript -e 'remotes::install_github("dcgerard/updog")'
+RUN Rscript -e 'remotes::install_github("dcgerard/updog", ref="f1")'
 
 RUN Rscript -e 'remotes::install_github("Cristianetaniguti/onemap")'
 

.dockerfiles/stacks/Dockerfile

@@ -2,9 +2,9 @@ FROM ubuntu:18.04
 
 RUN apt update \
     && apt install -y  build-essential wget zlib1g-dev   \
-    && wget http://catchenlab.life.illinois.edu/stacks/source/stacks-2.62.tar.gz \
-    && tar xfvz stacks-2.62.tar.gz \
-    && cd stacks-2.62 \
+    && wget http://catchenlab.life.illinois.edu/stacks/source/stacks-2.64.tar.gz \
+    && tar xfvz stacks-2.64.tar.gz \
+    && cd stacks-2.64 \
     && ./configure \
     && make \
     && make install

pipelines/EmpiricalMaps/EmpiricalMaps.changelog.md

@@ -1,3 +1,12 @@
+# 1.2.5
+
+* more flexibility to choose the probability to be used in the HMM: 
+
+* new parameters:
+- global_errors: array with global errors to be tested
+- genoprob_error: boolean defining if software probabilities should be tested
+- genoprob_global_errors: array with global errors to be tested together with the software probabilities following: 1 - (1 - global error) x (1 - software error probability)
+
 # 1.2.4
 
 * runtimes adapted to run with Caper

pipelines/EmpiricalMaps/EmpiricalMaps.inputs.json

@@ -7,17 +7,24 @@
     "cross": "String",
     "multiallelics": "Boolean"
   },
+  "Maps.run_supermassa": "Boolean (optional, default = false)",
   "Maps.max_cores": "Int",
   "Maps.gatk_vcf_multi": "File? (optional)",
   "Maps.gatk_mchap": "String",
   "Maps.vcfs_counts_source": "Array[String]",
-  "Maps.filters": "String? (optional)",
-  "Maps.filt_segr": "String? (optional)",
-  "Maps.prob_thres": "Float? (optional)",
+  "Maps.run_gusmap": "Boolean (optional, default = false)",
+  "Maps.genoprob_error": "Boolean (optional, default = true)",
+  "Maps.prob_thres": "Float (optional, default = 0.8)",
   "Maps.ploidy": "Int",
   "Maps.vcfs_software": "Array[String]",
   "Maps.filter_noninfo": "Boolean",
+  "Maps.filters": "String? (optional)",
+  "Maps.global_errors": "Array[String] (optional, default = [\"0.05\"])",
+  "Maps.run_polyrad": "Boolean (optional, default = true)",
   "Maps.vcfs": "Array[File]",
-  "Maps.replaceADbyMissing": "String"
+  "Maps.replaceADbyMissing": "String",
+  "Maps.genoprob_global_errors": "Array[String] (optional, default = [\"0.05\"])",
+  "Maps.filt_segr": "String? (optional)",
+  "Maps.run_updog": "Boolean (optional, default = true)"
 }
 

pipelines/EmpiricalMaps/EmpiricalMaps.wdl

@@ -25,14 +25,17 @@ workflow Maps {
         Boolean run_polyrad = true
         Boolean run_gusmap = false
         Boolean filter_noninfo
-        String replaceADbyMissing
+        Boolean replaceADbyMissing 
         File? gatk_vcf_multi
         String gatk_mchap
         String? filters
         Int max_cores
         Int ploidy
-        Float? prob_thres
+        Float prob_thres = 0.8
         String? filt_segr
+        Array[String] global_errors = ["0.05"]
+        Boolean genoprob_error = true
+        Array[String] genoprob_global_errors = ["0.05"]
     }
 
     if (defined(filters)) {
@@ -42,7 +45,8 @@ workflow Maps {
                 vcfs_SNPCall_software = vcfs_software,
                 vcfs_Counts_source = vcfs_counts_source,
                 vcfs_GenoCall_software = range(length(vcfs_software)),
-                filters = filters
+                filters = filters,
+                chromosome = dataset.chromosome
         }
     }
 
@@ -84,7 +88,11 @@ workflow Maps {
                         multiallelics = dataset.multiallelics,
                         multiallelics_file = splitgeno.multiallelics,
                         max_cores = max_cores,
-                        ploidy = ploidy
+                        ploidy = ploidy,
+                        global_errors = global_errors,
+                        genoprob_error = genoprob_error,
+                        prob_thres = prob_thres,
+                        genoprob_global_errors = genoprob_global_errors
                 }
             }
 
@@ -102,7 +110,11 @@ workflow Maps {
                         multiallelics = dataset.multiallelics,
                         multiallelics_file = splitgeno.multiallelics,
                         max_cores = max_cores,
-                        ploidy = ploidy
+                        ploidy = ploidy,
+                        global_errors = global_errors,
+                        genoprob_error = genoprob_error,
+                        prob_thres = prob_thres,
+                        genoprob_global_errors = genoprob_global_errors
                 }
             }
 
@@ -120,7 +132,11 @@ workflow Maps {
                         multiallelics = dataset.multiallelics,
                         multiallelics_file = splitgeno.multiallelics,
                         max_cores = max_cores,
-                        ploidy = ploidy
+                        ploidy = ploidy,
+                        global_errors = global_errors,
+                        genoprob_error = genoprob_error,
+                        prob_thres = prob_thres,
+                        genoprob_global_errors = genoprob_global_errors
                 }
             }
 
@@ -151,8 +167,10 @@ workflow Maps {
                         multiallelics = dataset.multiallelics,
                         max_cores = max_cores,
                         multiallelics_file = splitgeno.multiallelics,
-                        multiallelics_mchap = gatk_vcf_multi,
-                        mchap = gatk_mchap
+                        global_errors = global_errors,
+                        genoprob_error = genoprob_error,
+                        prob_thres = prob_thres,
+                        genoprob_global_errors = genoprob_global_errors
                 }
             }
         }

pipelines/EmpiricalReads2Map/EmpiricalReads2Map.changelog.md

@@ -1,3 +1,16 @@
+# 1.4.2
+
+* more flexibility to choose the probability to be used in the HMM: 
+
+* new parameters:
+- global_errors: array with global errors to be tested
+- genoprob_error: boolean defining if software probabilities should be tested
+- genoprob_global_errors: array with global errors to be tested together with the software probabilities following: 1 - (1 - global error) x (1 - software error probability)
+
+# 1.4.1
+
+* Use BCFtools norm to left-align indel marker positions identified by GATK, STACKs and freebayes
+
 # 1.4.0
 
 * STACKs included

pipelines/EmpiricalReads2Map/EmpiricalReads2Map.inputs.json

@@ -10,16 +10,10 @@
   "EmpiricalReads.chunk_size": "Int",
   "EmpiricalReads.Maps.run_gusmap": "Boolean (optional, default = false)",
   "EmpiricalReads.gatk_mchap": "Boolean (optional, default = false)",
+  "EmpiricalReads.Maps.global_errors": "Array[String] (optional, default = [\"0.05\"])",
   "EmpiricalReads.rm_dupli": "Boolean (optional, default = false)",
   "EmpiricalReads.SNPCalling.max_ram": "Int",
-  "EmpiricalReads.dataset": {
-    "parent2": "String",
-    "name": "String",
-    "parent1": "String",
-    "chromosome": "String",
-    "cross": "String",
-    "multiallelics": "Boolean"
-  },
+  "EmpiricalReads.Maps.genoprob_error": "Boolean (optional, default = true)",
   "EmpiricalReads.ploidy": "Int (optional, default = 2)",
   "EmpiricalReads.n_chrom": "Int",
   "EmpiricalReads.SNPCalling.chunk_size": "Int",
@@ -43,6 +37,15 @@
   "EmpiricalReads.Maps.filter_noninfo": "Boolean",
   "EmpiricalReads.SNPCalling.pop_map": "File? (optional)",
   "EmpiricalReads.samples_info": "File",
+  "EmpiricalReads.dataset": {
+    "parent2": "String",
+    "name": "String",
+    "parent1": "String",
+    "chromosome": "String",
+    "cross": "String",
+    "multiallelics": "Boolean"
+  },
+  "EmpiricalReads.Maps.genoprob_global_errors": "Array[String] (optional, default = [\"0.05\"])",
   "EmpiricalReads.pair_end": "Boolean (optional, default = false)",
   "EmpiricalReads.run_gatk": "Boolean (optional, default = true)"
 }

pipelines/EmpiricalSNPCalling/EmpiricalSNPCalling.changelog.md

@@ -1,3 +1,7 @@
+# 1.4.1
+
+*  Use BCFtools norm to left-align indel marker positions identified by GATK, STACKs and freebayes
+
 # 1.4.0
 
 * STACKs included
@@ -47,4 +51,4 @@ This workflow requires:
 
 * Diploid or polyploid specie
 * Single-end reads
-* A reference genome
+* A reference genome

pipelines/EmpiricalSNPCalling/EmpiricalSNPCalling.wdl

@@ -92,6 +92,7 @@ workflow SNPCalling {
   if(run_stacks) {
     call stacks.StacksGenotyping {
       input:
+        references = references,
         bams = CreateAlignmentFromFamilies.bam,
         pop_map = pop_map,
         max_cores = max_cores

pipelines/SimulatedReads2Map/SimulatedReads2Map.changelog.md

@@ -1,3 +1,12 @@
+# 1.0.2
+
+* more flexibility to choose the probability to be used in the HMM: 
+
+* new parameters:
+- global_errors: array with global errors to be tested
+- genoprob_error: boolean defining if software probabilities should be tested
+- genoprob_global_errors: array with global errors to be tested together with the software probabilities following: 1 - (1 - global error) x (1 - software error probability)
+
 # 1.0.1
 
 * runtimes adapted to run with Caper

pipelines/SimulatedReads2Map/SimulatedReads2Map.inputs.json

@@ -1,5 +1,7 @@
 {
+  "SimulatedReads.SimulatedSingleFamily.global_errors": "Array[String]",
   "SimulatedReads.filters": "String? (optional)",
+  "SimulatedReads.SimulatedSingleFamily.genoprob_global_errors": "Array[String]",
   "SimulatedReads.hardfilters": "Boolean (optional, default = true)",
   "SimulatedReads.family": {
     "cmBymb": "Float? (optional)",
@@ -22,6 +24,8 @@
   "SimulatedReads.n_chrom": "Int",
   "SimulatedReads.chunk_size": "Int (optional, default = 5)",
   "SimulatedReads.max_cores": "Int",
+  "SimulatedReads.SimulatedSingleFamily.genoprob_error": "String",
+  "SimulatedReads.SimulatedSingleFamily.prob_thres": "Float",
   "SimulatedReads.number_of_families": "Int",
   "SimulatedReads.gatk_mchap": "Boolean (optional, default = false)",
   "SimulatedReads.global_seed": "Int",

subworkflows/SimulatedSingleFamily.wdl

@@ -27,6 +27,10 @@ workflow SimulatedSingleFamily {
     Boolean gatk_mchap
     Boolean hardfilters
     Int n_chrom
+    Float prob_thres
+    Array[String] global_errors
+    String genoprob_error
+    Array[String] genoprob_global_errors
   }
 
   call simulation.CreateAlignmentFromSimulation {
@@ -148,7 +152,11 @@ workflow SimulatedSingleFamily {
             depth = sequencing.depth,
             multiallelics = sequencing.multiallelics,
             multiallelics_file = splitgeno.multiallelics,
-            ploidy = ploidy
+            ploidy = ploidy,
+            global_errors = global_errors,
+            genoprob_error = genoprob_error,
+            prob_thres = prob_thres,
+            genoprob_global_errors = genoprob_global_errors
         }
 
         call genotyping.onemapMaps as supermassaMaps {
@@ -168,7 +176,11 @@ workflow SimulatedSingleFamily {
             depth = sequencing.depth,
             multiallelics = sequencing.multiallelics,
             multiallelics_file = splitgeno.multiallelics,
-            ploidy = ploidy
+            ploidy = ploidy,
+            global_errors = global_errors,
+            genoprob_error = genoprob_error,
+            prob_thres = prob_thres,
+            genoprob_global_errors = genoprob_global_errors
         }
 
         call genotyping.onemapMaps as polyradMaps {
@@ -188,7 +200,11 @@ workflow SimulatedSingleFamily {
             depth = sequencing.depth,
             multiallelics = sequencing.multiallelics,
             multiallelics_file = splitgeno.multiallelics,
-            ploidy = ploidy
+            ploidy = ploidy,
+            global_errors = global_errors,
+            genoprob_error = genoprob_error,
+            prob_thres = prob_thres,
+            genoprob_global_errors = genoprob_global_errors
         }
     }
 

subworkflows/genotyping_empirical.wdl

@@ -18,6 +18,10 @@ workflow onemapMapsEmp {
     Boolean multiallelics
     File? multiallelics_file
     Int ploidy
+    Float prob_thres 
+    Array[String] global_errors
+    Boolean genoprob_error
+    Array[String] genoprob_global_errors
   }
 
   call utilsR.ReGenotyping {
@@ -74,31 +78,32 @@ workflow onemapMapsEmp {
       parent1 = parent1,
       parent2 = parent2,
       multiallelics = multiallelics,
-      SNPCall_program = SNPCall_program
+      SNPCall_program = SNPCall_program,
+      global_errors = global_errors,
+      genoprob_error = genoprob_error,
+      prob_thres = prob_thres,
+      genoprob_global_errors = genoprob_global_errors,
+      GenotypeCall_program = GenotypeCall_program
   }
 
-  Array[String] methods                         = [GenotypeCall_program, GenotypeCall_program + "0.05"]
-  Array[File] objects                           = [SetProbs.probs_onemap_obj, SetProbs.globalerror_onemap_obj]
-  Array[Pair[String, File]] methods_and_objects = zip(methods, objects)
-
-  scatter (item in methods_and_objects) {
+  scatter (item in range(length(SetProbs.probs_onemap_obj))) {
        call utilsR.CheckDepths {
            input:
-              onemap_obj = item.right,
+              onemap_obj = SetProbs.probs_onemap_obj[item],
               vcfR_obj = SetProbs.vcfR_obj,
               parent1 = parent1,
               parent2 = parent2,
               SNPCall_program = SNPCall_program,
-              GenotypeCall_program = item.left,
+              GenotypeCall_program = SetProbs.probs_onemap_obj_names[item],
               CountsFrom = CountsFrom,
               max_cores = max_cores
        }
 
        call utilsR.FiltersReportEmp {
             input:
-              onemap_obj = item.right,
+              onemap_obj = SetProbs.probs_onemap_obj[item],
               SNPCall_program = SNPCall_program,
-              GenotypeCall_program = item.left,
+              GenotypeCall_program = SetProbs.probs_onemap_obj_names[item],
               CountsFrom = CountsFrom,
               chromosome = chromosome
         }
@@ -107,7 +112,7 @@ workflow onemapMapsEmp {
           input:
             sequence_obj = FiltersReportEmp.onemap_obj_filtered,
             SNPCall_program = SNPCall_program,
-            GenotypeCall_program = item.left,
+            GenotypeCall_program = SetProbs.probs_onemap_obj_names[item],
             CountsFrom = CountsFrom,
             max_cores = max_cores
           }