Merge pull request #39 from zjnolen/rmtrans

Add in ability to remove transitions easily from config

.test/config/config.yaml

@@ -137,6 +137,7 @@ params:
     mindepthind: 1
     minind_dataset: 0.75 # as a proportion of N
     minind_pop: 0.5 # as a proportion of N
+    rmtrans: true
     extra: "" # goes to all ANGSD runs
     extra_saf: "" # goes to all -doSaf runs
     extra_beagle: "" # goes to all -doGlf 2 runs

config/config.yaml

@@ -136,6 +136,7 @@ params:
     mindepthind: 1
     minind_dataset: 0.75 # as a proportion of N
     minind_pop: 0.75 # as a proportion of N
+    rmtrans: false
     extra: "" # goes to all ANGSD runs
     extra_saf: "" # goes to all -doSaf runs
     extra_beagle: "" # goes to all -doGlf 2 runs

docs/config.md

@@ -451,6 +451,10 @@ or a pull request and I'll gladly put it in.
     - `maxdepth:` When calculating individual depth, sites with depth higher
       than this will be binned to this value. Should be fine for most to leave
       at `1000`. (integer, [docs](http://www.popgen.dk/angsd/index.php/Depth))
+    - `rmtrans:` Removes transitions using ANGSD, effectively removing them
+      from downstream analyses. This is useful for removing DNA damage from
+      analyses, and will automatically set the appropriate ANGSD flags (i.e.
+      using `-noTrans 1` for SAF files and `-rmTrans 1` for Beagle files.)
     - `mindepthind:` Individuals with sequencing depth below this value at a
       position will be treated as having no data at that position by ANGSD.
       ANGSD defaults to 1 for this. Note that this can be separately set for

workflow/rules/3.1_safs.smk

@@ -44,6 +44,7 @@ rule angsd_doSaf_pop:
         extra_saf=config["params"]["angsd"]["extra_saf"],
         mapQ=config["mapQ"],
         baseQ=config["baseQ"],
+        trans=get_trans,
         minind=get_minind,
         mininddp=config["params"]["angsd"]["mindepthind"],
         out=lambda w, output: os.path.splitext(output.arg)[0],
@@ -55,8 +56,8 @@ rule angsd_doSaf_pop:
         angsd -doSaf 1 -bam {input.bam} -GL {params.gl_model} -ref {input.ref} \
             -nThreads {threads} {params.extra} {params.minind} -minMapQ {params.mapQ} \
             -minQ {params.baseQ} -sites {input.sites} -anc {input.anc} \
-            -setMinDepthInd {params.mininddp} {params.extra_saf} -rf {input.regions} \
-            -out {params.out} &> {log}
+            -noTrans {params.trans} {params.extra_saf} -rf {input.regions} \
+            -setMinDepthInd {params.mininddp} -out {params.out} &> {log}
         """
 
 
@@ -140,6 +141,7 @@ rule angsd_doSaf_sample:
         mindepthind=config["params"]["angsd"]["mindepthind_heterozygosity"],
         mapQ=config["mapQ"],
         baseQ=config["baseQ"],
+        trans=get_trans,
         out=lambda w, output: os.path.splitext(output.arg)[0],
     resources:
         runtime="120m",
@@ -149,5 +151,5 @@ rule angsd_doSaf_sample:
             -nThreads {threads} {params.extra} -minMapQ {params.mapQ} \
             -minQ {params.baseQ} -sites {input.sites} -anc {input.anc} \
             -setMinDepthInd {params.mindepthind} {params.extra_saf} \
-            -out {params.out}) &> {log}
+            -noTrans {params.trans} -out {params.out}) &> {log}
         """

workflow/rules/3.2_beagles.smk

@@ -38,6 +38,7 @@ rule angsd_doGlf2:
         minmaf=config["params"]["angsd"]["min_maf"],
         majmin=config["params"]["angsd"]["domajorminor"],
         counts=get_docounts,
+        trans=get_trans,
         nind=lambda w: len(get_samples_from_pop(w.population)),
         minind=get_minind,
         mininddp=config["params"]["angsd"]["mindepthind"],
@@ -52,7 +53,8 @@ rule angsd_doGlf2:
             -SNP_pval {params.snp_pval} -nThreads {threads} {params.extra} \
             -minMapQ {params.mapQ} -minQ {params.baseQ} -sites {input.sites} \
             -anc {input.anc} {params.extra_beagle} -rf {input.regions} {params.minind} \
-            -setMinDepthInd {params.mininddp} {params.counts} -out {params.out} &> {log}
+            -setMinDepthInd {params.mininddp} {params.counts} -rmTrans {params.trans} \
+            -out {params.out} &> {log}
         """
 
 

workflow/rules/common.smk

@@ -661,6 +661,13 @@ def get_docounts(wildcard):
     return ""
 
 
+# Determine whether transitions should be removed based on user configuration
+def get_trans(wildcards):
+    if config["params"]["angsd"]["rmtrans"]:
+        return 1
+    return 0
+
+
 # ngsLD