画相关性分析图的韦恩图部分:
使用VennDiagram包
最普通的功能是venn.diagram()
,输入的是向量(最多五个集合),会自动生成韦恩图,有许多参数可以选择。也有较详细的功能,如两个集合比较、三个集合比较的功能。与venn.diagram
不同的是,draw.pairwise.venn
(两个集合比较)功能输入的参数是集合1的数量、集合2的数量以及共有元素的数量。在相关性分析中,我们想找出差异表达和差异甲基化基因的相关性,所以属于两个集合比较。除了画图,还需要让用户知道是哪些基因是overlap的,所以可以用专门针对两个集合比较的功能,顺便找出overlap的基因。准备输入的数据是之前算出来的有差异表达和差异甲基化的genelist,命名为de_genelist和dm_genelist。为了找出它们共有的基因,要运用到intersect()
函数,这个函数是取两个集合之间的交集的,集合中可以是数值、字符串等。拓展:并集union()
、找不同setdiff(x,y)
取x中与y不同的元素、判断相同setequal()
只有完全相同的时候才返回TRUE。
先试一下这个函数的功能:
venn.plot <- draw.pairwise.venn(100, 70, 30, c("First", "Second"));
grid.draw(venn.plot);
grid.newpage();
结果如图:
可以看到基本的信息都有了,但是可以把它做得美观一些
venn.plot <- draw.pairwise.venn(area1 = 100,
area2 = 70,
cross.area = 30,
category = c("First", "Second"),
fill = c("#F29B05","#A1D490"),
ext.text = TRUE,
ext.percent = c(0.1,0.1,0.1),
ext.length = 0.6,
label.col = rep("gray10",3),
lwd = 0,
cex = 2,
lty = "blank",
alpha = rep(0.3, 2),
fontface = rep("bold",3),
fontfamily = rep("sans",3),
cat.cex = 1.5,
cat.fontface = rep("plain",2),
cat.fontfamily = rep("sans",2),
cat.pos = c(0, 0),
print.mode = c("raw","percent"))
grid.draw(venn.plot)
经过对参数的调整,这样就比较美观了。
进入正题:
所以取两个基因集之间的交集是:
overlap <- intersect(de_genelist, dm_genelist)
画韦恩图:
venn.plot <- draw.pairwise.venn(area1 = length(de_genelist),
area2 = length(dm_genelist),
cross.area = length(overlap),
category = c("Different Expression Genes", "Different Methylation Genes"), #集合的名字
fill = c("#F29B05","#A1D490"), #集合填充的颜色
ext.text = TRUE, #是否在韦恩图的某一部分太小时把其中标签放在外面
ext.percent = c(0.1,0.1,0.1), #在某一部分小于这个值时才把标签放在外面
ext.length = 0.6, #连接外置标签和其指示部分的线长比例
label.col = rep("gray10",3), #每部分标签的颜色
lwd = 0, #集合外围的线的宽度
cex = 2, #标签字体大小
lty = "blank", #集合外围的线的样式,这里是无
alpha = rep(0.3, 2), #透明度,默认是0.5
fontface = rep("bold",3), #标签的字体样式,这里是加粗
fontfamily = rep("sans",3), #标签的字体,默认timesnewroman
cat.cex = 1.5, #集合名字的大小
cat.fontface = rep("plain",2), #集合名字的样式
cat.fontfamily = rep("sans",2), #集合名字的字体
cat.pos = c(0, 0), #集合名字的位置
print.mode = c("raw","percent")) #标签以数字展示还是百分比展示,这里是上数字下百分比
grid.draw(venn.plot) #画图
完成
散点图:
目标是做到上面的图
cor_plot <- function(CoxExpPlotData,gene1,gene2,cormethod="spearman"){
library("ggplot2")
# cor_df = data.frame(CoxExpPlotData,
# gene1 = CoxExpPlotData[,gene1],
# gene2 = CoxExpPlotData[,gene2])
pcutoff=0.05
CoxTest = cor.test(CoxExpPlotData[, gene1], CoxExpPlotData[, gene2], method = cormethod)
if(cormethod == "pearson"){
plottitle <- paste0(cormethod,", R = ", signif(CoxTest$estimate[[1]], 3), "\nP value = ", signif(CoxTest$p.value, 4), sep = "")
}else if(cormethod == "spearman"){
plottitle <- paste0(cormethod,", rho = ", signif(CoxTest$estimate[[1]], 3), "\nP value = ", signif(CoxTest$p.value, 4), sep = "")
}else if (cormethod == "kendall"){
plottitle <- paste0(cormethod,", tau = ", signif(CoxTest$estimate[[1]], 3), "\nP value = ", signif(CoxTest$p.value, 4), sep = "")
}
print(paste("pvalue =", round(CoxTest$p.value, 4)))
if(CoxTest$p.value <= pcutoff){
CoxExpPoint = ggplot(data = CoxExpPlotData, aes_string(x= gene1, y = gene2))+theme_classic()+
geom_point(size = 2, color = "gray36")+
annotate("text",x=-Inf,y=Inf,label=plottitle,hjust=-.2,vjust=2)+
stat_smooth(method = lm, se = FALSE, colour = "#191970")+
ylab(paste(gene2, "Exp.", sep = " "))+xlab(paste(gene1, "Exp.", sep = " "))+
scale_y_continuous(expand = c(0, 0))+scale_x_continuous(expand = c(0, 0))+
theme(plot.title = element_text(size = 15, angle = 0, face = "plain", colour = "black", hjust = 0.5, vjust = -0.5),
axis.title.x = element_text(size = 15, angle = 0, face = "plain", colour = "black"),
axis.title.y = element_text(size = 15, angle = 90, face = "plain", colour = "black"),
axis.text.x = element_text(size = 15, angle = 0, face = "plain", colour = "black"),
axis.text.y = element_text(size = 15, angle = 0, face = "plain", colour = "black"))}
return(CoxExpPoint)
}