remove markers in high ld with chr4 allele

quinlan-lab · Oct 25, 2023 · 49a4fbb · 49a4fbb
1 parent 8217e9d
commit 49a4fbb
Show file tree

Hide file tree

Showing 2 changed files with 9 additions and 3 deletions.
diff --git a/ihd/run_ihd_scan.py b/ihd/run_ihd_scan.py
@@ -108,10 +108,10 @@ def main(args):
     ld = squareform(r ** 2)
     np.fill_diagonal(ld, 1.)
 
-    # if specified, remove all markers with r2 > 0.5 with selected markers
+    # if specified, remove all markers with r2 > X with selected markers
     if args.adj_marker is not None:
         marker_idxs = geno_asint[geno_asint["marker"] == args.adj_marker].index.values
-        _, ld_markers = np.where(ld[marker_idxs] >= 0.5)
+        _, ld_markers = np.where(ld[marker_idxs] >= 0.1)
         geno_asint = geno_asint.iloc[geno_asint.index.difference(ld_markers)]
 
     print("Using {} genotypes that meet filtering criteria.".format(geno_asint.shape[0]))
@@ -190,6 +190,7 @@ def main(args):
         if chrom not in conf_int_chroms: continue
 
         chrom_genotype_matrix = chrom_df[samples].values
+        print (chrom_genotype_matrix.shape)
         # compute confidence intervals on the chromosome
         conf_int_lo, conf_int_hi, peak_markers = calculate_confint(
             spectra,

diff --git a/scripts/rules/run_ihd.smk b/scripts/rules/run_ihd.smk
@@ -1,10 +1,14 @@
+condition2markers = {"N": None, "D": "rs27509845", "B": "rs27509845"}
+
+
 rule run_ihd:
     input:
         singletons = "data/mutations/{cross}/annotated_filtered_singletons.condition_on_{condition}.csv",
         geno = "data/genotypes/{cross}.geno",
         config = "data/json/{cross}.json",
         py_script = "ihd/run_ihd_scan.py"
     output: "csv/{cross}.k{k}.genome.condition_on_{condition}.results.csv"
+    params: adj_marker = lambda wc: condition2markers[wc.condition]
     shell:
         """
         python {input.py_script} --mutations {input.singletons} \
@@ -14,7 +18,8 @@ rule run_ihd:
                                  -distance_method cosine \
                                  -permutations 10000 \
                                  -stratify_column true_epoch \
-                                 -threads 1 \
+                                 -threads 4 \
+                                 -adj_marker {params.adj_marker}
         """
 
 rule plot_ihd: