本帖最后由 anlan 于 2017-12-18 19:52 编辑
这里先放有道笔记链接，下面代码看不清可以看有道笔记的 去年的这个时候，我参加了生信技能树Jimmy办的线上R语言处理表达芯片教程，这也是我进入生信领域，第一个学习的完整流程。那时对生信的一些基础概念还是处于模糊的状态，只是略懂了点R语言，勉强的听完了整个教程。 后来似乎也再也没接触过芯片数据了，对于芯片数据的概念还停留在表达矩阵那块，以致有次在找公共数据时没找到表达矩阵，只找到了芯片的raw data，却因无法分析而无能为力（现在看来是由于在学习教程时并没有完全吃透）
基础知识 一组探针（或探针组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 如果是想获取GSE27114的RAW DATA，那么可以选择使用oligo包（Affymetrix 平台），主要用于读取Affymetrix平台的基因表达数据，处理为表达矩阵；如果是Illumina平台，教程则推荐使用beadarray或lumi 首先通过download.file函数下载raw data，至于下载URL如何获取有两个途径：
数据读入及标准化 将所有样本的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
首先我们先进 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 由于通常多个探针会对应一个基因，所以我们需要在后续分析中选其中最有代表性的一个探针的表达量来代表基因的表达水平 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
[size=0em]
| 
发表于 2017-12-17 22:33:19
第二个有道云链接挂了
楼主
|手机版|小黑屋|生信技能树
( 粤ICP备15016384号 ) 