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

1. R语言安装：官方网站/ 安装软件。2. 所需要的软件包： 2.1 affy数据处理相关的程序包在R中复制source(/biocLite.R) biocLite(affy) 2.2 热度图相关程序包 Gplots（）：install.packages(gplots)3. 获取基因表达数据 3.1 读取基因芯片数据（cel.files）the.filter - matrix(c(CEL file (*.cel), *.cel, All (*.*), *.*), ncol = 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为正则表达式查找函数(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)# 计算平均值，并做对数转换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.rma.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. 选取目的基因在/index.php上确定探针，选取数据；汇总到excel表格中，保存为csv格式。8. 热度图 cipk=read.csv(c:/users/suntao/desktop/TaCIPK affx arry log.csv) s(cipk)=cipk$genenamecipk source(/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)#必须，如果没有这句话会出错，原因不知-,-请高手指教imputed-impute.knn(as.matrix(rawexpr) ,k = 10, rowmax = 0.5, colmax = 0.8, maxp = 1500, rng.seed=)#impute.knn() 使用一个矩阵作为第一个参数，其他参数这里使用的是默认值write.table(imputed$data,file=imputed_data.txt)#write.table() 把数据保存在当前工作目录下的文件中，文件名用file= 指定，这一步不是必须的imputeddata-imputed$data#imputed$data是在R中储存imputed后的数据的矩阵现在在R里输入imputed，即填充好的数据矩阵，是不是NA值全都没了？关于impute package的详细Documentation在/packages/release/bioc/html/impute.html全部数据文件：/emanlee/R_bioconductor_genechip_data_process.zipfrom :/tag/r/用R和BioConductor进行基因芯片数据分析(三)：计算median接前一篇：/emanlee/archive/2012/12/05/.html我们已经知道要分析的数据对每个基因有3个重复测定值，经过缺失值填充后，每个基因都有3个可用值。这一步很简单，就是取这3个值的中位数，即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(imputeddata)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:/tag/r/用R和BioConductor进行基因芯片数据分析(三)：计算median接前一篇：/emanlee/archive/2012/12/05/.html我们已经知道要分析的数据对每个基因有3个重复测定值，经过缺失值填充后，每个基因都有3个可用值。这一步很简单，就是取这3个值的中位数，即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(imputeddata)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:/tag/r/用R和BioConductor进行基因芯片数据分析(五)：芯片间归一化接前一篇：用R和BioConductor进行基因芯片数据分析(四)：芯片内归一化上次进行了芯片内的归一化，但是我们的数据来自于10张芯片，为了让这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 in 1:dim(normed)1)for (j in 1:dim(normed)2)norm_logi,j-log(normedi,j)/log(2)然后利用函数进行芯片间归一化：source(/biocLite.R)biocLite(limma)library(limma)norm_log_btw_array-normalizeBetweenArrays(norm_log,method=scale)normalizeBetweenArrays()函数有许多方法，具体请看帮助。下面看看效果吧plot(density(norm_log,1),ylim=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和BioConductor进行基因芯片数据分析(六)：差异表达基因接前一篇：用R和BioConductor进行基因芯片数据分析(五)：芯片间归一化经过一系列的预处理，包括缺失值填充，中位数计算以及归一化，我们的数据终于可以用啦。下面我们就来分析一下new population和old population的个体是否有差异表达基因。判断一个基因是否差异表达有许多方法，最早使用的就是看log ratio的绝对值是否大于2,这种方法早已废弃。下一个想到的也许是t-test，诚然t-test可以统计地判断一个基因是否差异表达，但是对于有数千数万基因的芯片来说，会有很高的错误发现率（False Discovery Rate, 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).以下分析就使用Step-down maxT and minP multiple testing procedures，需要求助于Bioconductor的multtest package的mt.maxT()函数。具体用法可通过help(mt.maxT)查看。source(/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,B=)默认随机置换次数B=10000，对于microarray来说B应该比10000大很多，所以这里取B=以下是画图：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.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)adjpsPackages-“Install package from local zip file”选择package.zip文件Linux上安装R包（离线安装）:下载package.tar.gz文件在Shell终端（注意不是R）输入:sudo R CMD INSTALL package.tar.gz注意：需要sudo权限才能安装。否则提示：username is not in the sudoers file. This incident will be reported如何把sra格式转成fastq格式(fq格式)sra是NCBI 推出的存储高通量数据的格式，而平常我们工作用得多是fastq格式。如果需要把sra 转成fastq，从/Traces/sra/sra.cgi?cmd=show&f=software&m=software&s=software下载相应的软件。或者下载最新的source code，在服务器上用make 编译。然后使用如下命令行：sra_sdk-2.0.0rc1/linux/rel/gcc/x86_64/bin/fastq-dump -A SRR -D SRR.sra这样就可以很简单的把sra格式转成fastq格式了。REF:/s/blog_4055ao1mg.html/wuyu466/item/eb4363eac3baf37d29/Traces/sra/sra.cgi?view=software/s/blog_70b2bliee.html一、 R语言相关资料1、R语言/bbs/thread/2、R语言课件（复旦大学）/bbs/actions/archive/post/_1.html3、上传：GLM课件（R语言）/bbs/thread/4、利用R语言实现微阵列数据分析-聚类分析/bbs/actions/archive/post/_1.html5、Medline 文献挖掘的开放资源库-MedlineR/bbs/actions/archive/post/