搜索
查看: 3993|回复: 2

回顾基因表达芯片分析

[复制链接]

29

主题

131

帖子

1208

积分

金牌会员

Rank: 6Rank: 6

积分
1208
发表于 2017-12-17 22:33:19 | 显示全部楼层 |阅读模式
本帖最后由 anlan 于 2017-12-18 19:52 编辑

这里先放有道笔记链接,下面代码看不清可以看有道笔记的
去年的这个时候,我参加了生信技能树Jimmy办的线上R语言处理表达芯片教程,这也是我进入生信领域,第一个学习的完整流程。那时对生信的一些基础概念还是处于模糊的状态,只是略懂了点R语言,勉强的听完了整个教程。
后来似乎也再也没接触过芯片数据了,对于芯片数据的概念还停留在表达矩阵那块,以致有次在找公共数据时没找到表达矩阵,只找到了芯片的raw data,却因无法分析而无能为力(现在看来是由于在学习教程时并没有完全吃透
有幸不久之前看到一篇推荐的芯片教程Differential gene expression analysis I (microarray data),看了后感觉受益不少,然后再从头看了一遍Jimmy的教程,粗略了解Affymetrix平台芯片该如何进行前期预处理以及获取表达矩阵

基础知识
一组探针(或探针组probe set)来自于一个基因,通常是由11对探针组成,每对探针是由匹配探针(PM)和错配探针(MM)组成,叫做探针对(Probe Pair),样本RNA则叫做靶标。芯片中MM探针作用是检测非特异杂交信号,理论上,MM只有非特异杂交信号,而不会有特异性杂交。MM的信号值永远小于其对应的PM信号值,因此可以做一个PM-MM或PM/MM即可去除背景噪声的影响

数据下载
以芯片教程中的GSE27114作为实例进行后续操作。在GSE27114里有两个study,一个是比较旧的Affymetrix平台数据(GPL1261):Affymetrix Mouse Genome 430 2.0 Array(Affymetrix 3’ end arrays),另一个是相对较新的平台数据(GPL6246):Affymetrix Mouse Gene 1.0 ST Array(Affymetrix arrays Whole Transcript Gene Arrays),下面先以GPL1261为例:
在数据下载之前还需要先了解自己需要哪一个芯片数据,也就是说是需要芯片原始数据(CEL文件)还是表达矩阵数据(matrix)
如果是只想获得芯片数据已经做完标准化等预处理过的表达矩阵,那么可以使用GEOquery包的getGEO函数,表达矩阵文件的压缩包也会下载到R的工作目录下
library(GEOquery) gset <- getGEO("GSE27114", destdir = ".")
由于GSE27114里有2个study,所以我先看看有哪两个study
names(gset) ##[1] "GSE27114-GPL1261_series_matrix.txt.gz" "GSE27114-GPL6246_series_matrix.txt.gz"
由于下面会以GPL1261为例,所以我先取出第一个study,可以看到gset_one的数据是储存在ExpressionSet对象中,这个对象里包含了N多的信息,可以通过slotNames列出可以查看的信息,然后根据列出来的函数一一查看,比如assayData用于保存表达数据信息,experimentData用于保存实验设计信息
gset_one <- gset[[1]] class(gset_one) ##[1] "ExpressionSet" ##attr(,"package") ##[1] "Biobase"
接下来就可以通过exprs函数获取表达矩阵,有45101个探针,6个样本
exprSet <- exprs(gset_one) dim(exprSet) ##[1] 45101     6
还可以通过pData获取这个study的更多详细信息,比如样本信息、分组信息、rawdata信息,比如这个例子就有32列信息。。。
pdata <- pData(gset_one) dim(pdata) ##[1]  6 32
当然如果只想获取表达矩阵,并且其他信息并不想知晓(或者已经在GEO网站上查到),那么通过GEO网站手动下载表达矩阵也是可以的:点击GSE27114网页下方的Series Matrix File(s)即可进入ftp下载界面ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE27nnn/GSE27114/matrix/,而且这个ftp的URL可以根据你的GSE号自行进行修改,再配合脚本就可以批量下载N个GSE样本的表达矩阵,很好用。
如果是想获取GSE27114的RAW DATA,那么可以选择使用oligo包(Affymetrix 平台),主要用于读取Affymetrix平台的基因表达数据,处理为表达矩阵;如果是Illumina平台,教程则推荐使用beadarray或lumi
首先通过download.file函数下载raw data,至于下载URL如何获取有两个途径:
  • 通过上述的pdata里查看,有一列是各个样本的rawdata下载地址,然后写个循环即可批量下载数据
    pdata$supplementary_file
  • 通过GSE27114网站的最下方GSE27114_RAW.tar的,点击http即可下载所有文件的压缩包(网速慢或者不稳定可以首选这个,至少下载一个算一个,慢慢来。。。)


数据读入及标准化
将所有样本的raw data文件放置在R工作目录下(Microarray文件夹),然后用list.celfiles函数查看所有文件,再用read.celfiles函数读入R中,并将上述的pdata信息也一并写入
library(oligo) file_CELs <- list.celfiles("Microarray", listGzipped = TRUE, full.name = TRUE) rawAffy <- read.celfiles(filenames = file_CELs, phenoData = phenoData(gset_one), sampleNames = rownames(pData(gset_one)))
注:读入文件时会检查R中是否安装了该数据对应平台的注释包,比如这个例子需要这个pd.mouse430.2包
粗略看下,rawAffy数据中包含了芯片上所有探针的强度值信息
head(exprs(rawAffy), n = 2)
我们可以拿这个数据来进行质控,但是一般质控会在芯片公司那边就做了,所以我们拿公共数据的用的话,一般也不用进行质控了。所以我们接下来就对rawAffy数据进行标准化处理,因为每个样本芯片的绝对光密度是不一样的,不做处理的话无法将各个样本间进行数据分析,同时也是为了消除测量间的非实验差异
芯片常用的标准化方法也有好几种(如:MAS5),这里就不一一尝试了,就拿最常用的RMA(Robust Multichip Average algorithm),引用芯片教程,RMA标准化过程主要分为3步:
  • Background correction (removes array auto-fluorescence ) 背景矫正
    • Quantile normalization (makes all intensity distributions identical) 标准化
    • Probeset summarization (calculates one representative value per probeset) 汇总

由于rawAffy是ExpressionFeatureSet对象,则rma函数默认是probeset level
class(rawAffy) ##[1] "ExpressionFeatureSet" ##attr(,"package") ##[1] "oligoClasses" eset <- rma(rawAffy)
我们可以看看标准化后的结果

最后我们需要得到表达矩阵,也是用exprs函数,从结果上来看,跟之前getGEO得到的表达矩阵是一样的(说明用的是同一种标准化方法?以及相同的QC处理?)
exprSet <- exprs(eset) dim(exprSet) ##[1] 45101     6
这次以R语言处理表达芯片教程中的GSE号为例子,具体可以参看http://www.biotrainee.com/thread-566-1-1.html
首先我们先进GSE42872网站看一下这个芯片数据是什么平台的,如:Affymetrix Human Gene 1.0 ST Array [transcript (gene) version],这说明这个芯片是Affymetrix的一个human表达谱芯片,用来检测基因表达情况。
接下来就之前文章回顾基因表达芯片分析(数据处理)中同样的方法先对数据进行处理
用getGEO下载表达矩阵,其实这里不下载也没事,因为后续数据分析还是会从raw data开始,再到标准化后的表达矩阵。用slotNames等函数主要还是为了看下gset[[1]]这个对象的一些信息
gset <- getGEO("GSE42872",getGPL = FALSE) class(gset[[1]]) gset[[1]] slotNames(gset[[1]]) annotation(gset[[1]]) phenoData(gset[[1]]) featureData(gset[[1]]) assayData(gset[[1]])
然后获取表达矩阵和比较感兴趣的pdata
head(featureNames(gset[[1]])) pdata <- pData(gset[[1]]) exprSet_ori <- exprs(gset[[1]]) dim(exprSet_ori)
下面开始则是比较重要的步骤,从芯片的raw data数据开始到最后的差异表达分析

数据标准化
由于网速不太好,我还是选择人工下载再放置在data文件中,然后用read.celfiles读入
file_CELs <- list.celfiles("data", listGzipped = TRUE, full.name = TRUE) rawAffy <- read.celfiles(filenames = file_CELs, phenoData = phenoData(gset[[1]]), sampleNames = rownames(pData(gset[[1]])))
标准化这步跟之前的略有不同,区别在于我们要获取转录本的表达量,而不是probeset;所以在使用rma函数时略有不同(其实rma默认参数就是基于转录本的)
eset <- rma(rawAffy, target = "core") show(eset) dim(eset) featureData(eset) <- getNetAffx(eset, "transcript")
Select one probeset per gene
由于通常多个探针会对应一个基因,所以我们需要在后续分析中选其中最有代表性的一个探针的表达量来代表基因的表达水平
Differential gene expression analysis I (microarray data)教程中所介绍是使用IQR值最大的方法,也就是选择在所有样本中表达量变化最大的那个探针。还有一种方法(平时听得较多)是选用表达量最大的那个探针。不管哪种标准,其处理方法其实都是类似的,比如以IQR最大的方法为例,先计算每个探针的IQR值
iqrs <- apply(exprs(eset), 1, IQR)
然后使用genefilter包的findLargest函数,其实这步自己来写R代码来实现也可以的,嫌麻烦就用这个函数,并且这时还需要一个芯片对应的注释R包。比如这个芯片是Affymetrix Human Gene 1.0 ST Array,那么需要下载的注释包则是hugene10sttranscriptcluster.db;如果是小鼠的 Affymetrix Mouse Gene 1.0 ST Array,则对应的注释包为mogene10sttranscriptcluster.db
library(hugene10sttranscriptcluster.db) library(genefilter) prbs <- findLargest(featureNames(eset), testStat = iqrs, data = "hugene10sttranscriptcluster.db") eset2 <- eset[prbs, ] dim(eset2)
然后获取表达矩阵
exprSet <- exprs(eset)
如果想选择取探针中表达量最大的方法,只需要将iqrs向量变成每个探针在所有样本中表达量的平均值(或者总和)即可

Probe ID转化为symbol
这步是分析芯片数据中必不可少的,比如我们可以看下exprSet表达矩阵的探针ID
head(exprSet, 2) ##GSM1052615 GSM1052616 GSM1052617 GSM1052618 GSM1052619 GSM1052620 ##8144866   5.429810   5.628468   5.376701   6.054502   6.106508   5.975622 ##8066431   9.543157   9.406544   9.500876  10.150883  10.251234  10.006384
这个8144866数字就代表着一个探针,这对后续下游分析是无法识别的,所以要将其转化为更有识别性的ID,比如gene symbol,这里也需要用到芯片对应的注释R包
我们可以看下hugene10sttranscriptcluster.db这个包有哪些调用数据的函数
columns(hugene10sttranscriptcluster.db) ##[1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS"  ##[6] "ENTREZID"     "ENZYME"       "EVIDENCE"     "EVIDENCEALL"  "GENENAME"     ##[11] "GO"           "GOALL"        "IPI"          "MAP"          "OMIM"         ##[16] "ONTOLOGY"     "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"         ##[21] "PROBEID"      "PROSITE"      "REFSEQ"       "SYMBOL"       "UCSCKG"       ##[26] "UNIGENE"      "UNIPROT"
使用SYMBOL从注释包中提取出探针ID与symbol对应关系的数据框
probe2symbol <- toTable(hugene10sttranscriptclusterSYMBOL)
然后将exprSet中的探针ID转化为symbol
exprSet <- as.data.frame(exprSet) exprSet$probe_id <- rownames(exprSet) exprSet_symbol <- merge(probe2symbol, exprSet, by = "probe_id") dim(exprSet_symbol)
再把exprSet中symbol列转化为列名,并只保留各个样本表达量的列(去除掉其他不必要的列)
rownames(exprSet_symbol) <- exprSet_symbol$symbol  exprSet_symbol <- exprSet_symbol[, -c(1,2)]
差异表达分析
芯片数据的差异表达分析一般首选的就是limma包。当然limma现在不仅可以处理芯片数据,而且可以处理RNA-seq数据。如本文的例子,是两组样本单因素实验设计,代码就比较简单了
先设置分组信息(分组信息可以看pdata)
library(limma) condition <- factor(c(rep("control", 3), rep("vemurafenib", 3)), levels = c("control", "vemurafenib")) design <- model.matrix(~condition)
然后常规的差异分析,这里默认是使用limma-trend方法(如果是RNA-seq数据还可以使用voom方法)
fit <- lmFit(log2(exprSet_symbol), design) fit=eBayes(fit) output <- topTable(fit, coef=2,n=Inf)
最后统计下差异基因数目
sum(output$adj.P.Val<0.01) ##[1] 5274
GO/KEGG富集分析以及PPI
后续的GO/KEGG分析的方法跟RNA-seq数据,具体可参照之前的RNA-seq分析的博文转录组差异表达分析小实战(二)
PPI分析可以使用在线的STRING网站,也可以使用stringDB包,具体可以参照之前的一篇博文STRINGDB包的简单使用

[size=0em]​




上一篇:ChIP-seq基础入门学习
下一篇:转录组
回复

使用道具 举报

1

主题

43

帖子

463

积分

中级会员

Rank: 3Rank: 3

积分
463
发表于 2017-12-18 12:53:24 | 显示全部楼层
第二个有道云链接挂了
人生若只如初见!
回复 支持 反对

使用道具 举报

29

主题

131

帖子

1208

积分

金牌会员

Rank: 6Rank: 6

积分
1208
 楼主| 发表于 2017-12-18 19:54:10 | 显示全部楼层
y461650833y 发表于 2017-12-18 12:53
第二个有道云链接挂了

谢谢提醒,以重新修改链接~
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2019-9-22 06:26 , Processed in 0.042211 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.