两个检验给ceRNA锦上添花
前面提到过gdcRNAtools里面的ceRNA网络构建
构建鉴定ceRNA的标准有4个：
（1）lncRNA和mRNA之间共享的miRNA数量；
（2）lncRNA与mRNA的表达相关性；
（3）共享miRNA对lncRNA和mRNA的调控相似性；
（4）(sensitivity)灵敏度相关性
前两条在很多ceRNA网络构建的文章中提及，比如：
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6308055/
https://pubmed.ncbi.nlm.nih.gov/29082457/
当然也有大把的文章里没有提及这两个检验，做了是锦上添花，不做也没有太大问题。
落到具体的代码实现，就是计算一个是几何分布检验的p值，一个皮尔逊相关性p值。（我知道很多人要被名字吓跑了。但是我还是要说，别走！不难！）
gdcRNAtools里构建网络的函数需要将miRNA的表达矩阵一起输入，今天做的数据是没有miRNA表达矩阵的，其实不影响网络构建，只是不能用这个一步到位的函数啦。
超几何分布检验就是反应二者交集足够多不
关于超几何分布检验，放进R语言也无非就是一个函数。有一个讲的比较好的例子，参考自：https://blog.csdn.net/linkequa/article/details/88189582
设总共有29个人，其中11个吸烟者，18个非吸烟者，现从中随机抽取16个样本（在此实验中对应着肺癌病人），有10个是吸烟者，这样的事件是否显著？
p_value=phyper(10-1, #10个是吸烟者
11, #11个吸烟者
18, #18个非吸烟者
16, #随机抽取16个(癌症人数)
lower.tail=F)
p_value
## [1] 0.003135274
设总共有29个人，其中16个癌症患者，13个正常人，现从中随机选出11个人（对应着所有抽烟的人），有10个是癌症患者，这样的事件是否显著？
p_value=phyper(10-1, #10个是癌症患者
16, #16个癌症患者
13, #13个正常人
11, #随机选出11个人(吸烟人数)
lower.tail=F)
p_value
## [1] 0.003135274
emmm......想继续写点啥的，算了。
只需要学个小函数
我去找了gdcRNAtools 在github上的代码，发现超几何分布检验其实被写成函数了，只是没有export，因此只是相当于作者写的时候分段，而用户无法使用。
所以我将这个函数拿下来自己捯饬了一下，放进了tinyarray，下载1.3.3以上的版本可以用。使用起来很简单的！！！
1.输入数据
rm(list = ls())
library(tinyarray)
load("exp_target.Rdata")
ls()
## [1] "lnc_exp" "lnc_mi" "m_mi" "mRNA_exp"
lnc_exp[1:4,1:4]
## TCGA-BR-A4QL-01
## LINC00337 4.0000
## CDC42-IT1 2.5850
## LINC00853 5.2854
## ROR1-AS1 1.5850
## GTEX-QLQ7-0826-SM-447B3
## LINC00337 1
## CDC42-IT1 0
## LINC00853 0
## ROR1-AS1 1
## TCGA-VQ-A927-01
## LINC00337 3.8074
## CDC42-IT1 2.5850
## LINC00853 3.7004
## ROR1-AS1 4.1699
## GTEX-132AR-2426-SM-5IFFD
## LINC00337 0
## CDC42-IT1 2
## LINC00853 1
## ROR1-AS1 2
mRNA_exp[1:4,1:4]
## TCGA-BR-A4QL-01 GTEX-QLQ7-0826-SM-447B3
## CFAP74 5.4919 0.0000
## GABRD 4.3219 2.8074
## TP73 9.9092 3.3291
## HES2 10.6812 0.0000
## TCGA-VQ-A927-01 GTEX-132AR-2426-SM-5IFFD
## CFAP74 2.3219 2.0000
## GABRD 6.0444 2.8074
## TP73 7.8184 4.9773
## HES2 6.8826 2.0000
head(lnc_mi)
## geneName miRNAname
## 1211 MALAT1 hsa-miR-1-3p
## 1212 MALAT1 hsa-miR-1-3p
## 1258 UCA1 hsa-miR-1-3p
## 1303 SND1-IT1 hsa-miR-101-3p
## 1316 MALAT1 hsa-miR-101-3p
## 1336 TYMSOS hsa-miR-101-3p
head(m_mi)
## target_symbol mature_mirna_id
## 125 ACHE hsa-miR-212-3p
## 128 ADRA2A hsa-miR-26b-5p
## 129 ACAN hsa-miR-181a-5p
## 130 AGT hsa-miR-26b-5p
## 132 ALDH3B2 hsa-miR-124-3p
## 141 ALOX15 hsa-miR-665
lnc_mi是从数据库获取的lncRNA与miRNA对应关系，m_mi 是从数据库获取的mRNA和miRNA的对应关系。注意miRNA都应放在第二列的位置上，列名可以随意。
2.超几何分布检验
hyperoutput = hypertest(
lnc = rownames(lnc_exp),
pc = rownames(mRNA_exp),
lnctarget = lnc_mi,
pctarget = m_mi
)
hyperPValue 是超几何分布检验的p值，Counts是共享miRNA的数量，可以根据这两个指标来筛选。
k1 = hyperoutput$hyperPValue<0.05;table(k1)
## k1
## FALSE TRUE
## 3112 510
k2 = hyperoutput$Counts>2;table(k2)
## k2
## FALSE TRUE
## 2691 931
table(k1&k2)
##
## FALSE TRUE
## 3337 285
hysmall = hyperoutput[k1&k2,];dim(hysmall)
## [1] 285 9
3.相关性检验
直接选择符合超几何分布检验的基因和lncRNA来做，减少计算时间。plcortest函数，直接返回相关性p值<0.05的基因，以列表的形式组织，列表元素名字为lncRNA，元素内容为mRNA。
hylnc = unique(hysmall$lncRNAs)
hypc = unique(hysmall$Genes)
v = plcortest(lnc_exp[hylnc,],
mRNA_exp[hypc,])
length(v) #有多少个lncRNA找到了显著相关的mRNA
## [1] 19
str(v[1:3])
## List of 3
## $ TRPM2-AS : chr [1:169] "NUSAP1" "APOBEC3B" "GPRC5A" "CXCL5" ...
## $ TEX41 : chr [1:169] "NUSAP1" "APOBEC3B" "GPRC5A" "CXCL5" ...
## $ LINC00944: chr [1:169] "NUSAP1" "APOBEC3B" "GPRC5A" "CXCL5" ...
相关性显著与超几何分布检验显著的lncRNA就是v列表元素的名字。而两个都显著的mRNA，就需要取交集来实现了。
hc = lapply(hylnc,function(x){
intersect(hysmall[hysmall$lncRNAs==x,"Genes"],
v[[x]])
})
names(hc) = hylnc
hc = hc[sapply(hc,length)!=0]
str(hc[1:3])
## List of 3
## $ TRPM2-AS : chr [1:78] "NUSAP1" "GPRC5A" "ASPM" "MMP9" ...
## $ TEX41 : chr [1:77] "APOBEC3B" "CXCL5" "KIF18B" "KIF14" ...
## $ LINC00944: chr [1:3] "KLK10" "CREG2" "CCR6"
最后组合起来，成为一个数据框。可以统计经过这两种检验的三种RNA的个数了。
ce = hysmall[hysmall$lncRNAs %in% names(hc) & hysmall$Genes %in% unlist(hc),
c(1,2,ncol(hysmall))]
head(ce,3)
## lncRNAs Genes
## 17746 TRPM2-AS NUSAP1
## 3780 TEX41 APOBEC3B
## 17701 TRPM2-AS GPRC5A
## miRNAs
## 17746 hsa-miR-103a-3p,hsa-miR-107,hsa-miR-138-5p,hsa-miR-16-5p,hsa-miR-195-5p
## 3780 hsa-miR-103a-3p,hsa-miR-107,hsa-miR-16-5p,hsa-miR-195-5p
## 17701 hsa-miR-103a-3p,hsa-miR-107,hsa-miR-15a-5p,hsa-miR-15b-5p,hsa-miR-16-5p,hsa-miR-195-5p,hsa-miR-424-5p,hsa-miR-497-5p,hsa-miR-6838-5p
length(unique(ce$lncRNAs))
## [1] 18
length(unique(ce$Genes))
## [1] 169
length(unique(unlist(stringr::str_split(ce$miRNAs,","))))
## [1] 135
tinyarray包里现有的函数有不少啦，还会继续更新：
geo_download() : 提供geo编号，返回表达矩阵、临床信息表格和使用的平台编号。
get_deg() ：提供芯片表达矩阵、分组信息、探针注释，返回差异分析结果。
multi_deg() : 多个分组（最多5个）的差异分析
cor.full()和cor.one() :批量计算基因间的相关性
几个绘图函数：draw_heatmap,draw_volcano,draw_venn
trans_exp():将tcga或tcga+gtex数据进行基因id转换
t_choose():批量做单个基因的t检验
point_cut():批量计算生存分析最佳截点
surv_KM():批量做KM生存分析，支持用最佳截点分组
surv_cox():批量做cox生存分析，支持用最佳截点分组
hypertest():批量做mRNA和lncRNA的超几何分布检验
plcortest():批量做mRNA和lncRNA的相关性检验
顺便说一句，我上课时会和学员说：你运行一个代码报错了，不代表代码错了，在这里再补一句，你装不上一个包，也不代表这个包的作者写错了。解决安装R包的报错，是所有R语言学习者的必修课，会遇到很多坑，网络问题占大多数。github上的包可以下载下来，使用devtools::install_local进行本地安装，但仍然不能确保一次成功，需要解决问题的能力。
