R语言在基因芯片数据处理中的应用要点Word格式文档下载.docx
- 文档编号:19531868
- 上传时间:2023-01-07
- 格式:DOCX
- 页数:14
- 大小:262.59KB
R语言在基因芯片数据处理中的应用要点Word格式文档下载.docx
《R语言在基因芯片数据处理中的应用要点Word格式文档下载.docx》由会员分享,可在线阅读,更多相关《R语言在基因芯片数据处理中的应用要点Word格式文档下载.docx(14页珍藏版)》请在冰豆网上搜索。
-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="
symkey=FALSE,revC=FALSE,margins=c(10,10),denscol=tracecol,distfun=dist,hclustfun=hclust,dendrogram="
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("
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进行基因芯片数据分析(五):
芯片间归一化
用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)
dim(normed)[1]){
dim(normed)[2]){
norm_log[i,j]<
-log(normed[i,j])/log
(2)
}
然后利用函数进行芯片间归一化:
source("
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='
col='
green'
lines(density(norm_log_btw_array[,i]),type='
text(1.5,c(0.8,1.0),labels=c('
BEFOREnormalization'
'
AFTERnormalization'
),col=c('
black'
用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)查看。
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='
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'
))
text(1000,0.08,labels='
p=0.05'
可见调整后只有一个基因的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='
8){
points(sort(adjps$adjp[,i]),col=i)
因此可以考虑不这么严格的SAM(SignificanceAnalysisofMicroarrays)分析。
有兴趣的请看参考资料。
。
参考资料:
课堂讲义:
Differentialexpression.
Identifyingdifferentiallyexpressedgenes—notionsonmultipletestingandp-valueadjustments.
Dudoit,S.,Yang,Y.H.,Speed,T.P.,andCallow,M.J.(2002),
StatisticalmethodsforidentifyingdifferentiallyexpressedgenesinreplicatedcDNAmicroarrayexperiments,
Statistica
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 语言 基因芯片 数据处理 中的 应用 要点
![提示](https://static.bdocx.com/images/bang_tan.gif)