前面提到過gdcRNAtools里面的ceRNA網(wǎng)絡(luò)構(gòu)建
構(gòu)建鑒定ceRNA的標(biāo)準(zhǔn)有4個(gè)：
（1）lncRNA和mRNA之間共享的miRNA數(shù)量篡撵；
（2）lncRNA與mRNA的表達(dá)相關(guān)性追迟；
（3）共享miRNA對lncRNA和mRNA的調(diào)控相似性；
（4）(sensitivity)靈敏度相關(guān)性
前兩條在很多ceRNA網(wǎng)絡(luò)構(gòu)建的文章中提及，比如：
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6308055/
https://pubmed.ncbi.nlm.nih.gov/29082457/
當(dāng)然也有大把的文章里沒有提及這兩個(gè)檢驗(yàn)酗失，做了是錦上添花姐霍，不做也沒有太大問題。
落到具體的代碼實(shí)現(xiàn)盯荤，就是計(jì)算一個(gè)是幾何分布檢驗(yàn)的p值馋吗，一個(gè)皮爾遜相關(guān)性p值。（我知道很多人要被名字嚇跑了廷雅。但是我還是要說耗美，別走京髓！不難！）
gdcRNAtools里構(gòu)建網(wǎng)絡(luò)的函數(shù)需要將miRNA的表達(dá)矩陣一起輸入商架，今天做的數(shù)據(jù)是沒有miRNA表達(dá)矩陣的堰怨，其實(shí)不影響網(wǎng)絡(luò)構(gòu)建，只是不能用這個(gè)一步到位的函數(shù)啦蛇摸。
超幾何分布檢驗(yàn)就是反應(yīng)二者交集足夠多不
關(guān)于超幾何分布檢驗(yàn)备图，放進(jìn)R語言也無非就是一個(gè)函數(shù)。有一個(gè)講的比較好的例子赶袄，參考自：https://blog.csdn.net/linkequa/article/details/88189582
設(shè)總共有29個(gè)人揽涮，其中11個(gè)吸煙者，18個(gè)非吸煙者饿肺，現(xiàn)從中隨機(jī)抽取16個(gè)樣本（在此實(shí)驗(yàn)中對應(yīng)著肺癌病人）蒋困，有10個(gè)是吸煙者，這樣的事件是否顯著敬辣？
p_value=phyper(10-1, #10個(gè)是吸煙者
11, #11個(gè)吸煙者
18, #18個(gè)非吸煙者
16, #隨機(jī)抽取16個(gè)(癌癥人數(shù))
lower.tail=F)
p_value
## [1] 0.003135274
設(shè)總共有29個(gè)人雪标，其中16個(gè)癌癥患者，13個(gè)正常人溉跃，現(xiàn)從中隨機(jī)選出11個(gè)人（對應(yīng)著所有抽煙的人）村刨，有10個(gè)是癌癥患者，這樣的事件是否顯著撰茎？
p_value=phyper(10-1, #10個(gè)是癌癥患者
16, #16個(gè)癌癥患者
13, #13個(gè)正常人
11, #隨機(jī)選出11個(gè)人(吸煙人數(shù))
lower.tail=F)
p_value
## [1] 0.003135274
emmm......想繼續(xù)寫點(diǎn)啥的嵌牺，算了。
只需要學(xué)個(gè)小函數(shù)
我去找了gdcRNAtools 在github上的代碼龄糊，發(fā)現(xiàn)超幾何分布檢驗(yàn)其實(shí)被寫成函數(shù)了逆粹，只是沒有export，因此只是相當(dāng)于作者寫的時(shí)候分段绎签，而用戶無法使用枯饿。
所以我將這個(gè)函數(shù)拿下來自己捯飭了一下，放進(jìn)了tinyarray诡必，下載1.3.3以上的版本可以用奢方。使用起來很簡單的！０质妗蟋字！
1.輸入數(shù)據(jù)
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是從數(shù)據(jù)庫獲取的lncRNA與miRNA對應(yīng)關(guān)系，m_mi 是從數(shù)據(jù)庫獲取的mRNA和miRNA的對應(yīng)關(guān)系扭勉。注意miRNA都應(yīng)放在第二列的位置上鹊奖，列名可以隨意。
2.超幾何分布檢驗(yàn)
hyperoutput = hypertest(
lnc = rownames(lnc_exp),
pc = rownames(mRNA_exp),
lnctarget = lnc_mi,
pctarget = m_mi
)
hyperPValue 是超幾何分布檢驗(yàn)的p值涂炎，Counts是共享miRNA的數(shù)量忠聚，可以根據(jù)這兩個(gè)指標(biāo)來篩選设哗。
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.相關(guān)性檢驗(yàn)
直接選擇符合超幾何分布檢驗(yàn)的基因和lncRNA來做，減少計(jì)算時(shí)間两蟀。plcortest函數(shù)网梢，直接返回相關(guān)性p值<0.05的基因，以列表的形式組織赂毯，列表元素名字為lncRNA战虏，元素內(nèi)容為mRNA。
hylnc = unique(hysmall$lncRNAs)
hypc = unique(hysmall$Genes)
v = plcortest(lnc_exp[hylnc,],
mRNA_exp[hypc,])
length(v) #有多少個(gè)lncRNA找到了顯著相關(guān)的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" ...
相關(guān)性顯著與超幾何分布檢驗(yàn)顯著的lncRNA就是v列表元素的名字党涕。而兩個(gè)都顯著的mRNA烦感，就需要取交集來實(shí)現(xiàn)了。
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"
最后組合起來膛堤，成為一個(gè)數(shù)據(jù)框手趣。可以統(tǒng)計(jì)經(jīng)過這兩種檢驗(yàn)的三種RNA的個(gè)數(shù)了骑祟。
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包里現(xiàn)有的函數(shù)有不少啦回懦，還會(huì)繼續(xù)更新：
geo_download() : 提供geo編號，返回表達(dá)矩陣次企、臨床信息表格和使用的平臺編號。
get_deg() ：提供芯片表達(dá)矩陣潜圃、分組信息缸棵、探針注釋，返回差異分析結(jié)果谭期。
multi_deg() : 多個(gè)分組（最多5個(gè)）的差異分析
cor.full()和cor.one() :批量計(jì)算基因間的相關(guān)性
幾個(gè)繪圖函數(shù)：draw_heatmap,draw_volcano,draw_venn
trans_exp():將tcga或tcga+gtex數(shù)據(jù)進(jìn)行基因id轉(zhuǎn)換
t_choose():批量做單個(gè)基因的t檢驗(yàn)
point_cut():批量計(jì)算生存分析最佳截點(diǎn)
surv_KM():批量做KM生存分析堵第，支持用最佳截點(diǎn)分組
surv_cox():批量做cox生存分析，支持用最佳截點(diǎn)分組
hypertest():批量做mRNA和lncRNA的超幾何分布檢驗(yàn)
plcortest():批量做mRNA和lncRNA的相關(guān)性檢驗(yàn)
順便說一句隧出，我上課時(shí)會(huì)和學(xué)員說：你運(yùn)行一個(gè)代碼報(bào)錯(cuò)了踏志，不代表代碼錯(cuò)了，在這里再補(bǔ)一句胀瞪，你裝不上一個(gè)包针余，也不代表這個(gè)包的作者寫錯(cuò)了。解決安裝R包的報(bào)錯(cuò)凄诞，是所有R語言學(xué)習(xí)者的必修課圆雁，會(huì)遇到很多坑，網(wǎng)絡(luò)問題占大多數(shù)帆谍。github上的包可以下載下來伪朽，使用devtools::install_local進(jìn)行本地安裝，但仍然不能確保一次成功汛蝙，需要解決問題的能力烈涮。
