ceRNA 最近被玩的很多，构建的方法很多，这个是我这几天探索觉得比较好用也比较省力的，分享给大家。
我已经把代码 数据 文件打包好了，在生信星球公众号后台回复ce657即可获得。
0.R包安装
options(BioC_mirror="http://mirrors.tuna.tsinghua.edu.cn/bioconductor/")
options("repos" = c(CRAN="http://mirrors.cloud.tencent.com/CRAN/"))
options(download.file.method = 'libcurl')
options(url.method='libcurl')
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
if(!require(multiMiR))BiocManager::install("multiMiR",ask = F,update = F)
library(multiMiR)
1.从mRNA得到miRNA
输入数据的获得方式：使用某个癌症的RNA-seq数据经过差异分析、PPI网络等的筛选，得到几个关键的mRNA，保存在9hubgenes.txt。
x = read.table("9hubgenes.txt",stringsAsFactors = F)$V1;x
## [1] "MMP9" "CXCL8" "ACTB" "TGB1" "STAT1" "TOP2A" "CDK1" "GNMT" "ABAT"
我翻阅了很多网页工具和教程，发现multiMiR这个R包很优秀，结合了14个数据库，其中包括了有实验方法验证互作关系的mirTarbase。详见：[microRNAs靶基因数据库哪家强]https://mp.weixin.qq.com/s/n_UncYeGIQFLneTMK2rTXQ
这个工具可以从mRNA得到miRNA，也可以从miRNA得到mRNA，在这里使用前者，table = ’validated’这个参数是默认的，写上是为了和table = ’predicted’区分开，两者对应的数据库源不同，前者是经实验验证的，数量会更少一些，也更可靠一些。
gene2mir <- get_multimir(org = 'hsa',
target = x,
table = 'validated',
summary = TRUE,
predicted.cutoff.type = 'n',
predicted.cutoff = 500000)
## Searching mirecords ...
## Searching mirtarbase ...
## Searching tarbase ...
ez = gene2mir@data[gene2mir@data$database=="mirtarbase",];dim(ez)
## [1] 397 10
这里大材小用一下，只选了3个数据库中的一个，实验验证也分几个等级,最严谨的是Luciferase reporter assay，只选出它：
table(ez$support_type)
##
## Functional MTI Functional MTI (Weak) Non-Functional MTI
## 43 353 1
ez = ez[stringr::str_detect(ez$experiment,
"Luciferase reporter assay"),];dim(ez)
## [1] 25 10
miRNAs = unique(ez$mature_mirna_id)
2.从miRNA得到lncRNA
我查了一下相关的文献，lncRNA 和miRNA互作的数据库使用比较多的有三个：mircode，starbase，mirnet，我都看了一下，最后感觉starbase表现最好。
我在网页上戳戳戳了半天，发现starbase只能一次搜索一个miRNA，我疯了。想要下载它的lncRNA - miRNA 互作数据自己探索，没有在网页上找到直接下载的按钮，但找到了关于API的说明里有：
截图出自：http://starbase.sysu.edu.cn/tutorialAPI.php
在linux命令行用curl，写对筛选要求就可以获取数据啦。全部lncRNA - miRNA 互作数据的获取代码是：
curl 'http://starbase.sysu.edu.cn/api/miRNATarget/?assembly=hg19&geneType=lncRNA&miRNA=all&clipExpNum=0°raExpNum=0&pancancerNum=0&programNum=0&program=None&target=all&cellType=all' > starBaseV3_hg19_CLIP-seq_lncRNA_all.txt &
下载的这句代码来自：https://www.jianshu.com/p/b7e4830c0b01
有了这个数据，读进R语言就可以随便玩耍啦。
starbase = data.table::fread("starBaseV3_hg19_CLIP-seq_lncRNA_all.txt");dim(starbase)
## [1] 71952 17
很多文献里会把GENECODE里没有注释的lncRNA去掉，这个也很好操作，anno.Rdata来自于genecodev22版本的gtf文件，参考：https://mp.weixin.qq.com/s/bGoUbLuBdPteo-oG8ckMVw
在今天的资料里这个anno.Rdata文件直接提供，可以反复使用的。
load("anno.Rdata")
p1 = starbase$geneName %in% lnc_anno$gene_name;table(p1)
## p1
## FALSE TRUE
## 43924 28028
starbase = starbase[p1,];dim(starbase)
## [1] 28028 17
lnc_mi = starbase[starbase$miRNAname %in% miRNAs,]
length(unique(lnc_mi$miRNAname))
## [1] 19
colnames(lnc_mi)
## [1] "miRNAid" "miRNAname" "geneID" "geneName" "geneType"
## [6] "chromosome" "start" "end" "strand" "clipExpNum"
## [11] "degraExpNum" "RBP" "merClass" "miRseq" "align"
## [16] "targetSeq" "pancancerNum"
对starbase数据库提供的数据列名的解释里，有两个比较重要的：
clipExpNum ：The number of CLIP-seq experiments ; pancancerNum : Number of Cancer types (Pan-Cancer) that miRNA-target show anti-correlation relationships (pearson correlation: r<0, p-value<0.05).
所以pancerNum是几 就意味着这对miRNA-lncRNA在多少种癌症中表达量负相关，省掉了很多计算。
数据库的tutorial里面还提到了一个比较严格的筛选标准： CLIP evidence (>=5), degradome evidence (>=1), Pan-Cancer (>=10), program number (>=5) and predicted program (None).
degradome evidence 的限制条件，加上和不加数量有差别。
p2 = lnc_mi$pancancerNum >10 &lnc_mi$clipExpNum>4;table(p2)
## p2
## FALSE TRUE
## 927 15
p3 = lnc_mi$pancancerNum >10 &lnc_mi$clipExpNum>4 & lnc_mi$degraExpNum >0;table(p3)
## p3
## FALSE TRUE
## 939 3
lnc_mi$geneName[p2]
## [1] "MALAT1" "MALAT1" "SNHG1" "KCNQ1OT1" "ZFAS1"
## [6] "GAS5" "FGD5-AS1" "PITPNA-AS1" "LRRC75A-AS1" "LRRC75A-AS1"
## [11] "H19" "OIP5-AS1" "TUG1" "TUG1" "HAGLR"
lnc_mi$geneName[p3]
## [1] "SNHG1" "PITPNA-AS1" "LRRC75A-AS1"
后面是网络的可视化，可以在cytoscape里完成啦。
网友评论