##输入文件:
#01矩阵,每行一个基因,每列一个sample。01矩阵代表基因在这个sample里是否变异。
#
#pathway.txt,基因所在的pathway;
#
# input_cli.csv,信息,例如人种、病史、治疗方式等。
#
# 不同文件之间,基因名、sampleID必须一致。
library(ComplexHeatmap)
#读取基因变异数据
mygene <- read.csv("1matrix01.csv", header = T,row.names = 1, as.is = T)
#mygene<-read.table("clipboard",header = T,row.names = 1,as.is = T)
#查看部分数据的格式
mygene[1:3,1:4]
#用"mut"替换1,用""替换0
mygene[mygene == 1] <- "mut"
mygene[mygene == 0] <- ""
#读取基因所在的pathway
mypathway <- read.csv("2pathway.csv", header = T, row.names = 1, as.is = T)
head(mypathway)
#把pathway添加到mygene的最后一列
mygene$pathway <- mypathway[as.character(rownames(mygene)),]
#读取临床信息
mytype <- read.csv("3clib.csv")
head(mytype)
mytype$V2
#查看有多少种临床信息类型
unique(mytype$type_A)
unique(mytype$type_B)
unique(mytype$type_C)
##设置变异的颜色、形状
#用fill = 设置mutation以及背景用什么颜色
alter_fun = list(
background = function(x, y, w, h) {
grid.rect(x, y, w-unit(0.5, "mm"), h-unit(0.5, "mm"), gp = gpar(fill = NA, col = NA)) #不要背景色
},
mut = function(x, y, w, h) {
grid.rect(x, y, w-unit(0.5, "mm"), h-unit(0.5, "mm"), gp = gpar(fill = "#008000", col = NA)) #mut是绿色
})
#mutation bar plot的颜色,跟瀑布图一致
col = c("mut" = "#008000")
#定义足够多的颜色
mycol <- c("#223D6C","#D20A13","#FFD121","#088247","#11AA4D","#58CDD9","#7A142C","#5D90BA","#431A3D","#91612D","#6E568C","#E0367A","#D8D155","#64495D","#7CC767")
###从临床数据中删除没有变异的sample
#参数remove_empty_columns = TRUE会删掉没有变异的sample,我们需要把临床数据中相应的sample也删掉。
#为了提取瀑布图中的sample,修改了oncoprint函数,保存为oncoPrint_plus.R,搜索#Ya#查看修改的位置。
#注意:在后面所有画图命令中,column_order和remove_empty_columns这两个参数,都要跟p1的参数一致。
#提取瀑布图中画出的sample
svg("1.onco.svg",width = 30,height = 15)
p1 <- oncoPrint(mygene[1:(ncol(mygene)-1)], get_type = function(x) x,
alter_fun = alter_fun, col = col,
remove_empty_columns = TRUE, #删除没有突变的sample
#column_order = NULL, #不按突变频率给sample排序
#row_order = NULL, #不按突变频率给基因排序
row_split = mygene$path)
p1
dev.off()
matrix <- p1@matrix
sampleOrder <- data.frame(p1@column_order)
rownames(sampleOrder) <- p1@column_names_param$labels
sampleOrder$oriOrder <- row.names(sampleOrder)
sampleOrder <- sampleOrder[order(as.numeric(sampleOrder[,1]),decreasing=F),]
rownames(mytype) <- mytype$sample
#在临床数据中,只保留瀑布图中画出的sample
mytype <- mytype[sampleOrder$oriOrder,]
mytype
##现在,这个mytype就跟瀑布图中的sample一致了。
mytype[3]
#画临床数据的heatmap
my_annotation = HeatmapAnnotation(df = data.frame(mytype[3]),
col = list(Chr = c("chr01" = mycol[5],"chr06" = mycol[2], "chr02" = mycol[4], "chr11" = mycol[1],"chr07" = mycol[3],"chr08" = mycol[6])))
library(ggplot2)
mytype[2]
#画瀑布图
png("2.lubby.png")
p2 <- oncoPrint(mygene[1:(ncol(mygene)-1)], get_type = function(x) x,
alter_fun = alter_fun, col = col,
remove_empty_columns = TRUE,#删除没有突变的sample
#column_order = NULL, #不按突变频率给sample排序
#row_order = NULL, #不按突变频率给基因排序
show_pct = FALSE, #左侧不显示百分比
bottom_annotation = my_annotation,#把临床信息画在下面
row_split = mygene$path, #按照pathway分开画
show_heatmap_legend = FALSE) #不显示突变的图例
p2
dev.off()
我放一张结果图
只能放一半
网友评论