【R语言高级详细教程】表达谱芯片的聚类分析.doc

【在国际上，R软件的应用是数据分析的主流发展趋势之一，但我发现在国内R软件的使用远不如SPSS、SAS等软件那么流行。为推广R软件的使用，本博客将陆续推出“R高级教程”系列专辑，希望对生命科学领域的科技工作者有少许帮助.】 通常来讲，对于一般的统计分析，基于傻瓜式操作的SPSS（PASW）软件已经足够，但在涉及个性化要求很高的复杂数据处理时，SPSS就开始显得力不从心，这时必须依赖功能更为强大的SAS等软件。以前在自己的科研过程中分析数据多用SPSS、SAS等。在统计遗传和基因组学领域，SAS可以处理很多问题，但与此同时，SAS实现复杂问题过于麻烦，很多问题SAS也不是首选。后来开始运用R环境中的各种免费统计包，特别是Bioconductor的系列分析包，我发觉非常适合生命科学领域的研究者。R有很多优点：（1）免费，不需要去寻找破解版，不用担心版权问题，使用非常方便； （2）功能非常强大，单个包的功能比较有限，但多个包组合起来使用则功能无比强大，远胜于SPSS、SAS等； （3）源代码开放，稍作修改后就能满足个性化的复杂统计分析，满足个性化需求是R的最大特点之一； （4）程序阅读容易，再加上参考学习资料很多，上手比较容易，提高也不是很难，根据个人经验，要比SAS高级阶段的进阶容易许多； （5）国际同行高度认同R，我发现很多专用软件都开发了软件的R版，今后R将是数据分析的主流发展方向。 R软件的安装、基本使用等初级教程就不谈了，随便在官方网站找个学习资料就搞定了。“R系列”专辑拟推出中级、高级分析教程。今天推出基因表达谱芯片的聚类分析专题。 本专题示例芯片数据来自GEO数据库中检索号为GSE11787的Affymetrix芯片的CEL文件，共6个CEL文件，3个正常对照组，3个HPS刺激组，为免疫器官脾脏的表达数据。（一）原始数据的读入、RNA降解评估和标准化pd rawAffyData summary(exprs(rawAffyData) deg plotAffyRNAdeg(deg, col=c(1,2,3,4,5,6) eset summary(exprs(eset) op cols boxplot(rawAffyData,col=cols,names=1:6, main = unnormalized.data)boxplot(data.frame(exprs(eset) ,names=1:6, main = normalization.data, col=blue, border=brown)par(op)（二）聚类分析 原始数据读入，经AffyBatch目标转成ExpressionSet目标后，为提高后续分析（如差异表达基因的检测）的统计功效，往往需要进一步经过Detection Call Filter和IQR filter等过滤（“基因芯片数据的特异性过滤与非特异性过滤”将在另一专题里专门讨论）。 需要说明的是，常规做法是先筛选出差异表达基因，然后只用差异表达基因进行聚类分析（本示例直接用了过滤后的数据集，聚类图的效果稍差一点）。（1）样本聚类dd diag(dd) dd.row row.ord library(latticeExtra)legend lp plot(lp)（2）二维聚类source(/tgirke/Documents/R_BioCond/My_R_Scripts/my.colorFct.R)mydatamydatascale hr hc heatmap.2(mydata, Rowv=as.dendrogram(hr), Colv=as.dendrogram(hc), col=redgreen(75), scale=row, ColSideColors=heat.colors(length(hc$labels), RowSideColors=heat.colors(length(hr$labels), trace=none, key=T) 上述聚类图一般和论文里的聚类图有点不同，聚类的模式不太直观，你也可以用下面的语句进行更直观的作图：mycl mycolhc - sample(rainbow(256);mycolhc myc2 - cutree(hc, h=max(hc$height)/1.5); mycolhr - sample(rainbow(256); mycolhr heatmap(mydatascale, Rowv=as.dendrogram(hr), Colv=as.dendrogram(hc), col=my.colorFct(), scale=row, ColSideColors=mycolhr, RowSideColors=mycolhc)（3）MantelCorrs聚类程序kmeans.result x=exprs(eset2)DistMatrices.result MantelCorrs.result permuted.pval ClusterLists ClusterGenes h=hclust(dist(MantelCorrs.result)plot(h)【注：除Bioconductor图标外，所有图片均为软件实际运行所得】应学生及个别博友的要求，尽管专业博文点击率和反应均很差，但在去San Diego参加PAG会议之前，还是抽时间给出【R高级教程】的第二专题。专题一给出了聚类分析的示例，本专题主要谈在表达谱芯片分析中如何利用Bioconductor鉴定差异表达基因。 鉴定差异表达基因是表达谱芯片分析pipeline中必须的分析步骤。差异表达基因分析是根据表型协变量（分类变量）鉴定组间差异表达，它属于监督性分类的一种。在鉴定差异表达基因以前，一般需要对表达值实施非特异性过滤（在机器学习框架下属于非监督性分类），因为适当的非特异性过滤可以提高差异表达基因的检出率、甚至是功效。R分析差异表达基因的library有很多，但目前运用最广泛的Bioconductor包是limma。 本专题示例依然来自GEO数据库中检索号为GSE11787 的Affymetrix芯片的数据，数据介绍参阅专题一。library(limma)design - model.matrix( -1+factor(c(1,1,1, 2,2,2) 这个是根据芯片试验设计，对表型协变量的水平进行design，比如本例中共有6张芯片，前3张为control对照组，后3张芯片为实验处理组，用1表示对照组，用2表示处理组。其他试验设计同理，比如2*2的因子设计试验，如果每个水平技术重复3次，那么可以表示为：design colnames(design) fit contrast.matrix fit fit2 fit2 resultssummary(results)vennCounts(results)vennDiagram(results)比较遗憾的是，目前limma自带的venn作图函数不能做超过3维的高维venn图，只能画出3个圆圈的venn图，即只能同时对三个coef进行venn作图。上面的venn图只有一个coef，太简单了。下面是一个由本实验室芯片数据得出的三个coef的venn图例：heatDiagram(results,fit2$coef) 红色为control中（与LPS相比）的高表达基因，绿色为control中（与LPS相比）的低表达基因，x轴的数字表示差异表达基因在eset2中所处的位置。xwrite.table(x, file=limma.xls, s=F, sep=t) 将结果写入limma.xls文件中，内容包括AveExpr值（比较组间绝对值的平均差异值）、logFC值（差异倍数）、t值、P值、q值（即adj.P.Val值）和B值。一般logFC值、P值、q值和AveExpr值用来作为差异表达的判断标准，比如差异倍数在2倍以上、绝对差异表达在10以上、P值小于0.01等。在Excel文件中，根据各项判断标准排序，可以很容易地得到差异表达基因列表，这个列表可以用来进行后续的分析，如GO注释、基因网络绘制等。 专题一中提到实际研究中，一般只用差异表达基因进行聚类分析，在R中，根据差异表达结果过滤表达值很简单（具体的值可以依据芯片数据的实际情况设定，比如P值可以设宽松点0.05、logFC的绝对值