一个函数实现基因内具有多种突变类型的热图的绘制

??我们平常多见的基因突变热图是一个基因一个格子，一种突变类型，但实际上在同一个病人中，同一个基因往往具有多种突变类型，因此传统的热图绘制工具并不能满足我们绘图的需要。应研究需要，本人自己写了一个热图绘制函数，内部调用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()

原文地址：https://www.cnblogs.com/ywliao/p/8991346.html

时间： 2024-07-31 20:41:34

一个函数实现基因内具有多种突变类型的热图的绘制的相关文章

【C语言】编写一个函数，传入a,b两个int类型的变量，返回两个值的最大公约数。

/*编写一个函数,传入a,b两个int类型的变量,返回两个值的最大公约数. 例如:输入传入(0 , 5)函数返回5,传入(10 , 9)函数返回1,传入(12 , 4)函数返回4 */ #include <stdio.h> int yue(int a,int b) { int temp; int n; if (a>b) { temp=a; a=b; b=temp; } n=a; if(a==0) return b; else while(n!=0) { if( a%n==0 &&

自己写一个swap函数交换任意两个相同类型元素的值 对空指针的使用 字节大小的判断（二）了解原理

验证的代码: 1 #include <stdio.h> 2 int main(){ 3 4 char c = 'z'; 5 int num = ('Z' << 24) + (c << 16) + ('A' << 8) + 'a'; 6 7 printf("'Z'=0x%x c=0x%x 'A'=0x%x 'a'=0x%x\n", 'Z', c, 'A', 'a'); 8 printf("num=0x%x\n", num

Winform开发框架中实现同时兼容多种数据库类型处理

在很多应用系统里面,虽然一般采用一种数据库运行,但是由于各种情况的需要,可能业务系统会部署在不同类型的数据库上,如果开发的系统能够很方便支持多种数据库的切换,那可以为我们减少很多烦恼,同时提高系统的适应性和强壮型.还有一种情况,由于业务数据库的不断膨胀或者方便数据库的切割隔离,有时候也会把不同的业务数据库进行分拆,如权限提供数据库,客户关系管理数据库,工作流程数据库,企业营运数据库等等,因此在一个系统里面,同时使用2个或者以上的数据库的情况也是有的. 在我较早期的一篇随笔<Winform开发框架

Winform开发框架中实现多种数据库类型切换以及分拆数据库的支持 - 伍华聪

在很多应用系统里面,虽然一般采用一种数据库运行,但是由于各种情况的需要,可能业务系统会部署在不同类型的数据库上,如果开发的系统能够很方便支持多种数据库的切换,那可以为我们减少很多烦恼,同时提高系统的适应性和强壮型.还有一种情况,由于业务数据库的不断膨胀或者方便数据库的切割隔离,有时候也会把不同的业务数据库进行分拆,如权限提供数据库,客户关系管理数据库,工作流程数据库,企业营运数据库等等,因此在一个系统里面,同时使用2个或者以上的数据库的情况也是有的.针对这两种情况,本文介绍在我的Winform开

利用接口实现多种数据库类型的灵活更换

当存在可能要更换数据库类型的时候,要考虑两个问题: 一,不同类型的数据库命名空间不一样,用到的函数名也不一样,尽管很相似: 二,有些SQL语句在不一样的数据库之间是不通用的! 那么要在更换数据库类型的时候,如何做到尽量少受因为上面两点而造成的影响呢?!利用接口,可以将第一点的影响降到最低!至于第二点, 可以将数据的操作尽量在数据库服务器端实现,这点的探讨不在本文研究范围内,在此也就不多说!下面主要就利用接口将第一点的影响降到最低这一点进行详细说 明. 首先,明确两点: 1, 无论哪种数据库类型,

实现一个函数clone，使JavaScript中的5种主要的数据类型（包括Number、String、Object、Array、Boolean）进行值复制

实现一个函数clone,可以对JavaScript中的5种主要的数据类型(包括Number.String.Object.Array.Boolean)进行值复制. 1 /** 对象克隆 2 * 支持基本数据类型及对象 3 * 递归方法 */ 4 function clone(obj) { 5 var o; 6 switch (typeof obj) { 7 case "undefined": 8 break; 9 case "string": o = obj + &q

一个函数只能返回一个值？笑话

函数的define:若有两个变量x和y,对于x的每一个值,y都有唯一确定的值与它对应,则y与x有函数关系.一般用表示.其中x叫做自变量,y叫做因变量. 唯一的值?对啊,C语言函数中也是只能return 一个的而已,虽然你能写多个return 到C语言的函数中,可是见到第一个return 后,后面的操作都是不执行的.读者可以试试!!那么,我们怎样从一个C语言函数中得到多个值呢?我们可以用些技术吧???对的,真的可以!!!①可以利用指针,指针里放着变量的地址,如果我们改变它的指向,不就可以改变它的值

一定间隔时间下重复执行一个函数的几个方法

如果有个操作,我们需要过一会儿再做,或者每隔一段时间就要做一次.可以有很多种做法. 独立线程 是的,对.NET Framework本身一知半解的程序员才会使用这种方案.不过,现实中这个方案其实并不少见. public static void Repeat(this Action action, TimeSpan interval) { new Thread(new ThreadStart(() => { while(true) { Thread.Sleep(interval); action()

编写一个函数，可以左旋字符串中k个字符

题目: 实现一个函数,可以左旋字符串中的k个字符. 例如: abcdef左旋一个字符得到bcdefa abcdef左旋两个字符得到cdefab 题目分析: 对于这个问题,可以用很多种方法求解,这里介绍两种方法: 算法一: 左旋字符串的k个字符,我们可以先将剩下的n-k个字符移动最前面,然后将左旋的k个字符移动到字符串的最后面的位置,这就完成了字符串左旋k个字符的功能.其中,对于n-k个字符的移动,可以利用while循环,将后面的字符逐一进行拷贝.下面是具体的程序: #define _CRT_SE