2019-04-01

画相关性分析图的韦恩图部分：
使用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();

结果如图：

image.png

可以看到基本的信息都有了，但是可以把它做得美观一些

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)

image.png

经过对参数的调整，这样就比较美观了。

进入正题：
所以取两个基因集之间的交集是：

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) #画图

完成
散点图：

image.png

目标是做到上面的图

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)
}