一個函數實現基因內具有多種突變類型的熱圖的繪制
阿新 • • 發佈：2018-05-07
TP 必須 長度 部分 point 顯示效果 heatmap points bar
??我們平常多見的基因突變熱圖是一個基因一個格子，一種突變類型，但實際上在同一個病人中，同一個基因往往具有多種突變類型，因此傳統的熱圖繪制工具並不能滿足我們繪圖的需要。應研究需要，本人自己寫了一個熱圖繪制函數，內部調用image 進行熱圖的繪制， barplot進行直方圖繪制， 用data.table進行數據處理。對於一個基因內多種突變類型如何表現出來的問題， 這個函數先采用image將初步的熱圖繪制出來，再使用points，以方塊形式將第二種突變，第三種突變依次添加， 在添加的同時方塊位置稍為移動並且伴隨著大小的略微縮小，以實現更好的顯示效果，最多能在一個熱圖格子上表示四種突變。
??函數如下， 需要安裝並加載data.table 1.10.4， 加載RColorBrewer
my_heatmap <- function(vr, pal = c("#F2F2F2",colorRampPalette(c("blue", "white", "red"))(5)[c(1,2)],"#F2F2F2",colorRampPalette(c("blue", "white", "red"))(5)[c(4,5)],brewer.pal(n = 8, name ="Accent")[c(1,4,6,8,2,3,5,7)],"#E31A1C","#6A3D9A"),type = c("DEL","LOSS","NEUTRAL","GAIN","AMPL","nonsynonymous SNV","synonymous SNV","intronic","stopgain","nonframeshift deletion","splicing", "frameshift deletion","UTR3","frameshift insertion","UTR5"), order_gene = T, order_patient = T, hist_plot = T, legend_dist = 0.4, col_text_cex = 1, sub_gene= NULL,heatmap_mar = c(5,17,1,2), heatmap_oma=c(0.2,0.2,0.2,0.2),heatmap_mex=0.5, legend_mar = c(1,0,4,1),xlab_adj=1, order_omit=c("NoMut","NEUTRAL")) { if((length(pal) - length(type)) !=1 ){stop("Pal must be one longer than type, because first one pal is col for no mutation")} if(!is.null(sub_gene)){ pal_dt <- data.table(pal, type=c("NoMut",type)) vr <- vr[Gene %in% sub_gene,] type <- pal_dt[type %in% unique(vr$Type),type] pal <- c(pal[1],pal_dt[type, on="type"][,pal]) }else{ pal_dt <- data.table(pal, type=c("NoMut",type)) type <- pal_dt[type %in% unique(vr$Type),type] pal <- c(pal[1],pal_dt[type, on="type"][,pal]) } dt <- unique(vr[,.(Gene,Type,Patient)]) if(is.null(type)){type <- data.table(table(vr$Type))[order(-N),V1]} dt$Type <- factor(dt$Type, levels = type) gene <- dt[!Type %in% order_omit,.(N=length(unique(Patient))),by=Gene][order(N),Gene] dt$Gene <- factor(dt$Gene, levels = gene) patient <- data.table(table(vr[!Type %in% order_omit,]$Patient))[order(-N),V1] dt$Patient <- factor(dt$Patient, levels = c(patient, setdiff(unique(dt$Patient),patient))) if(order_gene & !order_patient){setkey(dt, "Gene")} if(!order_gene & order_patient){setkey(dt, "Patient")} if(order_gene & order_patient){setkey(dt, "Gene","Patient","Type")} n <- length(unique(dt$Type)) dt$Gene_Patients <- paste(dt$Gene, dt$Patient) dt_inf <- dt[,.N,by=.(Gene, Patient)] max_mut_num <- max(dt_inf$N) dt[,Mut_num:=seq_len(.N),by=.(Patient,Gene)] #main plot, heatmap using first mutation type dt1 <- copy(dt) dt1[Mut_num !=1, Type:=NA] dc <- data.frame(dcast(dt1, Patient ~ Gene, value.var = "Type", fun.aggregate = function(x)(x[!is.na(x)][1]))) rownames(dc)<- dc[,1] data_matrix<-data.matrix(dc[,-1]) data_matrix[is.na(data_matrix)] <- 0 pal=pal breaks<-seq(-1,10,1) if(!hist_plot){ layout(matrix(data=c(1,2), nrow=1, ncol=2), widths=c(8,2), heights=c(1,1)) par(mar=heatmap_mar, oma=heatmap_oma, mex=heatmap_mex) }else if(hist_plot){ layout(matrix(c(2,4,1,3),2,2,byrow=TRUE), widths=c(3,1), heights=c(1,3), TRUE) par(mar=heatmap_mar) } image(x=1:nrow(data_matrix),y=1:ncol(data_matrix), z=data_matrix,xlab="",ylab="",breaks=breaks, col=pal[1:11],axes=FALSE) #sub plot, points using other mutation type add_plot <- function(dt, i){ dt1 <- copy(dt) dt1[Mut_num != i, Type:=NA] dc <- data.frame(dcast(dt1, Patient ~ Gene, value.var = "Type", fun.aggregate = function(x){ifelse(length(x) >1,x[!is.na(x)][1],factor(NA))})) rownames(dc)<- dc[,1] data_matrix <- data.matrix(dc[,-1]) xy <- which(data_matrix !=0, arr.ind = T) #apply(xy, 1, function(x)points(x[1], x[2],pch=15, cex=2.5 -0.5*i, col=pal[data_matrix[x[1],x[2]]+1])) apply(xy, 1, function(x)points(x[1]-0.6+i*0.25, x[2],pch=15, cex=1.2 - i*0.08, col=pal[data_matrix[x[1],x[2]]+1])) } ploti <- data.frame(i=2:max_mut_num) apply(ploti, 1, function(i){print(add_plot(dt, i))}) text(x=1:nrow(data_matrix)+0.1, y=par("usr")[1] - xlab_adj, srt = 90, adj = 0.5, labels = rownames(data_matrix), xpd = TRUE, cex=col_text_cex) axis(2,at=1:ncol(data_matrix),labels=colnames(data_matrix), col="white",las=1, cex.lab=0.1) abline(h=c(1:ncol(data_matrix))+0.5,v=c(1:nrow(data_matrix))+0.5, col="white",lwd=2,xpd=F) #title("Correlation between genes",line=8,adj=0) if(hist_plot){ #bar plot par(mar=c(0,2+0.5,3,heatmap_mar[4]-0.9)) patient_dt <- dt[,.N,by=.(Patient,Type)] mt <- data.frame(dcast(patient_dt, Type ~ Patient, value.var = "N")) data_matrix <- data.matrix(mt[,-1]) rownames(data_matrix) <- mt[,1] tryCatch(data_matrix <- data_matrix[setdiff(type, order_omit), patient], error = function(e){print("type argument or your patient name format(include "-" and so on )")}) data_matrix[is.na(data_matrix)] <- 0 omit_idx <- NULL for(i in order_omit){omit_idx <- c(omit_idx,1+which(type == i))} barplot(data_matrix, col=pal[-c(1,omit_idx)],space=0,border = "white",axes=T,xlab="",ann=F, xaxt="n") par(mar=c( heatmap_mar[1]-2 , 0.8, heatmap_mar[3]+2.2, 3),las=1) gene_dt <- dt[,.N,by=.(Gene,Type)] mt <- data.frame(dcast(gene_dt, Type ~ Gene, value.var = "N")) data_matrix <- data.matrix(mt[,-1]) rownames(data_matrix) <- mt[,1] gene <- gsub("ATM,", "ATM.", gene) tryCatch(data_matrix <- data_matrix[setdiff(type, order_omit), gene], error = function(e){print("type argument or check your gene name format(please not include "-" and so on)")}) data_matrix[is.na(data_matrix)] <- 0 barplot(data_matrix, col=pal[-c(1,omit_idx)],space=0,border = "white",axes=T,xlab="", ann=F, horiz = T, yaxt="n") } #add legend par(mar=legend_mar) plot(3, 8, axes=F, ann=F, type="n") ploti <- data.frame(i=1:length(type)) if(!hist_plot){ tmp <- apply(ploti, 1, function(i){print(points(2, 10+(length(type)-i)*legend_dist, pch=15, cex=2, col=pal[i+1]))}) tmp <- apply(ploti, 1, function(i){print(text(3, 10+(length(type)-i)*legend_dist, labels = type[i],pch=15, cex=1, col="black"))}) } if(hist_plot){ tmp <- apply(ploti, 1, function(i){print(points(2, 5+(length(type)-i)*legend_dist, pch=15, cex=0.9, col=pal[i+1]))}) tmp <- apply(ploti, 1, function(i){print(text(2.8, 5+(length(type)-i)*legend_dist, labels = type[i],pch=15, cex=0.9, col="black"))}) } }
描述：
??繪制一個基因可以同時顯示多種突變類型的熱圖，輸入三列的data table數據框， 列名分別是Gene，Type， 和 Patient，輸出熱圖， 還可以在熱圖上方和右方添加突變的直方圖。
用法：
my_heatmap(vr, pal = c("#F2F2F2",colorRampPalette(c("blue", "white", "red"))(5)[c(1,2)],"#F2F2F2",colorRampPalette(c("blue", "white", "red"))(5)[c(4,5)],brewer.pal(n = 8, name ="Accent")[c(1,4,6,8,2,3,5,7)],"#E31A1C","#6A3D9A"),type = c("DEL","LOSS","NEUTRAL","GAIN","AMPL","nonsynonymous SNV","synonymous SNV","intronic","stopgain","nonframeshift deletion","splicing", "frameshift deletion","UTR3","frameshift insertion","UTR5"),order_gene = T, order_patient = T, hist_plot = T, legend_dist = 0.4, col_text_cex = 1,xlab_adj=1, sub_gene= NULL,heatmap_mar = c(5,17,1,2), heatmap_oma=c(0.2,0.2,0.2,0.2),heatmap_mex=0.5, legend_mar = c(1,0,4,1), order_omit=c("NEUTRAL"))
參數:
vr： 含有變異數據的數據框，共三列，列名分別是Gene， Type， Patient；
pal： 色板，向量，需要根據數據框中突變類型的數量進行自定義，需要比突變類型多一種顏色作為背景色，背景色放在第一位；
type：色板相對應的突變類型，向量，type必須等於或者多於數據中所出現的所有類型；默認使用拷貝數的四種突變類型加上拷貝數中性再加上annovar中所有突變類型；自定義設置時長度要比色板少1；
order_gene： 默認T， 對基因按照突變的病人數目進行排序；
order_patient：默認T，對病人按照突變的基因數目進行排序；
hist_plot：默認T，在上方和右方加上對應的直方圖；
legend_dist：默認0.4，調整圖例之間相互的距離，一般需要自行調整；
col_text_cex：調整病人名稱的大小，默認1；
xlab_adj：調整病人名稱與熱圖之間的距離；
sub_gene：只選擇部分基因進行畫圖，需要給基因名的向量，並且基因需要在數據中存在，默認NULL
heatmap_mar：mar參數，調整熱圖前後左右的邊緣長度，默認c(5,17,1,2)
heatmap_oma：oma參數，調整熱圖前後左右的外邊緣長度，默認c(0.2,0.2,0.2,0.2)
mex：調整熱圖的mex參數，用於描繪繪圖邊緣的坐標，默認0.5
legend_mar：legend的mar參數，調整圖例的位置，默認c(1,0,4,1)
order_omit：排序時忽略的變異類型，這些突變類型在直方圖中也會被過濾，默認c("NEUTRAL")，如果不存在"NEUTRAL"這種突變類型，也可以保持默認參數。
細節：
- 運行前需要加載data.table1.10.4， RColorBrewer；
- 如果要繪制帶有直方圖的熱圖，因為圖片尺寸過大，因此要使用pdf函數並要給足夠大的寬度和長度；
- 默認使用的是annovar註釋的突變類型；
- 因為繪制時影響熱圖和直方圖對齊的因素太多，很難通過調節相應的mar，mex，oma參數達到較好的效果，因此推薦快速畫出個大概後，再使用inkscape或adobe進行排版對齊
- 如果熱圖中突變類型的點過小，可以減小pdf文件的寬度和長度。
使用例子：
#without hist plot
pdf("~/project/PE/fromws02/PE/cnv_plot/heatmap_cnv_mut.pdf", height=12, width = 12)
my_heatmap(vr, heatmap_mar = c(17,17,1,2),hist_plot = F, legend_dist=0.1, xlab_adj = 1.2, order_patient = T, order_gene = T)
dev.off()
#with hist plot
pdf("~/project/PE/fromws02/PE/cnv_plot/heatmap_hist_cnv_mut.pdf", height=12, width = 12)
my_heatmap(vr, heatmap_mar = c(17,7,1,2),hist_plot = T, legend_dist=0.3, xlab_adj = 1.2, order_patient = T, order_gene = T)
dev.off()
#only a few gene
pdf("~/project/PE/fromws02/PE/cnv_plot/Assoc_CN1.pdf", height=2,width = 14)
my_heatmap(vr, heatmap_mar = c(7,17,1,2), sub_gene = c("CDKN2A", "GNAQ", "NOTCH1", "RB1", "SMAD4", "ABL1"),hist_plot = F,legend_dist=0.2, xlab_adj = 0.9, order_omit = "NEUTRAL")
dev.off()一個函數實現基因內具有多種突變類型的熱圖的繪制