R语言在基因芯片数据处理中的应用要点.docx
《R语言在基因芯片数据处理中的应用要点.docx》由会员分享，可在线阅读，更多相关《R语言在基因芯片数据处理中的应用要点.docx（14页珍藏版）》请在冰豆网上搜索。
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
升级会员 