R语言在基因芯片数据处理中的应用要点.docx

1、R语言在基因芯片数据处理中的应用要点1. R语言安装：官方网站http:/www.r-project.org/ 安装软件。2. 所需要的软件包： 2.1 affy数据处理相关的程序包在R中复制source(http:/bioconductor.org/biocLite.R) biocLite(affy) 2.2 热度图相关程序包 Gplots（）：install.packages(gplots)3. 获取基因表达数据 3.1 读取基因芯片数据（cel.files）the.filter - matrix(c(CEL file (*.cel), *.cel, All (*.*), *.*), nc

2、ol = 2, byrow = T)cel.files - choose.files(caption = Select CEL files, multi = TRUE, filters = the.filter, index = 1)raw.data - ReadAffy(filenames = cel.files) 3.2 sampleNames(raw.data)ang #先看看原样品名称的规律 3.3 pat - .*MT-(0-9A-Z+).* #样品名称查找的正则表达式 sampleNames(raw.data) - gsub(pat, 1, cel.files) #gsub为正则表

3、达式查找函数(samples - sampleNames(raw.data)pData(raw.data)$treatment - rep(c(0h, 1h, 24h, 7d), each = 2)#确定样品重复数使用rma方法进行预处理：eset.rma - rma(raw.data)# Background correcting# Normalizing# Calculating Expression4.计算基因表达量emat.rma.log2 - exprs(eset.rma)class(emat.rma.log2)head(emat.rma.log2, 1)# 计算平均值，并做对数转换

4、results.rma - data.frame(emat.rma.log2, c(1, 3, 5, 7) + emat.rma.log2, c(2, 4, 6, 8)/2)# 计算表达量差异倍数results.rma$fc.1h - results.rma, 2 - results.rma, 1results.rma$fc.24h - results.rma, 3 - results.rma, 1results.rma$fc.7d - results.rma, 4 - results.rma, 1head(results.rma, 2)5.T检验p.value = apply(emat.rm

5、a.log2,1, function(x)(t.test(x7:9, x10:12)$p.value)6.导出数据write.csv(results.rma,file=C:/users/suntao/desktop/data.csv) 7. 选取目的基因在http:/www.plexdb.org/index.php上确定探针，选取数据；汇总到excel表格中，保存为csv格式。8.热度图 cipk=read.csv(c:/users/suntao/desktop/TaCIPK affx arry log.csv) row.names(cipk)=cipk$genenamecipk source

6、(http:/bioconductor.org/biocLite.R) biocLite(impute)impute是专门用KNN法进行缺失值填充的R package:设置好当前工作目录( Windows是在R的菜单栏-文件-改变工作目录设置，Linux下用setwd()函数）然后在R控制台输入以下代码:library(impute)#导入impute packageraw-read.table(raw_data_3_replicates.txt,header=TRUE)rawexpr-raw,-1#移除第一列ID列if(exists(.Random.seed) rm(.Random.seed

7、)#必须，如果没有这句话会出错，原因不知-,-请高手指教imputed-impute.knn(as.matrix(rawexpr) ,k = 10, rowmax = 0.5, colmax = 0.8, maxp = 1500, rng.seed=362436069)#impute.knn() 使用一个矩阵作为第一个参数，其他参数这里使用的是默认值write.table(imputed$data,file=imputed_data.txt)#write.table() 把数据保存在当前工作目录下的文件中，文件名用file= 指定，这一步不是必须的imputeddata-imputed$dat

8、a#imputed$data是在R中储存imputed后的数据的矩阵现在在R里输入imputed，即填充好的数据矩阵，是不是NA值全都没了？关于impute package的详细Documentation在http:/bioconductor.fhcrc.org/packages/release/bioc/html/impute.html全部数据文件：from :用R和BioConductor进行基因芯片数据分析(三)：计算median接前一篇：我们已经知道要分析的数据对每个基因有3个重复测定值，经过缺失值填充后，每个基因都有3个可用值。这一步很简单，就是取这3个值的中位数，即median。方

9、法很多，在excel中可以用median函数;在R中以下代码进行操作：get_median-function(i,j)num_vec-c(imputeddatai*3-2,j,imputeddatai*3-1,j,imputeddatai*3,j)median(num_vec)#A simple function to calculate median value of three replicatesdimrow-(dim(imputeddata)1)/3mediandata-matrix(data = NA, nrow =dimrow, ncol = dim(imputeddata)2,

10、byrow = TRUE, dimnames = NULL)#Create a blank matrix to store median valuesfor (i in 1:dimrow)for (j in 1:dim(imputeddata)2)mediandatai,jdim(imputeddata)1 11571 20dim(mediandata)1 3857 20from:用R和BioConductor进行基因芯片数据分析(三)：计算median接前一篇：我们已经知道要分析的数据对每个基因有3个重复测定值，经过缺失值填充后，每个基因都有3个可用值。这一步很简单，就是取这3个值的中位数，

11、即median。方法很多，在excel中可以用median函数;在R中以下代码进行操作：get_median-function(i,j)num_vec-c(imputeddatai*3-2,j,imputeddatai*3-1,j,imputeddatai*3,j)median(num_vec)#A simple function to calculate median value of three replicatesdimrow-(dim(imputeddata)1)/3mediandata-matrix(data = NA, nrow =dimrow, ncol = dim(impute

12、ddata)2, byrow = TRUE, dimnames = NULL)#Create a blank matrix to store median valuesfor (i in 1:dimrow)for (j in 1:dim(imputeddata)2)mediandatai,jdim(imputeddata)1 11571 20dim(mediandata)1 3857 20from:用R和BioConductor进行基因芯片数据分析(五)：芯片间归一化接前一篇：用R和BioConductor进行基因芯片数据分析(四)：芯片内归一化上次进行了芯片内的归一化，但是我们的数据来自于1

13、0张芯片，为了让这10张芯片之间有可比性，需要进行芯片间归一化。具体原理就不介绍了。这里用到Bioconductor的一个package，叫做limma，以及其中的函数normalizeBetweenArrays()由于normalizeBetweenArrays()需要log intensity或log ratio作为输入，于是先进行log转化：#log transformationnorm_log-matrix(data = NA, nrow =dim(normed)1, ncol = dim(normed)2, byrow = TRUE, dimnames = NULL)for (i i

14、n 1:dim(normed)1)for (j in 1:dim(normed)2)norm_logi,j-log(normedi,j)/log(2)然后利用函数进行芯片间归一化：source(http:/bioconductor.org/biocLite.R)biocLite(limma)library(limma)norm_log_btw_array-normalizeBetweenArrays(norm_log,method=scale)normalizeBetweenArrays()函数有许多方法，具体请看帮助。下面看看效果吧plot(density(norm_log,1),ylim=

15、c(0,1.35),xlab=log intensity)for (i in 2:20)lines(density(norm_log,i),type=l)lines(density(norm_log_btw_array,1),type=l,col=green)for (i in 2:20)lines(density(norm_log_btw_array,i),type=l,col=green)text(1.5,c(0.8,1.0),labels=c(BEFORE normalization,AFTER normalization),col=c(black,green)用R和BioConduct

16、or进行基因芯片数据分析(六)：差异表达基因接前一篇：用R和BioConductor进行基因芯片数据分析(五)：芯片间归一化经过一系列的预处理，包括缺失值填充，中位数计算以及归一化，我们的数据终于可以用啦。下面我们就来分析一下new population和old population的个体是否有差异表达基因。判断一个基因是否差异表达有许多方法，最早使用的就是看log ratio的绝对值是否大于2,这种方法早已废弃。下一个想到的也许是t-test，诚然t-test可以统计地判断一个基因是否差异表达，但是对于有数千数万基因的芯片来说，会有很高的错误发现率（False Discovery Rate,

17、 FDR)，如果 p value 0.05,则10000个基因里有500个基因实际没有差异表达而被误认为是差异表达。因此t-test方法需要改进。于是 Westfall & Young (1993) 提出了Step-down maxT and minP multiple testing procedures，大意就是比较几个group间有没有差异基因表达，就通过随机置换这些group的标记，相当于随机互换group的成员，模拟一个空分布(null distribution），以此计算调整后的p value，这个方法可以极大的减小Family-wise Error Rate (FWER).以下分

18、析就使用Step-down maxT and minP multiple testing procedures，需要求助于Bioconductor的multtest package的mt.maxT()函数。具体用法可通过help(mt.maxT)查看。source(http:/bioconductor.org/biocLite.R)biocLite(multtest)library(multtest)classlabel-c(0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1,0,1)maxttt-mt.maxT(norm_log_btw_array,classlabel,

19、B=100000)默认随机置换次数B=10000，对于microarray来说B应该比10000大很多，所以这里取B=100000以下是画图：rawpsort(rawppp)170 1 0.0493 sort(rawppp)171 1 0.0502 170个raw p小于0.05abline(h=0.05,col=blue)text(1000,c(0.6,0.7),labels=c(raw p-value,adjusted p-value),col=c(black,red)text(1000,0.08,labels=p=0.05,col=blue)可见调整后只有一个基因的p value小于0.

20、05，而未调整的有170个基因的p value小于0.05，可以说虽然此方法降低了错误发现率，但是也导致了很高的False negative.此外可以考虑使用multtest package的mt.rawp2adjp()函数，这个函数可以通过”Bonferroni”, “Holm”, “Hochberg”, “SidakSS”, “SidakSD”, “BH”, “BY”等方法调整p value，不过对我们的数据来说都过于严格了。procs-c(Bonferroni,Holm,Hochberg,SidakSS,SidakSD,BH,BY)adjps-mt.rawp2adjp(rawp,proc

21、s)plot(sort(adjps$adjp,1),ylab=p value)for (i in 2:8)points(sort(adjps$adjp,i),col=i)abline(h=0.05,col=blue)text(1000,0.08,labels=p=0.05,col=blue)因此可以考虑不这么严格的SAM (Significance Analysis of Microarrays)分析。有兴趣的请看参考资料。参考资料：课堂讲义：Differential expression.Identifying differentially expressed genes notions on multiple testing and p-value adjustments.Dudoit, S., Yang, Y.H., Speed, T.P., and Callow, M.J. (2002),Statistical methods for identifying differentially expressed genes in replicated cDNA microarray experiments,Statistica