merip_reanalysis_scripts/plot_fig2.R at master · mywanuo/merip_reanalysis_scripts · GitHub

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
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
library(ggplot2)
library(ComplexHeatmap)
library(circlize) #for colorRamp2
#also requires ggpubr

args <- commandArgs(TRUE)
peakcaller <- args[1]
thresh <- args[2] #used only for file name in this script
exp.summary <- args[3]
fig <- args[4]

cluster <- TRUE

format.annotations <- function(peakcaller,thresh,exp.summary,fig){

  peaks.overlap <- read.csv(paste0('fig2/peak_overlaps_summary_',peakcaller,'_',thresh,'.txt'),sep=' ')
  peaks.overlap$exp1 <- gsub("_"," ",peaks.overlap$exp1)
  peaks.overlap$exp2 <- gsub("_"," ",peaks.overlap$exp2)
  peaks.overlap$cellline1 <- gsub("_"," ",peaks.overlap$cellline1)
  peaks.overlap$cellline2 <- gsub("_"," ",peaks.overlap$cellline2)
  overlap.annotations <- read.csv(exp.summary,sep=' ')
  overlap.annotations$cell_line <- gsub("_"," ",overlap.annotations$cell_line)
  overlap.annotations$label <- gsub("_"," ",overlap.annotations$label)
  if (fig == 'xiao'){
    overlap.annotations <- overlap.annotations[which(overlap.annotations$fig == "xiao"),]
    overlap.annotations$explabel <- gsub("_"," ",paste0(overlap.annotations$label))
  }else{
    overlap.annotations <- overlap.annotations[which(overlap.annotations$fig == "2"),]
    overlap.annotations$explabel <- gsub("_"," ",paste0(overlap.annotations$label," (",overlap.annotations$cell_line,")"))
  }

  rownames(overlap.annotations) <- paste(overlap.annotations$label,overlap.annotations$cell_line)
  if (fig == 'xiao'){ overlap.annotations$label <- regmatches(overlap.annotations$control, regexpr("[[:digit:]]+", overlap.annotations$control))}
  #write.table(overlap.annotations)
  peaks.overlap$explabel1 <- overlap.annotations[paste(peaks.overlap$exp1,peaks.overlap$cellline1),"explabel"]
  peaks.overlap$explabel2 <- overlap.annotations[paste(peaks.overlap$exp2,peaks.overlap$cellline2),"explabel"]

  peaks.overlap$exp1 <- as.character(peaks.overlap$exp1)
  peaks.overlap$exp2 <- as.character(peaks.overlap$exp2)
  x <- list()
  x$peaks.overlap <- peaks.overlap
  x$overlap.annotations <- overlap.annotations
  return(x)
}

#plot Fig 2a
plot.fig2a <- function(peaks.overlap,overlap.annotations,finame){
  for (cellline in unique(peaks.overlap$cellline1)){
    cell.peak.ol <- peaks.overlap[which(peaks.overlap$cellline1 == cellline & peaks.overlap$cellline2 == cellline),]
    if (dim(cell.peak.ol)[1] > 0){
      perc.overlap.mat <- matrix(data = NA,nrow=length(unique(cell.peak.ol$exp1)),ncol=length(unique(cell.peak.ol$exp2)))
      perc.overlap.mat <- as.data.frame(perc.overlap.mat)
      rownames(perc.overlap.mat) <- unique(cell.peak.ol$exp1)
      colnames(perc.overlap.mat) <- unique(cell.peak.ol$exp1)
      for (exp1 in cell.peak.ol$exp1){
        perc.overlap.mat[exp1,exp1] <- 100
        for (exp2 in cell.peak.ol[which(cell.peak.ol$exp1 == exp1),'exp2']){
          peak.overlap <- cell.peak.ol[which(cell.peak.ol$exp1 == exp1 & cell.peak.ol$exp2 == exp2),'peak_overlap'][1]
          total.peaks <- cell.peak.ol[which(cell.peak.ol$exp1 == exp1 & cell.peak.ol$exp2 == exp2),'total_peaks'][1]
          perc.overlap.mat[exp1,exp2] <- peak.overlap*100/total.peaks
        }
      }

      diag(perc.overlap.mat)=NA
      perc.overlap.sum <- as.vector(as.matrix(perc.overlap.mat))
      summary <- paste(summary(perc.overlap.sum))
      writeLines(paste(cellline,'1st quartile = ',summary[2],', median = ',summary[3],', 3rd quartile = ',summary[5]))
      writeLines(paste(cellline,'min = ',summary[1],', max = ',summary[6]))

      diag(perc.overlap.mat) = 100
      #could add the total #peaks per comparison
      pdf(paste0("fig2/fig2a_",cellline,"_peak_heatmap_",finame,".pdf"),width=3.5,height=3)
      ch = Heatmap(perc.overlap.mat,col=colorRamp2(c(0,100),c("white","#89043d")), cluster_rows = FALSE, cluster_columns = FALSE,
              column_names_side = "top",row_names_side = "left",name="% peaks\noverlapped",rect_gp = gpar(col = "white", lwd = 2),
              column_title = "Experiment 2",row_title = "Experiment 1",column_title_gp = gpar(fontsize=12),row_title_gp = gpar(fontsize=12))
      draw(ch, column_title=cellline,column_title_gp = gpar(fontsize=14))
      dev.off()
    }
  }
}

#plot Fig 2b
plot.fig2b <- function(peaks.overlap,overlap.annotations,finame,feature1="cell line",feature2="study",fig='2'){
  #write.table(peaks.overlap)
  #break()
  for (species in c("hg38")){
    cell.lines <- overlap.annotations[which(overlap.annotations$species == species & overlap.annotations$fig == fig),"cell_line"]
    cell.peak.ol <- peaks.overlap[which(peaks.overlap$cellline1 %in% cell.lines),]
    perc.overlap.mat <- matrix(data = NA,nrow=length(unique(cell.peak.ol$explabel1)),ncol=length(unique(cell.peak.ol$explabel2)))
    perc.overlap.mat <- as.data.frame(perc.overlap.mat)
    rownames(perc.overlap.mat) <- unique(cell.peak.ol$explabel1)
    colnames(perc.overlap.mat) <- unique(cell.peak.ol$explabel2)
    for (explabel1 in cell.peak.ol$explabel1){
      perc.overlap.mat[explabel1,explabel1] <- 100
      for (explabel2 in cell.peak.ol[which(cell.peak.ol$explabel1 == explabel1),'explabel2']){
        peak.overlap <- cell.peak.ol[which(cell.peak.ol$explabel1 == explabel1 & cell.peak.ol$explabel2 == explabel2),'peak_overlap'][1]
        total.peaks <- cell.peak.ol[which(cell.peak.ol$explabel1 == explabel1 & cell.peak.ol$explabel2 == explabel2),'total_peaks'][1]
        perc.overlap.mat[explabel1,explabel2] <- peak.overlap*100/total.peaks
      }
    }
  }

  row.order <- sort(overlap.annotations[which(overlap.annotations$explabel %in% rownames(perc.overlap.mat)),]$explabel)
  #write.table(row.order)
  perc.overlap.mat <- perc.overlap.mat[row.order,row.order]
  o.a <- overlap.annotations[which(overlap.annotations$explabel %in% row.order),]
  #write.table(o.a$explabel)
  rownames(o.a) <- o.a$explabel
  o.a <- o.a[row.order,]
  cl <- sort(unique(gsub("_"," ",o.a$cell_line)))
  cell.colours <- list(feature1= setNames(ggpubr::get_palette(c("#3f1a1c","#d78521","#e6af2e","#ede5a6"),length(cl)),cl))
  st <- unique(o.a$label)
  if (feature2 == "gestation time (w)"){f2.pal <- ggpubr::get_palette("Blues",length(st))
  }else{ f2.pal <- ggpubr::get_palette("lancet",length(st))}
  study.colours <- list(feature2= setNames(f2.pal,st))

  celllines = rowAnnotation(df = data.frame(feature1 = sub("_"," ",o.a$cell_line)),
                                col = cell.colours,annotation_legend_param = list(feature1 = list(title = feature1)),
                            show_annotation_name = FALSE)
  studies = rowAnnotation(df = data.frame(feature2 = o.a$label),col = study.colours,
                          annotation_legend_param = list(feature2 = list(title = feature2)),show_annotation_name = FALSE)

  pdf(paste0("fig2/fig2b_clustered_peak_heatmap_",finame,".pdf"),width=8,height=6)
  ch = Heatmap(perc.overlap.mat,col=colorRamp2(c(0,100),c("white","#89043d")), cluster_rows = TRUE, cluster_columns = TRUE,
               column_names_side = "top",row_names_side = "left",name="% peaks\noverlapped",rect_gp = gpar(col = "white", lwd = 2),
               column_title = "Experiment 2",row_title = "Experiment 1",column_title_gp = gpar(fontsize=12),row_title_gp = gpar(fontsize=12))
  draw(ch + celllines + studies, heatmap_legend_side = "right")
  dev.off()

  writeLines(paste0("fig2/fig2b_clustered_peak_heatmap_",finame,".pdf"))
  perc.overlap.sum <- as.vector(as.matrix(perc.overlap.mat))
  perc.overlap.sum <- perc.overlap.sum[which(perc.overlap.sum != 100)]
  summary(perc.overlap.sum)
}

x <- format.annotations(peakcaller,thresh,exp.summary,fig)
if (fig == "xiao"){
  plot.fig2b(x$peaks.overlap,x$overlap.annotations,paste0(fig,thresh),"tissue","gestation time (w)",fig)
}else{
  plot.fig2a(x$peaks.overlap,x$overlap.annotations,peakcaller)
  plot.fig2b(x$peaks.overlap,x$overlap.annotations,paste0(peakcaller,thresh))
}