高能推荐！批量在多个组织中找出跟你的分子最相关的基因。
在第一个课程中，我们像素级别还原了Nature的一张小图，
跟Nature一起学习TCGA,GTEx和CCLE数据库的使用
更重要的是，在实现这个图之前的数据处理过程，我们强行执行了数据调整的三大步
1.基因名称转换
2.行列转置
3.添加分组信息
1.已经实现的技能
这样我们就实现了输入任意两个基因，返回其在多个组织中相关性数据的功能。
比如你看到一个分子TMED3，想看看他跟诺奖明星低氧诱导分子(HIF1A)之间的关系
只要三步就可以实现：
第一步加载数据
load(file = "pancancer_mRNA_exprSet.Rdata")
data = mRNA_exprSet
test <- data[1:10,1:10]
第二步加载函数运算
他会返回一个清洁数据框，
### 写一个函数
pancor <- function(gene1,gene2,data=data){
data1 <- split(data,data$type)
do.call(rbind,lapply(data1, function(x){
dd <- cor.test(as.numeric(x[,gene1]),as.numeric(x[,gene2]),method="pearson")
data.frame(type=x$type[1],cor=dd$estimate,p.value=dd$p.value )
}))
}
### 运算
plotdf <- pancor("TMED3","HIF1A",data)
第三步画图
pancorplot(plotdf)
图马上就出来了。
心里还是有一些希望的，因为有颜色就代表p值是小于0.05的
这两个分子至少在一些组织中是有相关性的，只是有的高(红色)有的低(蓝色)。
2.目前的困境
我想看看在我的乳腺癌中有没有相关性呢？
作图函数pancorplot支持指定肿瘤类型(如果是GTEx数据，支持的是组织类型)
pancorplot(plotdf,"BRCA")
这回彻底悲凉了，没有相关性。竖线是p值0.05。
当然，如果会编程，这些都不是个事情。我们可以让2万个基因都跟HIF1A算一次
然后返回所有的结果，让我来选。
这样就可以找到，在所有组织中相关性都很好的基因GNB1
当然，我们也可以找到在大部分组织中都负相关的分子ABHD16B
你要知道，其实总体而言，两个基因相关性为正的概率远远大于相关性为负的，这是从两万个数据中统计得出的。
这样一来，真的是天高任鸟飞啦。
3.解决方案
方法显而易见，就是写个批量的循环来做这个啊
因为我们可以算一个，就可以算多个。
修改上面的函数，让他返回一个总结后的结果，也就是一个基因一行，告诉我p值有意义的有多少个，多少是正的，多少是负的。
data1 = split(data,data$type)
### 加载函数
pancorBatch <- function(gene1,gene2,data1){
pancordata = do.call(rbind,lapply(data1, function(x){
dd <- cor.test(as.numeric(x[,gene1]),as.numeric(x[,gene2]),method="pearson")
data.frame(type=x$type[1],cor=dd$estimate,p.value=dd$p.value )
}))
subdata = pancordata[pancordata$p.value < 0.05 ,]
data.frame(gene1,
count= nrow(subdata),
positive= sum(subdata$cor > 0),
engitive= sum(subdata$cor < 0))
}
测试一下
pancorBatch("GNB1","HIF1A",data1)
返回结果：
gene1 count positive engitive
1 GNB1 33 33 0
而且我发现时间大概是0.03秒，两万个基因大概就是600秒，10分钟就搞定了。
不过我们可以用并行的方法来加快速度
library(future.apply)
plan(multiprocess)
正式运行
yourgene <- "HIF1A"
genelist <- setdiff(colnames(data),c("sample","TCGA_id","type","subtype",yourgene))
system.time(pancorfilter <- do.call(rbind,future_lapply(genelist,pancorBatch,yourgene,data1)))
运算过程中，我的小电脑被榨干了每一滴性能。
最终大概2分钟搞定。
用户 系统 流逝
28.08 47.30 139.03
返回的数据看一样吧
GBN1基因赫然在列。
我为什么没用新冠相关的分子ACE2来作为例子呢。
因为总有人要拿他来发文章，我就不展现了。
4.使用场景扩展
这样一个功能，对于任何手上有个基因的科研工作者都可以用起来。
但还有更好的场景
1.非编码RNA: 比如，你有个lncRNA，但是不知道有什么功能，可以用这里的方法来看一下。
在单个组织中我已经写过了。
单基因批量相关性分析的妙用
2.miRNA：假如你有一个miRNA在手上，那就更好了，因为miRNA和mRNA的相关性正好吻合他的作用机制。
我们可以先预测miRNA的靶基因
然后再用批量相关性去筛选负相关的。在泛组织层面得出的结论更加可靠
3.甲基化位点: 如果你能够学会数据整理的方法，把甲基化数据也整合到一起。
那么你就可以实现，任意一个甲基化位点的注释。
4.转录因子筛选: 你的基因的启动子预测出一大堆转录因子，然后用你的基因和转录因子求相关性，是不是结果更可信呢，拿出其中的3-5个去做过表达实验，只要能让你的基因表达升高或者降低，就是很好的开题，而且一定可以做出结果。
好了，这次就讲这么多，代码已经给大家展示了，
其实这里用到的就是lapply和do.call的功能，可以看看这个帖子，自己就可以写出来的。
视频小教程_R语言中的批量操作
如果在文章课题中用到的的话可以给我说一声哈，我也高兴高兴。
该部分的project文件和讲解视频会更新到产品交流群，为什么要录视频，因为有一些人迷恋我的声音(看图说话)
还有几十个没有加入的，麻烦看到联系一下我，以后所有的更新都在群里，没办法单独联系。