External peaks

src/AOS/atac.py

@@ -118,6 +118,13 @@
     show_default=True,
     help='Pseudocount to add to the count matrix prior to differential calling.'
 )
+@click.option(
+    '--peakset',
+    required=False,
+    default=None,
+    show_default=True,
+    help='Include an external peak file (bed format).'
+)
 def main(bamdir,
         outputdir,
         gtf,
@@ -133,7 +140,8 @@ def main(bamdir,
         upstreamuro,
         downstreamuro,
         featureuro,
-        pseudocount
+        pseudocount,
+        peakset
         ):
     # Init
     pf = Preflight(**locals())

src/AOS/helper.py

@@ -4,6 +4,7 @@
 import seaborn as sns
 import matplotlib.pyplot as plt
 from itertools import repeat
+import shutil
 
 def idx_to_mit(_f):
     count = 0
@@ -234,39 +235,42 @@ def tsv_to_bed(tsv, bed, grint):
             index=False
         )
 
-def peak_boundaries(peaks, genomefa, of):
-    chromdic = {}
-    with open(genomefa) as f:
-        for line in f:
-            if line.startswith('>'):
-                header = str(line.strip().replace('>', '').split(' ')[0])
-                chromdic[header] = 0
-            else:
-                chromdic[header] += len(line.strip())
-    bedlis = []
-    peakchange = 0 
-    with open(peaks) as f:
-        for line in f:
-            chrom = str(line.strip().split()[0])
-            start = int(line.strip().split()[1])
-            end = int(line.strip().split()[2])
-            if end > chromdic[chrom]:
-                bedlis.append(
-                    [chrom, start, chromdic[chrom]]
-                )
-                peakchange += 1
-            else:
-                bedlis.append(
-                    [chrom, start, end]
-                )
-    print("Changed {} peaks.".format(peakchange))
-    beddf = pd.DataFrame(bedlis)
-    beddf.to_csv(
-        of,
-        sep='\t',
-        header=False,
-        index=False
-    )
+def peak_boundaries(peaks, genomefa, peakset, of):
+    if not peakset:
+        chromdic = {}
+        with open(genomefa) as f:
+            for line in f:
+                if line.startswith('>'):
+                    header = str(line.strip().replace('>', '').split(' ')[0])
+                    chromdic[header] = 0
+                else:
+                    chromdic[header] += len(line.strip())
+        bedlis = []
+        peakchange = 0 
+        with open(peaks) as f:
+            for line in f:
+                chrom = str(line.strip().split()[0])
+                start = int(line.strip().split()[1])
+                end = int(line.strip().split()[2])
+                if end > chromdic[chrom]:
+                    bedlis.append(
+                        [chrom, start, chromdic[chrom]]
+                    )
+                    peakchange += 1
+                else:
+                    bedlis.append(
+                        [chrom, start, end]
+                    )
+        print("Changed {} peaks.".format(peakchange))
+        beddf = pd.DataFrame(bedlis)
+        beddf.to_csv(
+            of,
+            sep='\t',
+            header=False,
+            index=False
+        )
+    else:
+        shutil.copyfile(peakset, of)
 
 def PCA_colors(config):
     colors = [

src/AOS/preflight.py

@@ -27,7 +27,8 @@ def __init__(
         upstreamuro,
         downstreamuro,
         featureuro,
-        pseudocount
+        pseudocount,
+        peakset
     ):
         def retabspath(_p):
             if _p:
@@ -50,7 +51,8 @@ def retabspath(_p):
             'fna': retabspath(genomefasta),
             'motif': retabspath(motifs),
             'samplesheet': retabspath(samplesheet),
-            'comparison': retabspath(comparison)
+            'comparison': retabspath(comparison),
+            'peakset': retabspath(peakset)
         }
         self.vars = {
             'fragsize': fragsize,

src/AOS/rules/peaks.smk

@@ -103,10 +103,11 @@ rule peakbounds:
   output:
     'peakset/peaks.bed'
   params:
-    config['files']['fna']
+    fna = config['files']['fna'],
+    peakset = config['files']['peakset']
   threads: 1
   run:
-    peak_boundaries(input[0], params[0], output[0])
+    peak_boundaries(input[0], params.fna, params.peakset, output[0])
 
 rule uropa:
   input: