update mapping plotting scripts

jonassibbesen · Jun 28, 2022 · 39f7d9b · 39f7d9b
1 parent e14771e
commit 39f7d9b
Show file tree

Hide file tree

Showing 4 changed files with 24 additions and 12 deletions.
diff --git a/R/plot_coverage_correlation.R b/R/plot_coverage_correlation.R
@@ -49,7 +49,8 @@ coverage_data$Method <- recode_factor(coverage_data$Method,
                                                        "star" = "STAR",
                                                        "map_fast" = "vg map",
                                                        "mpmap" = "vg mpmap", 
-                                                       "star_alleleseq" = "Diploid reference (STAR)")
+                                                       "star_alleleseq" = "Diploid reference (STAR)",
+                                                       "star_alleleseq_levio" = "Diploid reference (STAR, LevioSAM)")
 
 coverage_data$Reference = recode_factor(coverage_data$Reference, 
                                         "1kg_nonCEU_af001_gencode100" = "Spliced pangenome graph",
@@ -107,15 +108,15 @@ for (reads in unique(coverage_data_pb_mq_corr_main$Reads)) {
 
 coverage_data_pb_mq_corr_personal <- coverage_data_pb_mq_corr %>%
   filter(Reads == "ENCSR000AED_rep1") %>%
-  filter(Method == "STAR" | ((Method == "vg mpmap" | Method == "Diploid reference (STAR)") & Reference == "Spliced personal graph/reference"))
+  filter(Method == "STAR" | ((Method == "vg mpmap" | Method == "Diploid reference (STAR)" | Method == "Diploid reference (STAR, LevioSAM)") & Reference == "Spliced personal graph/reference"))
 
 for (reads in unique(coverage_data_pb_mq_corr_personal$Reads)) {
 
   coverage_data_pb_mq_corr_personal_reads <- coverage_data_pb_mq_corr_personal %>%
     filter(Reads == reads) %>%
     rename(MapQ = Threshold)
 
-  plotIsoSeqCorrelationBenchmark(coverage_data_pb_mq_corr_personal_reads, wes_cols[c(2,4,5)], paste("plots/real_corr/real_r2_cov_corr_personal_", reads, sep = ""))
+  plotIsoSeqCorrelationBenchmark(coverage_data_pb_mq_corr_personal_reads, wes_cols[c(2,4,5,6)], paste("plots/real_corr/real_r2_cov_corr_personal_", reads, sep = ""))
 }
 
 ########

diff --git a/R/plot_mapping_bias.R b/R/plot_mapping_bias.R
@@ -92,6 +92,7 @@ coverage_data_mq_bias$Method = recode_factor(coverage_data_mq_bias$Method,
                                              "mpmap" = "vg mpmap",
                                              "mpmap_multi10" = "vg mpmap",
                                              "star_alleleseq" = "Diploid reference (STAR)",
+                                             "star_alleleseq_levio" = "Diploid reference (STAR, LevioSAM)",
                                              "star_wasp" = "WASP (STAR)")
 
 coverage_data_mq_bias[coverage_data_mq_bias$Method == "WASP (STAR)",]$Reference <- "1kg_NA12878_gencode100"
@@ -109,6 +110,11 @@ for (reads in unique(coverage_data_mq_bias_debug$Reads)) {
     filter(Reads == reads)
 
   plotMappingBiasBenchmark(coverage_data_mq_bias_debug_reads, wes_cols, paste("plots/sim_bias/debug/vg_sim_r2_mapping_bias_debug_", reads, sep = ""), 20)
+
+  coverage_data_mq_bias_debug_reads$FacetCol <- coverage_data_mq_bias_debug_reads$var
+  coverage_data_mq_bias_debug_reads$FacetRow <- paste(coverage_data_mq_bias_debug_reads$Simulation, ",\n primary alignments", sep = "")
+
+  plotMappingBiasBinomBenchmark(coverage_data_mq_bias_debug_reads, wes_cols, paste("plots/sim_bias/debug/vg_sim_r2_mapping_bias_binom_a001_debug_", reads, sep = ""), 0.01, 20)
 }
 
 ########
@@ -117,6 +123,7 @@ coverage_data_mq_bias_main <- coverage_data_mq_bias %>%
   filter(Reads == "sim_vg_r2_ENCSR000AED_rep1_uni") %>%
   filter(Method != "WASP (STAR)") %>%
   filter(Method != "Diploid reference (STAR)") %>%
+  filter(Method != "Diploid reference (STAR, LevioSAM)") %>%
   filter(Method != "vg map (def)") %>%
   filter(Reference != "1kg_NA12878_gencode100") %>%
   filter(Reference != "1kg_NA12878_exons_gencode100")
@@ -139,7 +146,9 @@ for (reads in unique(coverage_data_mq_bias_main$Reads)) {
 
 coverage_data_mq_bias_binom <- coverage_data_mq_bias %>%
   filter(Reads == "sim_vg_r2_ENCSR000AED_rep1_uni") %>%
-  filter(Method != "vg map (def)") 
+  filter(Method != "vg map (def)") %>%
+  filter(Method != "Diploid reference (STAR)") %>%
+  filter(Method != "Diploid reference (STAR, LevioSAM)")
 
 coverage_data_mq_bias_binom$Reference = recode_factor(coverage_data_mq_bias_binom$Reference,
                                                  "1kg_nonCEU_af001_gencode100" = "Spliced pangenome graph",

diff --git a/R/plot_mapping_stats.R b/R/plot_mapping_stats.R
@@ -39,7 +39,8 @@ mapping_data$Method <- recode_factor(mapping_data$Method,
                                      "star" = "STAR",
                                      "map_fast" = "vg map", 
                                      "mpmap" = "vg mpmap",
-                                     "star_alleleseq" = "Diploid reference (STAR)")
+                                     "star_alleleseq" = "Diploid reference (STAR)",
+                                     "star_alleleseq_levio" = "Diploid reference (STAR, LevioSAM)")
 
 
 mapping_data$Reference = recode_factor(mapping_data$Reference, 
@@ -93,8 +94,8 @@ for (reads in unique(mapping_data_stats_main$Reads)) {
 
 mapping_data_stats_personal <- mapping_data_stats %>%
   filter(Reads == "ENCSR000AED_rep1") %>%
-  filter(Method == "STAR" | ((Method == "vg mpmap" | Method == "Diploid reference (STAR)") & Reference == "Spliced personal graph/reference"))
-
+  filter(Method == "STAR" | ((Method == "vg mpmap" | Method == "Diploid reference (STAR)" | Method == "Diploid reference (STAR, LevioSAM)") & Reference == "Spliced personal graph/reference"))
+  
 for (reads in unique(mapping_data_stats_personal$Reads)) {
 
   mapping_data_stats_personal_reads <- mapping_data_stats_personal %>%
@@ -105,7 +106,7 @@ for (reads in unique(mapping_data_stats_personal$Reads)) {
     add_row(Reads = reads, Method = "vg mpmap", Reference = "Spliced reference", count = 0, MapQ0 = 0, MapQ30 = 0, Filter = "MapQ0_frac", Frac = 0) %>%
     add_row(Reads = reads, Method = "vg mpmap", Reference = "Spliced reference", count = 0, MapQ0 = 0, MapQ30 = 0, Filter = "MapQ30_frac", Frac = 0)
 
-  mapping_data_stats_personal_reads$Method <- factor(mapping_data_stats_personal_reads$Method, levels = c("STAR", "vg mpmap", "Diploid reference (STAR)"))
+  mapping_data_stats_personal_reads$Method <- factor(mapping_data_stats_personal_reads$Method, levels = c("STAR", "vg mpmap", "Diploid reference (STAR)", "Diploid reference (STAR, LevioSAM)"))
 
   mapping_data_stats_personal_reads$Reference = recode_factor(mapping_data_stats_personal_reads$Reference,
                                                           "Spliced reference" = "Spliced\nreference",
@@ -118,7 +119,7 @@ for (reads in unique(mapping_data_stats_personal$Reads)) {
   mapping_data_stats_personal_reads$FacetCol <- "Real reads,\nprimary alignments"
   mapping_data_stats_personal_reads$FacetRow <- ""
 
-  plotMappingStatsBenchmarkWide(mapping_data_stats_personal_reads, wes_cols[c(2,4,5)], paste("plots/real_stats/real_r2_stats_bar_personal_", reads, sep = ""))
+  plotMappingStatsBenchmarkWide(mapping_data_stats_personal_reads, wes_cols[c(2,4,5,6)], paste("plots/real_stats/real_r2_stats_bar_personal_", reads, sep = ""))
 }
 
 ########
diff --git a/R/plot_vg_sim_overlap.R b/R/plot_vg_sim_overlap.R
@@ -83,7 +83,7 @@ overlap_data %>% group_by(Reads, Reference, Method) %>%
   print(n = 100)
 
 overlap_data_prim1 <- overlap_data %>%
-  filter(Method == "hisat2_multi10" | Method == "star_multi10" | Method == "star_wasp" | Method == "star_alleleseq" | Method == "mpmap_multi10") %>%
+  filter(Method == "hisat2_multi10" | Method == "star_multi10" | Method == "star_wasp" | Method == "star_alleleseq" | Method == "star_alleleseq_levio" | Method == "mpmap_multi10") %>%
   mutate(MapQ = PrimaryMapq) %>%
   mutate(Overlap = PrimaryOverlap) %>%
   mutate(Eval = "Primary")
@@ -111,6 +111,7 @@ overlap_data$Method <- recode_factor(overlap_data$Method,
                                      "mpmap" = "vg mpmap",
                                      "mpmap_multi10" = "vg mpmap",
                                      "star_alleleseq" = "Diploid reference (STAR)",
+                                     "star_alleleseq_levio" = "Diploid reference (STAR, LevioSAM)",
                                      "star_wasp" = "WASP (STAR)")
 
 overlap_data[overlap_data$Method == "WASP (STAR)",]$Reference <- "1kg_NA12878_gencode100"
@@ -352,7 +353,7 @@ for (reads in unique(overlap_data_main_multi_indel$Reads)) {
 
 overlap_data_personal_prim <- overlap_data %>%
   filter(Reads == "sim_vg_r2_ENCSR000AED_rep1_uni") %>%
-  filter((Method == "STAR" & Reference == "Spliced reference") | ((Method == "vg mpmap" | Method == "Diploid reference (STAR)") & Reference == "Spliced personal graph/reference")) %>%
+  filter((Method == "STAR" & Reference == "Spliced reference") | ((Method == "vg mpmap" | Method == "Diploid reference (STAR)" | Method == "Diploid reference (STAR, LevioSAM)") & Reference == "Spliced personal graph/reference")) %>%
   filter(Eval == "Primary") %>%
   mutate(FacetCol = paste(FacetCol, ", primary alignments", sep = ""))
 
@@ -361,7 +362,7 @@ for (reads in unique(overlap_data_personal_prim$Reads)) {
   overlap_data_personal_prim_reads <- overlap_data_personal_prim %>%
     filter(Reads == reads)
 
-  plotRocBenchmarkMapQ(overlap_data_personal_prim_reads, wes_cols[c(2,4,5)], "Reference", paste("plots/sim_overlap/vg_sim_r2_overlap_personal_ovl", overlap_threshold, "_", reads, "_primary", sep = ""))
+  plotRocBenchmarkMapQ(overlap_data_personal_prim_reads, wes_cols[c(2,4,5,6)], "Reference", paste("plots/sim_overlap/vg_sim_r2_overlap_personal_ovl", overlap_threshold, "_", reads, "_primary", sep = ""))
 }
 
 ########