more tweaks

rules/LepAnchor.smk

@@ -357,13 +357,13 @@ rule mareymap_data:
     agp = expand("10_Anchoring/agp/chr.{lgs}.agp", lgs = lg_range)
   output: 
     mareydata = "11_MareyMaps/marey.data.gz",
-    midpoints = "11_MareyMaps/marey.midpoints_data.gz"
+    sexavg = "11_MareyMaps/marey.data.sexavg.gz"
   log: "11_MareyMaps/missing_scaffolds.txt"
   message: 
     """
     Creating Marey map interval data
     first points in uncertainty intervals  | {output.mareydata}
-    midpoints in uncertainty intervals     | {output.midpoints}  
+    midpoints in uncertainty intervals     | {output.sexavg}  
     """
   params:
     chrom = lg
@@ -377,20 +377,24 @@ rule mareymap_data:
     for c in $(seq 1 {params.chrom})
     do
       awk -vn=$c '($3==n)' {input.lift} | awk -f software/LepAnchor/scripts/liftover.awk 10_Anchoring/agp/chr.$c.agp - | awk -vm=1 '(/LG/ && NR>=4){{if (NF>4) s=0.5; else s=1;print $1"\t"$2"\t"$3"\t"m"\t"s*($4+$5)}}' | gzip
-    done > {output.midpoints} 2> /dev/null
+    done > {output.sexavg} 2> /dev/null
     """
 
 rule mareymaps:
   input:
     data = "11_MareyMaps/marey.data.gz",
+    sexavg = "11_MareyMaps/marey.data.sexavg.gz",
     agp = expand("10_Anchoring/agp/chr.{lgs}.agp", lgs = lg_range)
   output: 
     indiv_plots = expand("11_MareyMaps/LG.{lgs}.mareymap.png", lgs = lg_range),
     summary = "11_MareyMaps/LepAnchor.mareymaps.pdf",
-    sequential = "11_MareyMaps/LepAnchor.sequentialmaps.pdf"
+    sequential = "11_MareyMaps/LepAnchor.sequentialmaps.pdf",
+    SAsummary = "11_MareyMaps/LepAnchor.sexavg.mareymaps.pdf",
+    SAsequential = "11_MareyMaps/LepAnchor.sexavg.sequentialmaps.pdf"
   message: "Creating Marey Maps"
   shell: 
     """
     Rscript software/LepAnchor/scripts/plot_marey.R {input.data} 10_Anchoring/agp
-    Rscript scripts/LA_summary.r {input.data}
+    Rscript scripts/LASummary.r {input.data}
+    Rscript scripts/LASummarySexAvg.r {input.sexavg}
     """

scripts/LASummary.r

@@ -0,0 +1,68 @@
+#! /usr/bin/env Rscript
+
+suppressMessages(library(tidyverse, warn.conflicts = FALSE, quietly = TRUE))
+
+args <- commandArgs(trailingOnly = TRUE)
+# args[1] = mareydata file
+
+allmaps <- suppressMessages(read_tsv(gzfile(args[1]), col_names = FALSE)) %>%
+  select(X3, X2, X5, X6) %>%
+  rename(lg = X3, Mb = X2, male = X5, female = X6) %>%
+  mutate(Mb = Mb/1000000) %>%
+  group_by(lg) %>%
+  mutate(marker = seq_along(Mb)) %>%
+  pivot_longer(c(male, female), names_to = "sex", values_to = "cM" )
+
+allmaps %>%
+  ggplot(aes(x = Mb, y = cM, color = sex)) +
+  geom_point(size = 0.6)  +
+  labs(
+      title = "Marker Positions in Anchored+Oriented Assembly",
+      subtitle = "The relative marker positions vs their position in the chromosome/LG",
+      x = "Physical Position (Mbp)",
+      y = "Genetic Position (cM)"
+      ) +
+  theme(
+        legend.position = "top",
+        legend.margin = margin(0, 0, 0, 0),
+        legend.box.margin = margin(-5, -5, -5, -5),
+        legend.title = element_blank(),
+        legend.text = element_text(size = 13),
+        axis.text.x = element_text(colour = "black", size = 11),
+        axis.text.y = element_text(colour = "black", size = 11),
+        strip.text = element_text(size = 13),
+        axis.title.x = element_text(size = 13, margin = margin(t = 5, r = 0, b = 0, l = 0)),
+        axis.title.y = element_text(size = 13, margin = margin(t = 0, r = 5, b = 0, l = 0))
+        ) +
+  facet_wrap(~lg, ncol = 4, scales = "free")
+
+outfile <- paste0(dirname(args[1]), "/LepAnchor.mareymaps.pdf")
+savedims <- length(unique(allmaps$lg)) * 2.5
+
+ggsave(outfile, width = 8.5, height = savedims/4, units = "in")
+
+allmaps %>%
+  ggplot(aes(x = marker, y = cM)) +
+  geom_point(size = 0.6, alpha = 0.5)  +
+  labs(
+      title = "Relative Marker Positions within Linkge Groups",
+      subtitle =  "The distance of sequential markers from each other in the linkage maps",
+      x = "Marker Number",
+      y = "Genetic Position (cM)"
+      ) +
+  theme(
+        legend.position = "top",
+        legend.margin = margin(0, 0, 0, 0),
+        legend.box.margin = margin(-5, -5, -5, -5),
+        legend.title = element_blank(),
+        legend.text = element_text(size = 13),
+        axis.text.x = element_text(colour = "black", size = 11),
+        axis.text.y = element_text(colour = "black", size = 11),
+        strip.text = element_text(size = 13),
+        axis.title.x = element_text(size = 13, margin = margin(t = 5, r = 0, b = 0, l = 0)),
+        axis.title.y = element_text(size = 13, margin = margin(t = 0, r = 5, b = 0, l = 0))
+        ) +
+  facet_wrap(lg ~ sex, ncol = 2, scales = "free_x")
+
+outfile2 <- paste0(dirname(args[1]), "/LepAnchor.sequentialmaps.pdf")
+ggsave(outfile2,  width = 8.5, height = savedims, units = "in", limitsize = FALSE)

scripts/LASummarySexAvg.r

@@ -0,0 +1,68 @@
+#! /usr/bin/env Rscript
+
+suppressMessages(library(tidyverse, warn.conflicts = FALSE, quietly = TRUE))
+
+args <- commandArgs(trailingOnly = TRUE)
+# args[1] = mareydata file
+# args[2] = marey midpoint data (optional)
+
+allmaps <- suppressMessages(read_tsv(gzfile(args[1]), col_names = FALSE)) %>%
+  select(X3, X2, X5) %>%
+  rename(lg = X3, Mb = X2, cM = X5) %>%
+  mutate(Mb = Mb/1000000) %>%
+  group_by(lg) %>%
+  mutate(marker = seq_along(Mb))
+
+allmaps %>%
+  ggplot(aes(x = Mb, y = cM)) +
+  geom_point(size = 0.6, color = "dodgerblue")  +
+  labs(
+      title = "Marker Positions in Anchored+Oriented Assembly",
+      subtitle = "The relative marker positions vs their position in the chromosome/LG",
+      x = "Physical Position (Mbp)",
+      y = "Genetic Position (cM)"
+      ) +
+  theme(
+        legend.position = "top",
+        legend.margin = margin(0, 0, 0, 0),
+        legend.box.margin = margin(-5, -5, -5, -5),
+        legend.title = element_blank(),
+        legend.text = element_text(size = 13),
+        axis.text.x = element_text(colour = "black", size = 11),
+        axis.text.y = element_text(colour = "black", size = 11),
+        strip.text = element_text(size = 13),
+        axis.title.x = element_text(size = 13, margin = margin(t = 5, r = 0, b = 0, l = 0)),
+        axis.title.y = element_text(size = 13, margin = margin(t = 0, r = 5, b = 0, l = 0))
+        ) +
+  facet_wrap(~lg, ncol = 4, scales = "free")
+
+outfile <- paste0(dirname(args[1]), "/LepAnchor.sexavg.mareymaps.pdf")
+savedims <- length(unique(allmaps$lg)) * 2.5
+
+ggsave(outfile, width = 8.5, height = savedims/4, units = "in")
+
+allmaps %>%
+  ggplot(aes(x = marker, y = cM)) +
+  geom_point(size = 0.6, alpha = 0.5, color = "dodgerblue")  +
+  labs(
+      title = "Relative Marker Positions within Linkge Groups",
+      subtitle =  "The distance of sequential markers from each other in the linkage maps",
+      x = "Marker Number",
+      y = "Genetic Position (cM)"
+      ) +
+  theme(
+        legend.position = "top",
+        legend.margin = margin(0, 0, 0, 0),
+        legend.box.margin = margin(-5, -5, -5, -5),
+        legend.title = element_blank(),
+        legend.text = element_text(size = 13),
+        axis.text.x = element_text(colour = "black", size = 11),
+        axis.text.y = element_text(colour = "black", size = 11),
+        strip.text = element_text(size = 13),
+        axis.title.x = element_text(size = 13, margin = margin(t = 5, r = 0, b = 0, l = 0)),
+        axis.title.y = element_text(size = 13, margin = margin(t = 0, r = 5, b = 0, l = 0))
+        ) +
+  facet_wrap(~lg, ncol = 4, scales = "free_x")
+
+outfile2 <- paste0(dirname(args[1]), "/LepAnchor.sexavg.sequentialmaps.pdf")
+ggsave(outfile2,  width = 8.5, height = savedims/5, units = "in", limitsize = FALSE)