Skip to content

GitLab

Switch branch/tag
Find file Normal viewHistoryPermalink
Compare_variants_AF_v3.R 5.41 KB
Newer Older
Matteo Barcella's avatar
typos    
Matteo Barcella committed
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122 
Compare_variants_AF_v3 <- function(dblist, sampleid_1, sampleid_2, 
                                my.width, my.height, my.res, 
                                tag = "demotag",
                                outfolder = "Outfold/"){
  
  library(ggplot2)
  library(plotly)
  library(ggrepel)
  library(reshape2)
  library(cowplot)
  library(openxlsx)
  
  dir.create(path = outfolder, recursive = T)
  
  db <- dblist
  
  Common_vars <- intersect(x = db[[sampleid_1]]$variantkey, y = db[[sampleid_2]]$variantkey)
  
  fulldf <- rbind.data.frame(x = db[[sampleid_1]], y = db[[sampleid_2]])
  
  fulldf_annot <- subset.data.frame(x = fulldf, select = c("variantkey","ID","HGVS_P","Effect","Impact", 
                                                           "dbNSFP_CADD_phred", "dbNSFP_PROVEAN_pred",
                                                           "dbNSFP_MutationTaster_pred", "dbNSFP_SIFT_pred",
                                                           "dbNSFP_Polyphen2_HVAR_pred"))
  fulldf_annot <- fulldf_annot %>% distinct(variantkey,.keep_all = TRUE)
  
  sampleid_1_sub <- subset.data.frame(x = db[[sampleid_1]], subset = variantkey %in% Common_vars, 
                                      select = c("variantkey","AF","DP","FILTER","Gene"))
  colnames(sampleid_1_sub) <- c("variantkey","AF_GFP_H","DP_GFP_H","FILTER_GFP_H","Gene")
  
  sampleid_2_sub <- subset.data.frame(x = db[[sampleid_2]], subset = variantkey %in% Common_vars, 
                                      select = c("variantkey","AF","DP","FILTER","Gene"))
  colnames(sampleid_2_sub) <- c("variantkey","AF_GFP_L","DP_GFP_L","FILTER_GFP_L","Gene")
  
  # merge datasets
  
  common_df <- merge.data.frame(x = sampleid_1_sub, y = sampleid_2_sub, by = c('variantkey','Gene'))
  common_df_melted <- melt(common_df, id.vars = c("variantkey","Gene", "DP_GFP_H", "DP_GFP_L","FILTER_GFP_H","FILTER_GFP_L"))
  
  # plotting common variants AF
  
  common_df_tmp <- common_df
  colnames(common_df_tmp) = c("variantkey","Gene","AF_GFP_H","DP_GFP_H","FILTER_GFP_H","AF_GFP_L","DP_GFP_L","FILTER_GFP_L")
  
  common_df_melted_tmp <- melt(common_df_tmp, id.vars = c("variantkey","Gene", "DP_GFP_H", "DP_GFP_L","FILTER_GFP_H","FILTER_GFP_L"))
  
  g_point <- ggplot(data = common_df_tmp, mapping = aes(x = AF_GFP_H, y = AF_GFP_L)) + geom_point() + xlab(label = "AF_GFP_H") + ylab(label = "AF_GFP_L")
  g_density <- ggplot(data = common_df_melted_tmp, mapping = aes(value, fill = variable)) + geom_density() + 
    facet_wrap(facets = ~ variable, nrow = 2) + theme(legend.position = "none", axis.title = element_blank())
  
  title <- ggdraw() +
    draw_label(
      paste0("\tAllele Fraction\n", sampleid_1, " vs ", sampleid_2),
      fontface = 'bold',
      x = 0,
      hjust = -0.5
    ) +
    theme(
      # add margin on the left of the drawing canvas,
      # so title is aligned with left edge of first plot
      plot.margin = margin(1, 1, 1, 1)
    )

  aa <- plot_grid(plotlist = list(g_density, g_point), ncol = 2)
  
  png(filename = paste0(outfolder,"AF_Comparison_",sampleid_1,"_",sampleid_2,"_common_",tag,".png"), 
      width = my.width, height = my.height, units = "in", res = my.res)
  print(plot_grid(
    title, 
    aa,
    ncol = 1,
    # rel_heights values control vertical title margins
    rel_heights = c(0.1, 1)
  ))
  dev.off()
  
  # Merging
  
   union_vars <- union(x = db[[sampleid_1]]$variantkey, y = db[[sampleid_2]]$variantkey)
   
   sampleid_1_full <- subset.data.frame(x = db[[sampleid_1]], select = c("variantkey","AF","DP","FILTER","Gene"))
   colnames(sampleid_1_full) <- c("variantkey","AF_GFP_H","DP_GFP_H","FILTER_GFP_H","Gene")
  # 
   sampleid_2_full <- subset.data.frame(x = db[[sampleid_2]], select = c("variantkey","AF","DP","FILTER","Gene"))
   colnames(sampleid_2_full) <- c("variantkey","AF_GFP_L","DP_GFP_L","FILTER_GFP_L","Gene")
  # 
  union_df <- merge.data.frame(x = sampleid_1_full, y=sampleid_2_full, by = c("variantkey","Gene"), all = T)
 
  union_df$AF_GFP_H[is.na(union_df$AF_GFP_H)] <- 0
  union_df$AF_GFP_L[is.na(union_df$AF_GFP_L)] <- 0

  union_df$DP_GFP_H[is.na(union_df$DP_GFP_H)] <- 0
  union_df$DP_GFP_L[is.na(union_df$DP_GFP_L)] <- 0
  
  colnames(union_df) <- c("variantkey","Gene","AF_GFP_H","DP_GFP_H","FILTER_GFP_H","AF_GFP_L","DP_GFP_L","FILTER_GFP_L")
  
  # Plotting Both

  png(filename = paste0(outfolder,"AF_Comparison_",sampleid_1,"_",sampleid_2,"_union_",tag,".png"), width = 12, height = 9, units = "in", res = 200)
  print(ggplot(data = union_df, mapping = aes(x = AF_GFP_H, y = AF_GFP_L, label = Gene)) + 
          geom_point(size = 0.5) + #geom_text(size=2) + 
          ggtitle(label = paste0("AF distribution ", tag), subtitle = paste0(sampleid_1, " vs ", sampleid_2)) +
          theme(plot.subtitle = element_text(hjust = 0.5))) 
  dev.off()

  # retrieve info from db for complete annotation 
  
  common_annot <- merge.data.frame(x = common_df, y = fulldf_annot, by = 'variantkey', all.x = T, sort = F)
  union_annot <- merge.data.frame(x = union_df, y = fulldf_annot, by = 'variantkey', all.x = T, sort = F)
    
  resultdata <- list(common = common_df, commonannot = common_annot, 
                     union = union_df, unionannot = union_annot,
                     comvars = fulldf)
  
  write.xlsx(x = resultdata, file = paste0(outfolder,"Comparison_",sampleid_1,"_",sampleid_2,"_", tag, ".xlsx"), asTable = T)
  
  saveRDS(object = resultdata, file = paste0(outfolder,"Comparison_",sampleid_1,"_",sampleid_2,"_", tag, ".rds"))
  
  return(resultdata)
}