从mRNA开始的ceRNA网络构建

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的说明里有：[图片上传中…(image-3103bb-1593013367975-1)]

截图出自：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&degraExpNum=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文件，参考：[图片上传中…(image-7f516a-1593013367975-0)]

在今天的资料里这个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里完成啦。