R语言在基因芯片数据处理中的应用要点.docx
《R语言在基因芯片数据处理中的应用要点.docx》由会员分享,可在线阅读,更多相关《R语言在基因芯片数据处理中的应用要点.docx(14页珍藏版)》请在冰豆网上搜索。
![R语言在基因芯片数据处理中的应用要点.docx](https://file1.bdocx.com/fileroot1/2023-1/7/7b6766b2-eb83-442e-aa60-be49a8cb2476/7b6766b2-eb83-442e-aa60-be49a8cb24761.gif)
R语言在基因芯片数据处理中的应用要点
1.R语言安装:
官方网站 http:
//www.r-project.org/安装软件。
2.所需要的软件包:
2.1affy数据处理相关的程序包
在R中复制source("http:
//bioconductor.org/biocLite.R")
biocLite("affy")
2.2热度图相关程序包
Gplots():
install.packages("gplots")
3.获取基因表达数据
3.1读取基因芯片数据(cel.files)
the.filter<-matrix(c("CELfile(*.cel)","*.cel","All(*.*)","*.*"),ncol=2,byrow=T)
cel.files<-choose.files(caption="SelectCELfiles",multi=TRUE,filters=the.filter,index=1)
raw.data<-ReadAffy(filenames=cel.files)
3.2sampleNames(raw.data)ang#先看看原样品名称的规律
3.3pat<-".*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)
##Backgroundcorrecting
##Normalizing
##CalculatingExpression
4.计算基因表达量
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[,1]
results.rma$fc.24h<-results.rma[,3]-results.rma[,1]
results.rma$fc.7d<-results.rma[,4]-results.rma[,1]
head(results.rma,2)
5.T检验
p.value=apply(emat.rma.log2,1,function(x)(t.test(x[7:
9],x[10:
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/TaCIPKaffxarrylog.csv")
row.names(cipk)=cipk$genename
cipk<-cipk[,-1]
cipk_matrix=data.matrix(cipk)
library(gplots)
heatmap.2(cipk_matrix,Rowv=FALSE,Colv=FALSE,col=greenred(75),key=TRUE,keysize=0.8,trace="none",density.info="none",symkey=FALSE,revC=FALSE,margins=c(10,10),denscol=tracecol,distfun=dist,hclustfun=hclust,dendrogram="none",symm=FALSE)
heatmap.2颜色选择函数col=colorRampPalette(c("black","red"))
用R和BioConductor进行基因芯片数据分析
(二):
缺失值填充
以下分析用到的数据可以在这里( )下载,这个数据来自关于基因对蝴蝶迁移性的研究,样本是20个蝴蝶个体,其中10个是当地固有个体(old),另外10个是新迁入的个体(new),old和new个体两两随机配对,分别用不同颜色染料(波长分别为555和647nm)标记后,在同一张基因芯片上杂交;此外,每个基因在每张芯片上都重复点样3次,因此此数据是有3个replicates及10张芯片的双通道芯片。
数据是样点的信号强度值,没有经过标准化处理的。
拿到数据你会看到许多”NA”,这是因为我把缺失的空白值替换成NA了,以便用R进行缺失值填充。
说到缺失值填充,通常有3种方法:
A.用此基因的平均表达值填充;如果有多张重复芯片,可以取不同芯片上的平均值;对于时间序列芯片,可以通过插值法填充。
此方法很简单,也比较常用,但是效果不及下面2种方法
B.基于SVD(即单值分解)方法的填充:
简单地讲,此方法是通过描述基因表达谱的几个基本模式来对缺失值进行填充。
C.基于KNN(最近邻)方法的填充:
此方法是寻找和有缺失值的基因的表达谱相似的其他基因,通过这些基因的表达值(依照表达谱相似性加权)来填充缺失值。
KNN法是这3种方法里效果最好的,因此对本数据的缺失值用KNN法填充。
对以上3种方法的比较,这篇paper提供了清晰的说明:
Troyanskaya,O.,Cantor,M.,Sherlock,G.,Brown,P.,Hastie,T.,Tibshirani,R.,Botstein,D.,andAltman,R.B.(2001), MissingvalueestimationmethodsforDNAmicroarrays, Bioinformatics 17(6):
520-525. 其中关于KNN的介绍如下:
TheKNN-basedmethodselectsgeneswithexpressionprofilessimilartothegeneofinteresttoimputemissingvalues.IfweconsidergeneAthathasonemissingvalueinexperiment1,thismethodwouldfindKothergenes,whichhaveavaluepresentinexperiment1,withexpressionmostsimilartoAinexperiments2–N(whereNisthetotalnumberofexperiments).Aweightedaverageofvaluesinexperiment1fromtheKclosestgenesisthenusedasanestimateforthemissingvalueingeneA.
Intheweightedaverage,thecontributionofeachgeneisweightedbysimilarityofitsexpressiontothatofgeneA.
下面分析数据
首先要安装 R (http:
//www.r-project.org/)
然后下载安装叫做impute的package
>source("http:
//bioconductor.org/biocLite.R")
>biocLite("impute")
impute是专门用KNN法进行缺失值填充的Rpackage:
设置好当前工作目录(Windows是在R的菜单栏->文件->改变工作目录…设置,Linux下用setwd()函数)
然后在R控制台输入以下代码:
library(impute)
#导入imputepackage
raw<-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=362436069)
#impute.knn()使用一个矩阵作为第一个参数,其他参数这里使用的是默认值
write.table(imputed$data,file='imputed_data.txt')
#write.table()把数据保存在当前工作目录下的文件中,文件名用file=’‘指定,这一步不是必须的
imputeddata<-imputed$data
#imputed$data是在R中储存imputed后的数据的矩阵
现在在R里输入imputed,即填充好的数据矩阵,是不是NA值全都没了?
关于imputepackage的详细Documentation在http:
//bioconductor.fhcrc.org/packages/release/bioc/html/impute.html
全部数据文件:
from:
用R和BioConductor进行基因芯片数据分析(三):
计算median
接前一篇:
我们已经知道要分析的数据对每个基因有3个重复测定值,经过缺失值填充后,每个基因都有3个可用值。
这一步很简单,就是取这3个值的中位数,即median。
方法很多,在excel中可以用median函数;
在R中以下代码进行操作:
get_median<-function(i,j){
num_vec<-c(imputeddata[i*3-2,j],imputeddata[i*3-1,j],imputeddata[i*3,j])
median(num_vec)
}
#Asimplefunctiontocalculatemedianvalueofthreereplicates
dimrow<-(dim(imputeddata)[1])/3
mediandata<-matrix(data=NA,nrow=dimrow,ncol=dim(imputeddata)[2],byrow=TRUE,dimnames=NULL)
#Createablankmatrixtostoremedianvalues
for(iin1:
dimrow){
for(jin1:
dim(imputeddata)[2]){
mediandata[i,j]<-get_median(i,j)
}
}
#Assignmedianvalueusingthefunctionget_median()
现在我们得到了中位数的数据,储存在mediandata对象里,行数是缺失值填充数据imputeddata的1/3,doublecheck一下:
> dim(imputeddata)
[1]1157120
> dim(mediandata)
[1]385720
from:
用R和BioConductor进行基因芯片数据分析(三):
计算median
接前一篇:
我们已经知道要分析的数据对每个基因有3个重复测定值,经过缺失值填充后,每个基因都有3个可用值。
这一步很简单,就是取这3个值的中位数,即median。
方法很多,在excel中可以用median函数;
在R中以下代码进行操作:
get_median<-function(i,j){
num_vec<-c(imputeddata[i*3-2,j],imputeddata[i*3-1,j],imputeddata[i*3,j])
median(num_vec)
}
#Asimplefunctiontocalculatemedianvalueofthreereplicates
dimrow<-(dim(imputeddata)[1])/3
mediandata<-matrix(data=NA,nrow=dimrow,ncol=dim(imputeddata)[2],byrow=TRUE,dimnames=NULL)
#Createablankmatrixtostoremedianvalues
for(iin1:
dimrow){
for(jin1:
dim(imputeddata)[2]){
mediandata[i,j]<-get_median(i,j)
}
}
#Assignmedianvalueusingthefunctionget_median()
现在我们得到了中位数的数据,储存在mediandata对象里,行数是缺失值填充数据imputeddata的1/3,doublecheck一下:
> dim(imputeddata)
[1]1157120
> dim(mediandata)
[1]385720
from:
用R和BioConductor进行基因芯片数据分析(五):
芯片间归一化
接前一篇:
用R和BioConductor进行基因芯片数据分析(四):
芯片内归一化
上次进行了芯片内的归一化,但是我们的数据来自于10张芯片,为了让这10张芯片之间有可比性,需要进行芯片间归一化。
具体原理就不介绍了。
这里用到Bioconductor的一个package,叫做limma,以及其中的函数normalizeBetweenArrays()
由于normalizeBetweenArrays()需要logintensity或logratio作为输入,于是先进行log转化:
#logtransformation
norm_log<-matrix(data=NA,nrow=dim(normed)[1],ncol=dim(normed)[2],byrow=TRUE,dimnames=NULL)
for(iin1:
dim(normed)[1]){
for(jin1:
dim(normed)[2]){
norm_log[i,j]<-log(normed[i,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=c(0,1.35),xlab='logintensity')
for(iin2:
20){
lines(density(norm_log[,i]),type='l')
}
lines(density(norm_log_btw_array[,1]),type='l',col='green')
for(iin2:
20){
lines(density(norm_log_btw_array[,i]),type='l',col='green')
}
text(1.5,c(0.8,1.0),labels=c('BEFOREnormalization','AFTERnormalization'),col=c('black','green'))
用R和BioConductor进行基因芯片数据分析(六):
差异表达基因
接前一篇:
用R和BioConductor进行基因芯片数据分析(五):
芯片间归一化
经过一系列的预处理,包括缺失值填充,中位数计算以及归一化,我们的数据终于可以用啦。
下面我们就来分析一下newpopulation和oldpopulation的个体是否有差异表达基因。
判断一个基因是否差异表达有许多方法,最早使用的就是看logratio的绝对值是否大于2,这种方法早已废弃。
下一个想到的也许是t-test,诚然t-test可以统计地判断一个基因是否差异表达,但是对于有数千数万基因的芯片来说,会有很高的错误发现率(FalseDiscoveryRate,FDR),如果pvalue<0.05,则10000个基因里有500个基因实际没有差异表达而被误认为是差异表达。
因此t-test方法需要改进。
于是Westfall&Young(1993)提出了Step-downmaxTandminPmultipletestingprocedures,大意就是比较几个group间有没有差异基因表达,就通过随机置换这些group的标记,相当于随机互换group的成员,模拟一个空分布(nulldistribution),以此计算调整后的pvalue,这个方法可以极大的减小Family-wiseErrorRate(FWER).
以下分析就使用Step-downmaxTandminPmultipletestingprocedures,需要求助于Bioconductor的multtestpackage的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,B=100000)
默认随机置换次数B=10000,对于microarray来说B应该比10000大很多,所以这里取B=100000
以下是画图:
rawp<-maxttt$rawp[order(maxttt$index)]
plot(sort(rawp),type='p',col=1,ylim=c(-0.05,1.00),ylab='pvalue')
lines(sort(maxttt$adjp),type='p',col='red')
#minadj-p:
sort(maxttt$adjp)[1]0.0163
#rawp:
>sort(rawppp)[170][1]0.0493>sort(rawppp)[171][1]0.0502170个rawp小于0.05
abline(h=0.05,col='blue')
text(1000,c(0.6,0.7),labels=c('rawp-value','adjustedp-value'),col=c('black','red'))
text(1000,0.08,labels='p=0.05',col='blue')
可见调整后只有一个基因的pvalue小于0.05,而未调整的有170个基因的pvalue小于0.05,可以说虽然此方法降低了错误发现率,但是也导致了很高的Falsenegative.
此外可以考虑使用multtestpackage的mt.rawp2adjp()函数,这个函数可以通过”Bonferroni”,“Holm”,“Hochberg”,“SidakSS”,“SidakSD”,“BH”,“BY”等方法调整pvalue,不过对我们的数据来说都过于严格了。
procs<-c("Bonferroni","Holm","Hochberg","SidakSS","SidakSD","BH","BY")
adjps<-mt.rawp2adjp(rawp,procs)
plot(sort(adjps$adjp[,1]),ylab='pvalue')
for(iin2:
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(SignificanceAnalysisofMicroarrays)分析。
有兴趣的请看参考资料。
。
参考资料:
课堂讲义:
Differentialexpression. Identifyingdifferentiallyexpressedgenes—notionsonmultipletestingandp-valueadjustments.
Dudoit,S.,Yang,Y.H.,Speed,T.P.,andCallow,M.J.(2002), StatisticalmethodsforidentifyingdifferentiallyexpressedgenesinreplicatedcDNAmicroarrayexperiments, Statistica