updating AF comparison script

WES/AF_comparison.R

@@ -9,6 +9,7 @@ inputdata <- list()
 for (mysample in c("12_27_1","12_27_2", "10_395_3", "10_395_4")) {
   
   inputdata[["rawcall"]][[mysample]] <- filter_variants_v3(infile = paste0(mysample, "_dbsnp_cc_cnc_dbnsfp_refseq_light_raw_filtered.fields.txt"), sampleid = mysample, dp = 50)
+  inputdata[["passonly"]][[mysample]] <- filter_variants_v3(infile = paste0(mysample, "_dbsnp_cc_cnc_dbnsfp_refseq_light.fields.txt"), sampleid = mysample, dp = 50)
 }
 
 source(file = "Compare_variants_AF_v3.R")
@@ -69,6 +70,48 @@ plot_grid(AF_density_raw[["12-27_miRNA_L"]], AF_density_raw[["12-27_miRNA_H"]],
           nrow = 2)
 dev.off()
 
+## PLOTTING PASS-ONLY
+
+tmpdf_1_pass <- subset.data.frame(x = Comparison_12_27$passonly$unionannot, select = c("AF_GFP_H"), subset = AF_GFP_H > 0)
+tmpdf_2_pass <- subset.data.frame(x = Comparison_12_27$passonly$unionannot, select = c("AF_GFP_L"), subset = AF_GFP_L > 0)
+
+tmpdf_3_pass <- subset.data.frame(x = Comparison_10_395$passonly$unionannot, select = c("AF_GFP_H"), subset = AF_GFP_H > 0)
+tmpdf_4_pass <- subset.data.frame(x = Comparison_10_395$passonly$unionannot, select = c("AF_GFP_L"), subset = AF_GFP_L > 0)
+
+AF_density_pass <- list()
+
+AF_density_pass[["12-27_miRNA_L"]] <- ggplot(data = tmpdf_1_pass,
+                                        mapping = aes(AF_GFP_H)) + ylim(0,300) +
+  geom_histogram(bins = 50, fill = alpha("#003366",0.8)) + xlab("miRNA_low") +
+  theme(axis.text = element_text(size = 14), legend.title = element_blank()) +
+  ggtitle("12-27_miRNA_L")
+
+AF_density_pass[["12-27_miRNA_H"]] <- ggplot(data = tmpdf_2_pass,
+                                        mapping = aes(AF_GFP_L)) + ylim(0,300) +
+  geom_histogram(bins = 50, fill = alpha("#fc2803",0.8))+ xlab("miRNA_high") +
+  theme(axis.text = element_text(size = 14), legend.title = element_blank()) +
+  ggtitle("12-27_miRNA_H")
+
+AF_density_pass[["10-395_miRNA_L"]] <- ggplot(data = tmpdf_3_pass,
+                                         mapping = aes(AF_GFP_H)) + ylim(0,300) +
+  geom_histogram(bins = 50, fill = alpha("#003366",0.8)) + xlab("miRNA_low") +
+  theme(axis.text = element_text(size = 14), legend.title = element_blank()) +
+  ggtitle("10-395_miRNA_L")
+
+AF_density_pass[["10-395_miRNA_H"]] <- ggplot(data = tmpdf_4_pass,
+                                         mapping = aes(AF_GFP_L)) + ylim(0,300) +
+  geom_histogram(bins = 50, fill = alpha("#fc2803",0.8)) + xlab("miRNA_high") +
+  theme(axis.text = element_text(size = 14), legend.title = element_blank()) +
+  ggtitle("10-395_miRNA_H")
+
+png(filename = "AF_density_pass_variants.png", width = 8, height = 6, units = "in", res = 300)
+plot_grid(AF_density_pass[["12-27_miRNA_L"]], AF_density_pass[["12-27_miRNA_H"]],
+          AF_density_pass[["10-395_miRNA_L"]] , AF_density_pass[["10-395_miRNA_H"]],
+          nrow = 2)
+dev.off()
+
+# Histograms with PASS or clustered_events
+
 AF_density_mildfilt <- list()
 Mild_vars_12_27 <- Comparison_12_27$rawcall$unionannot
 
@@ -81,6 +124,7 @@ tmpdf_1_mild <- subset.data.frame(x = Mild_vars_12_27, select = c("AF_GFP_H"),
 tmpdf_2_mild <- subset.data.frame(x = Mild_vars_12_27, select = c("AF_GFP_L"), 
                                   subset = AF_GFP_L > 0 & (FILTER_GFP_H_simple == "PASS" | FILTER_GFP_L_simple == "PASS") & !HGVS_P == "")
 
+
 Mild_vars_10_395 <- Comparison_10_395$rawcall$unionannot
 Mild_vars_10_395$FILTER_GFP_H_simple <- ifelse(test = Mild_vars_10_395$FILTER_GFP_H == "PASS", "PASS","NOPASS")
 Mild_vars_10_395$FILTER_GFP_L_simple <- ifelse(test = Mild_vars_10_395$FILTER_GFP_L == "PASS", "PASS","NOPASS")
@@ -221,6 +265,8 @@ df_2 <- subset.data.frame(x = raw_vars_12_27_2, subset = variantkey %in% mykeys,
 colnames(df_2) <- c("variantkey","AF_GFP_L","DP_GFP_L","Gene", "Effect","FILTER")
 
 
+# modificare
+
 union_df <- merge.data.frame(x = df_1, y=df_2, by = c("variantkey","Gene","Effect"), all = T)
 
 
@@ -255,4 +301,3 @@ print(ggplot(data = union_df_exonic, mapping = aes(x = AF_GFP_H, y = AF_GFP_L, l
         geom_text_repel(size=3, show.legend = FALSE, segment.alpha = 0.4) +
         theme(plot.title = element_text(hjust = 0.5), legend.title = element_blank(), legend.position = "none"))
 dev.off()
-