数据挖掘生物信息学R语言源码

这张图科研人员必备,我们来看看Cancer Cell是怎么示范的

2019-10-11  本文已影响0人  9d760c7ce737

每周一的晚上7点,是我组织的线下文献精读会。每一个月,我自己挑选8篇文章,安排人员来读,而且要求文献至少准备20小时以上(自己看懂还要想办法让别人听懂,有很多工作准备工作)。

每次是两篇,大概都是1小时左右,昨天师姐分享了一篇GBM中m6A甲基化的文献。



如果回看我们当时的那张神图,会发现这里的讲的ALKBH5是擦除子,可以抹去RNA的m6A甲基化修饰。
作者就讲了ALKBH5通过抹去FOXM1的甲基化修饰避免其被降解,从而促进FOXM1表达的故事。



这篇文章使用了很多公共数据,我们下一次再讲。今天我要讲的是,高分文章也爱用基于TCGA数据的相关性分析,只不过本次有点蹊跷。

作者用免疫荧光发现ALKBH5和SOX2和nestin共定位。一个蛋白如果要发挥作用,要关注三个点,表达量,活性,位置。现在他们的位置靠在一起,有了相互作用的空间优势,那么表达量是否相关呢?作者就用了TCGA的GBM 转录组数据来测试了一下,发现效果不错。


你看这些公共数据放在那里,不用白不用,用了还有奇效。因为群里已经有三位朋友说,把Nature的那个图用在了自己的课题上。



但是,公共数据,你能用,别人也可以用,通过学习那个课程,我们可以轻易地自己画出相关性系数的图来测试一下。
加载数据

rm(list=ls())
load(file = "pancancer_mRNA_exprSet.Rdata")
data <- mRNA_exprSet
test <- data[1:10,1:10]

我们为了方便,把单个癌症中的相关性作图搞成了一个函数,ggcorplot

ggcorplot <- function(a,b,type="ALL"){
  corr_eqn <- function(x,y,digits=2) {
    test <- cor.test(x,y,type="pearson")
    paste(paste0("n = ",length(x)),
          paste0("r = ",round(test$estimate,digits),"(pearson)"),
          paste0("p.value= ",round(test$p.value,digits)),
          sep = ", ")
  }
  if(type=="ALL"){
    plot_df <- data[,c(a,b)]
  }else{
    plot_df <- data[data$type %in% type,c(a,b)]
  }
  
  names(plot_df) <- c("geneA","geneB")
  ggplot(plot_df,aes(geneA,geneB))+
    geom_point(col="#984ea3")+
    geom_smooth(method=lm, se=T,na.rm=T, fullrange=T,size=2,col="#fdc086")+
    geom_rug(col="#7fc97f")+
    theme_minimal()+
    xlab(paste0("Expression of ",a," (log2(TPM))"))+
    ylab(paste0("Expression of ",b," (log2(TPM))"))+
    labs(title = paste0(corr_eqn(plot_df$geneA,plot_df$geneB)))+
    theme(plot.title = element_text(hjust = 0.5),
          plot.margin = margin(1, 1, 1, 1, "cm"))
}

输入两个基因,不指定癌症就默认是全部癌症,如果指定就是特定癌症,我们来看看ALKBH5和SOX2在GBM中的表达

ggcorplot("ALKBH5","SOX2","GBM")

正相关,但是数据量比较少,我想到,胶质瘤在TCGA中有两种,还有一个是LGG,我们再试一下

ggcorplot("ALKBH5","SOX2","LGG")

没有相关性

我们也不能太武断,他的数据有可能是合并了GBM和LGG的,因为GEPIA就是这么做的,我们用GEPIA来测试一下




那如果我们把自己的这个数据合并一下呢,也很简单,多谢几个癌症名称就可以了

ggcorplot("ALKBH5","SOX2",c("GBM","LGG"))

结果差不多,但是我们现在知道这个结果谁的贡献比较大(如果单独下载TCGA的数据来分析,这个相关性系数是负数)。

如果用到我们课程里面的图,TCGA的这个数据可以这样表示。



这里我更新了代码,之前的竖线画在了0.5的位置,是我笔误,应该画在0.05。
把自己感兴趣的点标注出来,群里的孔老师昨天就实现了。我在这里让小于0.05的癌症全部变暗。
改变了x轴和y轴,留出空间给有可能增大的点。同时x轴起始设置在0.1,那么有三个点就不能出现在图上。所以我把起始设置在0.01开始,这样避免一些点不能显示在图上。这样的话,0.05这个竖线就会往中间偏移,视觉上更容易区分四个象限。

这些更新会在V2版本中更新给大家。
最后,我想说,这个技能在我心里重千斤,未来,1000个用户都是少的,因为真的是人人可用啊。

最后,课程正在拼团打折中。再次宣传。


如果不清楚这个图有多重要,就看这篇
跟Nature一起学习TCGA,GTEx和CCLE数据库的使用
我们还给这个技能拓展过1次
教程拓展:手上在研究的基因在各种组织,癌症,细胞系中的表达量。

在这之前,我还写过以下几个相关性的帖子,都可以看一看。
单基因批量相关性分析的妙用
又是神器!基于单基因批量相关性分析的GSEA
多个基因的相关性如何分析与展示?
未来,这些帖子都会在当前这个技能的加持下变得异常强大,敬请期待。

上一篇下一篇

猜你喜欢

热点阅读