GEO提取lncRNA、mRNA、miRNA的表达矩阵

生信喵

背景

从GEO下载芯片数据，从基因表达矩阵分别提取lncRNA、mRNA、miRNA的表达矩阵

这套代码使用三个相对独立的模块来解决上述需求，可根据自己的实际需要，灵活使用这三个模块：

【模块一】从GEO下载microarray数据，不包括测序数据。

【模块二】根据Ensembl的BioMart对基因的biotype注释，识别出哪些基因是lncRNA、哪些是mRNA、哪些是miRNA。

不仅限于芯片数据，只要有gene symbol或Ensembl ID就能识别。

如果没有gene symbol或Ensembl ID，或其他任何跟BioMart共有的ID，就需要做序列比对，本套代码不包含序列比对功能。

【模块三】根据基因名提取对应的表达矩阵。

不仅限于GEO数据，还适用于从TCGA下载的数据、你自己的测序数据。只要提供一个基因列表和一个表达矩阵文件。

**注：**如果想要miRNA成熟体的表达矩阵，最好的方式去找miRNA芯片数据，例如GSE113596，用【模块一】下载即可。如果用【模块二】和【模块三】从普通芯片数据里提取miRNA，那其实是miRNA前体的杂交信号。

应用场景

场景一：想计算lncRNA、miRNA、mRNA间的相关性，用来找miRNA、lncRNA的靶基因，或者用功能已知的mRNA推测lncRNA的功能，首先就要获得lncRNA、miRNA跟mRNA的表达矩阵。需要从GEO下载芯片数据，就从【模块一】开始；以TCGA表达矩阵作为输入，就从【模块二】开始。

场景二：打算做RNA-seq，正在设计实验，不知道该做哪种处理、选什么时间点。那就先看看别人的实验设计获得的数据效果如何，用【模块一】下载GEO的芯片数据。

场景三：TCGA里感兴趣的癌症类型样本量太少或没有，或者想拿多个来源的数据验证，就用【模块一】下载GEO的芯片数据。

这篇为了兼顾特殊情况，文字越写越多。如果你的数据属于大多数，就不用看那么多文字，只把GSE和GPL或者“easy_input.csv”文件替换成你自己的数据，运行代码就好。

环境设置

#使用国内镜像安装包
#options("repos"= c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
#options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
install.packages("BiocManager")
library(BiocManager)
install("GEOquery")
install("biomaRt")

library(Biobase)
library(GEOquery)
library(biomaRt)

Sys.setenv(LANGUAGE = "en") #显示英文报错信息
options(stringsAsFactors = FALSE) #禁止chr转成factor

【模块一】从GEO下载基于microarray的基因表达矩阵。

如果你获得了TCGA的基因表达矩阵，已经保存为 easy_input.csv的格式，不需要从GEO下载表达矩阵，就跳过这步，直接进入【模块二】。

• 第一步：到NCBI的GEO数据库查询你感兴趣的数据，https://www.ncbi.nlm.nih.gov/geo/
• 第二步：借助GEOquery，提取表达数据。

**注意：**把GSE替换成你想用的数据，切记同时替换成这套数据所用的GPL（哪个芯片平台），才能获得正确的探针组注释信息。**GPL去哪里找？**在https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36895页面搜platform，它右边就是GPL。

gset <- getGEO("GSE36895", GSEMatrix =TRUE, getGPL = TRUE, AnnotGPL = TRUE) 
#如果提示没有GPL***.anno.gz文件，就用下面这句
#gset <- getGEO("GSE36895", GSEMatrix =TRUE, getGPL = TRUE) 

if (length(gset) > 1) idx <- grep("GPL570", attr(gset, "names")) else idx <- 1
gset <- gset[[idx]]

#查看gset里的丰富信息
#str(gset)
#如果gset里没有gene symbol或其他任何跟ensembl共有的ID，就需要用序列比对的方式去找这段探针组所对应的基因名，才能提取lncRNA。

#提取表达矩阵
exprdf<-data.frame(exprs(gset))
dim(exprdf)

#默认样品名是geo_accession，例如GSM904985
#还可以提取title作为样品名，例如Normal cortex of patient 14
#colnames(exprdf)<-gset@phenoData@data$title

#保存到文件
#write.csv(exprdf,"not_easy_input.csv",quote = F)

现在你就获得了这套芯片的表达矩阵

这时，行名是探针组的名字，不好用，我们需要给它加上基因名。

你可能会遇到三种情况：

1. 像示例所用的芯片注释文件提供了gene symbol，我们就提取gene symbol。
2. 如果你用的芯片连gene symbol都没有，就去gset和【模块二】的listAttributes里找它俩共有的ID，然后替换下面的Gene symbol。
3. 如果没有任何共有ID，就需要你自己通过序列比对去给芯片做注释，找到序列对应的基因名，【模块一】对于你来说，到这里就暂停了。

exprdf$gsym<-gset@featureData@data$`Gene symbol`
#有时要用下面这行
#exprdf$gsym<-gset@featureData@data$GENE_SYMBOL

#删除没有gene symbol的探针组
exprdf<-exprdf[exprdf$gsym!="",]
dim(exprdf) # 45118

#有的探针组对应多个基因，用“///”分隔基因名，删掉这样的行
exprdf<-exprdf[!grepl("///", exprdf$gsym),]
dim(exprdf) # 42904

#有多个探针组对应同一个基因，取中值
#如果想取平均值，就把median改为mean
exprdf_uniq<-aggregate(.~gsym,exprdf,median)
dim(exprdf_uniq) # 20848

#现在就可以用gene symbol作为行名了
rownames(exprdf_uniq)<-exprdf_uniq$gsym
#删除gene symbol列
exprdf_uniq<-subset(exprdf_uniq,select = -gsym)

#保存所有基因的表达矩阵到文件
#write.csv(exprdf_uniq,file="gene_exp.csv",row.names = T,quote = F)
#此处仅保存前4个sample
write.csv(exprdf_uniq[,1:4],file="easy_input.csv",row.names = T,quote = F)

题外话：以“GSE36895”为例，在https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36895页面底部有个“Analyze with GEO2R”按钮，点击按钮，就进入了GEO2R界面，https://www.ncbi.nlm.nih.gov/geo/geo2r/，点击Define groups，输入分组信息，然后就能做简单的差异基因筛选和画图。其实都是用R实现的，点击选项卡里的R script，就能看到R代码。

【模块二】用Ensembl的biomaRt，识别出哪些基因是lncRNA、哪些是mRNA、哪些是miRNA。

输入文件

用【模块一】获得的 easy_input.csv作为输入，基因名是gene symbol；
有的文件基因名是ensembl ID；

或者你自己写的基因列表，至少包含第一列：基因名，可以是gene symbol，或者ensembl ID。

第二列开始是表达量，非必须。

exprdf_uniq<-read.csv("easy_input.csv",header = T,row.names = 1)
rownames(exprdf_uniq)[1:4]

从biomaRt提取gene biotype

这里用到biomaRt包，来自ensembl。

• 第一步，选择你要用的基因组版本。

此处用人类ensembl最新版本，如果想用旧的基因组版本或其他物种，需要按照注释修改host = 后面的参数，

点击链接，http://asia.ensembl.org/info/website/archives/assembly.html，查看genome assembly跟ensembl版本的对应关系

#用下面这行查看ensembl基因组版本跟host的对应关系
#listEnsemblArchives()

#运行下面两行，查看基因组
#mart = useMart('ensembl')
#listDatasets(mart)
#你需要哪个物种，就复制它在dataset列里的词，放在下面这行的`dataset = `参数里
ensembl <- useMart(biomart = "ENSEMBL_MART_ENSEMBL", 
                   dataset = "hsapiens_gene_ensembl", #人
                   #dataset = "mmusculus_gene_ensembl", #小鼠
                   #dataset = "rnorvegicus_gene_ensembl", #大鼠
                   #dataset = "dmelanogaster_gene_ensembl", #果蝇
                   host = "https://www.ensembl.org") 

#植物用下面三行，以拟南芥为例
mart <- useMart(biomart = 'plants_mart', host = "https://plants.ensembl.org")
listDatasets(mart)
ensembl <- useMart(biomart = "plants_mart",
                   dataset = "athaliana_eg_gene", 
                   host = "https://plants.ensembl.org") 

生信喵

• 第二步，提取biotype

我们需要用gene symbol或ensembl ID来提取gene_biotype

然后用gene biotype区分lncRNA、miRNA和mRNA

#查看Filters和Attributes提供了哪些信息
#listFilters(ensembl)
#listAttributes(ensembl)
#看到基因名gene symbol在第64行
listAttributes(ensembl)[64,]
#下面的“filters = ”要用到这个name

feature_info <- getBM(attributes = c("gene_biotype",
                                     #"transcript_biotype",#还可以提取transcript_biotype
                                     #如果基因名是gene symbol，就运行下面这行
                                     "hgnc_symbol"), 
                                     #如果基因名是ensembl ID，就运行下面这行
                                     #"ensembl_gene_id"),
                      #如果基因名是gene symbol，就运行下面这行
                      filters = "hgnc_symbol", #小鼠是mgi_symbol，大鼠是mgi_symbol
                      #如果基因名是ensembl ID，就运行下面这行
                      #filters = "ensembl_gene_id",
                      values = rownames(exprdf_uniq), mart = ensembl)

#有些芯片注释的gene symbol跟最新版本ensembl的基因名不一致，需要返回上一步，换比较老的版本。
#TCGA数据的ensembl ID跟最新版ensembl一致
if (nrow(exprdf_uniq) != nrow(feature_info)){
  #查看哪些基因名不一致
  library(dplyr)
  diffName<-setdiff(rownames(exprdf_uniq),feature_info[,2])
  length(diffName)
  head(diffName)
}

length(unique(feature_info$hgnc_symbol))
#有些gene symbol对应多个ensembl id，因此会有多个biotype，例如
feature_info[feature_info$hgnc_symbol == "ARHGEF26-AS1",]
#TCGA数据不会遇到这个问题，因为ensembl id跟gene_biotype是一一对应的关系

#把基因的biotype保存到文件
write.csv(feature_info[,c(2,1)],"gene_biotype.csv",quote = F,row.names = F)

识别lncRNA、mRNA和miRNA

对lncRNA的定义，可参考Vega的标准：
http://vega.archive.ensembl.org/info/about/gene_and_transcript_types.html

#查看gene biotype的类型
unique(feature_info$gene_biotype)

#此处定义protein_coding作为mRNA
mRNA <-"protein_coding"
#根据实际研究目的，调整定义为lncRNA的gene_biotype，此处根据Vega定义如下8种biotype为lncRNA
lncRNA <- paste("non_coding","3prime_overlapping_ncRNA","antisense","lincRNA","sense_intronic","sense_overlapping","macro_lncRNA","bidirectional_promoter_lncRNA",sep = "|")
#还可以定义miRNA
miRNA <-"miRNA"

#下面就是表达矩阵里的mRNA、lncRNA、miRNA及其数量，把每类基因的基因名保存到相应的文件里。
mRNA.list<-feature_info[grepl(mRNA, feature_info$gene_biotype),]
write.table(mRNA.list,"mRNA.list.txt",quote = F,row.names = F, col.names = F)
nrow(mRNA.list)

lncRNA.list<-feature_info[grepl(lncRNA, feature_info$gene_biotype),]
write.table(lncRNA.list,"lncRNA.list.txt",quote = F,row.names = F, col.names = F)
nrow(lncRNA.list)

miRNA.list<-feature_info[grepl(miRNA, feature_info$gene_biotype),]
write.table(miRNA.list,"miRNA.list.txt",quote = F,row.names = F, col.names = F)
nrow(miRNA.list)

【模块三】根据基因名提取对应的表达矩阵

输入文件

两个输入文件：

• 表达矩阵：easy_input.csv。第一列：基因名。可以是gene symbol（【模块一】的输出文件），或者ensembl_gene_id。第二列开始是表达量。每行一个基因，每列一个sample。
• 基因列表：mRNA.list.txt。【模块二】的输出文件。包含一列基因名，可以是gene symbol，或者ensembl_gene_id，要跟表达矩阵的基因名一致。

exprdf_uniq<-read.csv("easy_input.csv",header = T,row.names = 1)
head(exprdf_uniq)

#以miRNA基因列表为例
miRNA.list<-read.table("miRNA.list.txt")
head(miRNA.list)

提取表达矩阵，保存到文件

mRNA_expr <- exprdf_uniq[as.character(mRNA.list[,2]),]
write.csv(mRNA_expr,"mRNA_expr.csv",quote = F,row.names = T)
head(mRNA_expr)

lncRNA_expr <- exprdf_uniq[as.character(lncRNA.list[,2]),]
write.csv(lncRNA_expr,"lncRNA_expr.csv",quote = F,row.names = T)
head(lncRNA_expr)

miRNA_expr <- exprdf_uniq[as.character(miRNA.list[,2]),]
write.csv(miRNA_expr,"miRNA_expr.csv",quote = F,row.names = T)
head(miRNA_expr)