From 48c88309fabd41092573732283fdbf103df622ec Mon Sep 17 00:00:00 2001
From: Zachary Nolen <zachary.nolen@biol.lu.se>
Date: Fri, 29 Mar 2024 15:45:45 +0100
Subject: [PATCH] Add in ability to remove transitions easily from config

---
 .test/config/config.yaml       | 1 +
 config/config.yaml             | 1 +
 docs/config.md                 | 4 ++++
 workflow/rules/3.1_safs.smk    | 7 +++++--
 workflow/rules/3.2_beagles.smk | 3 ++-
 workflow/rules/common.smk      | 7 +++++++
 6 files changed, 20 insertions(+), 3 deletions(-)

diff --git a/.test/config/config.yaml b/.test/config/config.yaml
index 0a11c57..bd998bc 100644
--- a/.test/config/config.yaml
+++ b/.test/config/config.yaml
@@ -135,6 +135,7 @@ params:
     gl_model: 2 # gl_model - 1=SAMtools, 2=GATK, 3=SOAPsnp, 4=SYK
     maxdepth: 1000
     missing_data_mindepthind: 1
+    rmtrans: true
     extra: ""
     extra_saf: ""
     extra_beagle: ""
diff --git a/config/config.yaml b/config/config.yaml
index e8d5229..88e68a2 100644
--- a/config/config.yaml
+++ b/config/config.yaml
@@ -134,6 +134,7 @@ params:
     gl_model: 2 # gl_model - 1=SAMtools, 2=GATK, 3=SOAPsnp, 4=SYK
     maxdepth: 1000
     missing_data_mindepthind: 1
+    rmtrans: false
     extra: ""
     extra_saf: ""
     extra_beagle: ""
diff --git a/docs/config.md b/docs/config.md
index 799113e..37e14f4 100644
--- a/docs/config.md
+++ b/docs/config.md
@@ -450,6 +450,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.)
     - `extra:` Additional options to pass to ANGSD during genotype likelihood
       calculation at all times. This is primarily useful for adding BAM input
       filters. Note that `--remove_bads` and `-only_proper_pairs` are enabled
diff --git a/workflow/rules/3.1_safs.smk b/workflow/rules/3.1_safs.smk
index 7bd9af7..3895e37 100644
--- a/workflow/rules/3.1_safs.smk
+++ b/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,
         out=lambda w, output: os.path.splitext(output.arg)[0],
     resources:
         runtime=lambda wildcards, attempt: attempt * 180,
@@ -53,7 +54,8 @@ rule angsd_doSaf_pop:
         angsd -doSaf 1 -bam {input.bam} -GL {params.gl_model} -ref {input.ref} \
             -nThreads {threads} {params.extra} -minMapQ {params.mapQ} \
             -minQ {params.baseQ} -sites {input.sites} -anc {input.anc} \
-            {params.extra_saf} -rf {input.regions} -out {params.out} &> {log}
+            -noTrans {params.trans} {params.extra_saf} -rf {input.regions} \
+            -out {params.out} &> {log}
         """
 
 
@@ -137,6 +139,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",
@@ -146,5 +149,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}
         """
diff --git a/workflow/rules/3.2_beagles.smk b/workflow/rules/3.2_beagles.smk
index 6f69532..83972e7 100644
--- a/workflow/rules/3.2_beagles.smk
+++ b/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)),
         out=lambda w, output: os.path.splitext(output.arg)[0],
     threads: lambda wildcards, attempt: attempt
@@ -50,7 +51,7 @@ 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.counts} -out {params.out} &> {log}
+            -rmTrans {params.trans} {params.counts} -out {params.out} &> {log}
         """
 
 
diff --git a/workflow/rules/common.smk b/workflow/rules/common.smk
index e2c4daa..0861788 100644
--- a/workflow/rules/common.smk
+++ b/workflow/rules/common.smk
@@ -623,6 +623,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