Merge branch 'dev' of https://github.com/ktmeaton/plague-phylogeography…

… into dev

workflow/rules/alignment.smk

@@ -57,7 +57,7 @@ rule eager:
             --max_memory {resources.mem_mb}.MB \
             --max_time {resources.time_min}m \
 			{config[eager_other]} \
-            -resume;"
+            -resume; > {output.log_txt}"
         "{scripts_dir}/eager_cleanup.sh {results_dir} {wildcards.reads_origin} {wildcards.sample}; "
 
 # -----------------------------------------------------------------------------#

workflow/scripts/filter_sites.py

@@ -112,6 +112,7 @@
 num_samples = len(alignment_in[:, 1])
 logging.info("Number of samples: " + str(num_samples))
 # Counter to store number of singleton sites
+constant_sites = 0
 biallelic_singleton_sites = 0
 multiallelic_singleton_sites = 0
 pseudo_singleton_sites = 0
@@ -169,27 +170,30 @@
     # store the variant type (default singleton)
     site_type = "singleton"
 
-    # Remove singleton columns sites
+    # Determine Site Type
+    # Constant site format [0,0,0,x]
+    if count_list.count(0) == 3:
+        constant_sites += 1
     # Biallelic format [0,0,1,x]
-    if count_list.count(0) == 2 and count_list.count(1) == 1:
+    elif count_list.count(0) == 2 and count_list.count(1) >= 1:
         biallelic_singleton_sites += 1
-        if not keep_singleton:
-            continue
     # Multi allelic format [0,1,1,x]
-    elif count_list.count(0) == 1 and count_list.count(1) == 2:
+    elif count_list.count(0) == 1 and count_list.count(1) >= 2:
         multiallelic_singleton_sites += 1
-        if not keep_singleton:
-            continue
-    # Pseudo allelic format [0,1,x,x] or [1,1,x,x]
+    # Pseudo allelic format [0,1,x,x] or [1,1,x,x] or [1,x,x,x]
     elif count_list.count(1) > 0:
         pseudo_singleton_sites += 1
-        if not keep_singleton:
-            continue
     else:
         parsimony_informative_sites += 1
         site_type = "parsimony"
 
-    # Decide if this
+    # Decide whether to skip over this site
+    if site_type == "constant":
+        continue
+    if site_type == "singleton" and not keep_singleton:
+        continue
+
+    # Decide whether this site has enough informative nucleotides
     num_data = count_g + count_t + count_a + count_c
     site_prop = num_data / num_samples
     # Check if the amount of missing data passes user parameter
@@ -226,6 +230,13 @@
 total_singleton_sites = (
     biallelic_singleton_sites + multiallelic_singleton_sites + pseudo_singleton_sites
 )
+logging.info(
+    "Constant sites: "
+    + str(constant_sites)
+    + " ("
+    + str(int(constant_sites / alignment_in_len * 100))
+    + "%)"
+)
 logging.info(
     "Parsimony informative sites: "
     + str(parsimony_informative_sites)