GEO提取lncRNA、mRNA、miRNA的表达矩阵
<h1>背景</h1><p>从GEO下载芯片数据，从基因表达矩阵分别提取lncRNA、mRNA、miRNA的表达矩阵</p>
<p>这套代码使用三个相对独立的模块来解决上述需求，可根据自己的实际需要，灵活使用这三个模块：</p>
<p>【模块一】从GEO下载microarray数据，不包括测序数据。</p>
<p>【模块二】根据Ensembl的BioMart对基因的biotype注释，识别出哪些基因是lncRNA、哪些是mRNA、哪些是miRNA。</p>
<p>不仅限于芯片数据，只要有gene symbol或Ensembl ID就能识别。</p>
<p>如果没有gene symbol或Ensembl ID，或其他任何跟BioMart共有的ID，就需要做序列比对，本套代码不包含序列比对功能。</p>
<p>【模块三】根据基因名提取对应的表达矩阵。</p>
<p>不仅限于GEO数据，还适用于从TCGA下载的数据、你自己的测序数据。只要提供一个基因列表和一个表达矩阵文件。</p>
<p>**注：**如果想要miRNA成熟体的表达矩阵，最好的方式去找miRNA芯片数据，例如GSE113596，用【模块一】下载即可。如果用【模块二】和【模块三】从普通芯片数据里提取miRNA，那其实是miRNA前体的杂交信号。</p>
<h1>应用场景</h1>
<p>场景一：想计算lncRNA、miRNA、mRNA间的相关性，用来找miRNA、lncRNA的靶基因，或者用功能已知的mRNA推测lncRNA的功能，首先就要获得lncRNA、miRNA跟mRNA的表达矩阵。需要从GEO下载芯片数据，就从【模块一】开始；以TCGA表达矩阵作为输入，就从【模块二】开始。</p>
<p>场景二：打算做RNA-seq，正在设计实验，不知道该做哪种处理、选什么时间点。那就先看看别人的实验设计获得的数据效果如何，用【模块一】下载GEO的芯片数据。</p>
<p>场景三：TCGA里感兴趣的癌症类型样本量太少或没有，或者想拿多个来源的数据验证，就用【模块一】下载GEO的芯片数据。</p>
<blockquote>
<p>这篇为了兼顾特殊情况，文字越写越多。如果你的数据属于大多数，就不用看那么多文字，只把GSE和GPL或者“easy_input.csv”文件替换成你自己的数据，运行代码就好。</p>
</blockquote>
<h1>环境设置</h1>
<pre><code class="language-r">#使用国内镜像安装包
#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")
</code></pre>
<pre><code class="language-{r}">library(Biobase)
library(GEOquery)
library(biomaRt)
Sys.setenv(LANGUAGE = "en") #显示英文报错信息
options(stringsAsFactors = FALSE) #禁止chr转成factor
</code></pre>
<h1>【模块一】从GEO下载基于microarray的基因表达矩阵。</h1>
<p>如果你获得了TCGA的基因表达矩阵，已经保存为 <code>easy_input.csv</code>的格式，不需要从GEO下载表达矩阵，就跳过这步，直接进入【模块二】。</p>
<ul>
<li>第一步：到NCBI的GEO数据库查询你感兴趣的数据，https://www.ncbi.nlm.nih.gov/geo/</li>
<li>第二步：借助GEOquery，提取表达数据。</li>
</ul>
<p>**注意：**把GSE替换成你想用的数据，切记同时替换成这套数据所用的GPL（哪个芯片平台），才能获得正确的探针组注释信息。**GPL去哪里找？**在https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36895页面搜platform，它右边就是GPL。</p>
<pre><code class="language-{r,message=FALSE}">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[]
#查看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)
</code></pre>
<p>现在你就获得了这套芯片的表达矩阵</p>
<p>这时，行名是探针组的名字，不好用，我们需要给它加上基因名。</p>
<p>你可能会遇到三种情况：</p>
<ol>
<li>像示例所用的芯片注释文件提供了gene symbol，我们就提取gene symbol。</li>
<li>如果你用的芯片连gene symbol都没有，就去gset和【模块二】的listAttributes里找它俩共有的ID，然后替换下面的Gene symbol。</li>
<li>如果没有任何共有ID，就需要你自己通过序列比对去给芯片做注释，找到序列对应的基因名，【模块一】对于你来说，到这里就暂停了。</li>
</ol>
<pre><code class="language-{r}">exprdf$gsym<-gset@featureData@data$`Gene symbol`
#有时要用下面这行
#exprdf$gsym<-gset@featureData@data$GENE_SYMBOL
#删除没有gene symbol的探针组
exprdf<-exprdf
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)
</code></pre>
<p><strong>题外话：<strong>以“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，输入分组信息，然后就能做</strong>简单的差异基因筛选和画图</strong>。其实都是用R实现的，点击选项卡里的R script，就能看到R代码。</p>
<h1>【模块二】用Ensembl的biomaRt，识别出哪些基因是lncRNA、哪些是mRNA、哪些是miRNA。</h1>
<h2>输入文件</h2>
<p>用【模块一】获得的 <code>easy_input.csv</code>作为输入，基因名是gene symbol；<br />
有的文件基因名是ensembl ID；</p>
<p>或者你自己写的基因列表，至少包含第一列：基因名，可以是gene symbol，或者ensembl ID。</p>
<p>第二列开始是表达量，非必须。</p>
<pre><code class="language-{r}">exprdf_uniq<-read.csv("easy_input.csv",header = T,row.names = 1)
rownames(exprdf_uniq)
</code></pre>
<h2>从biomaRt提取gene biotype</h2>
<p>这里用到biomaRt包，来自ensembl。</p>
<ul>
<li>第一步，选择你要用的基因组版本。</li>
</ul>
<p>此处用人类ensembl最新版本，如果想用旧的基因组版本或其他物种，需要按照注释修改host = 后面的参数，</p>
<p>点击链接，http://asia.ensembl.org/info/website/archives/assembly.html，查看genome assembly跟ensembl版本的对应关系</p>
<pre><code class="language-{r}">#用下面这行查看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")
</code></pre>
<pre><code class="language-r">#植物用下面三行，以拟南芥为例
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")
</code></pre>
<ul>
<li>第二步，提取biotype</li>
</ul>
<p>我们需要用gene symbol或ensembl ID来提取gene_biotype</p>
<p>然后用gene biotype区分lncRNA、miRNA和mRNA</p>
<pre><code class="language-{r,message=FALSE}">#查看Filters和Attributes提供了哪些信息
#listFilters(ensembl)
#listAttributes(ensembl)
#看到基因名gene symbol在第64行
listAttributes(ensembl)
#下面的“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
#TCGA数据不会遇到这个问题，因为ensembl id跟gene_biotype是一一对应的关系
#把基因的biotype保存到文件
write.csv(feature_info[,c(2,1)],"gene_biotype.csv",quote = F,row.names = F)
</code></pre>
<h2>识别lncRNA、mRNA和miRNA</h2>
<p>对lncRNA的定义，可参考Vega的标准：<br />
http://vega.archive.ensembl.org/info/about/gene_and_transcript_types.html</p>
<pre><code class="language-{r}">#查看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
write.table(mRNA.list,"mRNA.list.txt",quote = F,row.names = F, col.names = F)
nrow(mRNA.list)
lncRNA.list<-feature_info
write.table(lncRNA.list,"lncRNA.list.txt",quote = F,row.names = F, col.names = F)
nrow(lncRNA.list)
miRNA.list<-feature_info
write.table(miRNA.list,"miRNA.list.txt",quote = F,row.names = F, col.names = F)
nrow(miRNA.list)
</code></pre>
<h1>【模块三】根据基因名提取对应的表达矩阵</h1>
<h2>输入文件</h2>
<p>两个输入文件：</p>
<ul>
<li>表达矩阵：easy_input.csv。第一列：基因名。可以是gene symbol（【模块一】的输出文件），或者ensembl_gene_id。第二列开始是表达量。每行一个基因，每列一个sample。</li>
<li>基因列表：mRNA.list.txt。【模块二】的输出文件。包含一列基因名，可以是gene symbol，或者ensembl_gene_id，要跟表达矩阵的基因名一致。</li>
</ul>
<pre><code class="language-{r}">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)
</code></pre>
<h3>提取表达矩阵，保存到文件</h3>
<pre><code class="language-{r}">mRNA_expr <- exprdf_uniq),]
write.csv(mRNA_expr,"mRNA_expr.csv",quote = F,row.names = T)
head(mRNA_expr)
lncRNA_expr <- exprdf_uniq),]
write.csv(lncRNA_expr,"lncRNA_expr.csv",quote = F,row.names = T)
head(lncRNA_expr)
miRNA_expr <- exprdf_uniq),]
write.csv(miRNA_expr,"miRNA_expr.csv",quote = F,row.names = T)
head(miRNA_expr)
</code></pre>
页:
[1]