R语言基础实战培训

R语言基础实战培训目录简介函数作图生信及案例解读包（package）学习技巧R简介第一章历史R是“GNUS”，一个能够自由有效地用于统计计算和绘图的语言和环境，它提供了广泛的统计分析和绘图技术，包括线性和非线性模型、统计检验、时间序列、分类、聚类等方法。我们更倾向于认为R是一个环境，在R环境里实现了很多经典的、现代的统计技术。

RobertGentleman

RossIhaka简介丰富的资源

涵盖了多种行业数据分析中几乎所有的方法。良好的扩展性

十分方便得编写函数和程序包，跨平台，可以胜任复杂的数据分析、绘制精美的图形。完备的帮助系统

每个函数都有统一格式的帮助，运行实例。GNU软件

免费、软件本身及程序包的源代码公开。跨平台

可在多种操作系统下运行，如Windows、MacOS、各种版本Linux等。特点用户需要对命令熟悉与代码打交道，需要记住常用命令。对于非计算机专业出身人员具有一定门槛。占用内存所有的数据处理在内存中进行，不适于处理超大规模的数据。运行速度稍慢即时编译，约相当于C语言的1/20。学习曲线陡峭隐藏在后面的统计知识很难在短时间内掌握。特点R使用R安装Rstudio使用简单会话类和对象R安装/R操作界面工具栏控制端口RStudioIDE的核心目标是集中为用户创造一个优雅并高效的使用R的环境。——J.J.AllaireRStudio是一种R语言的集成开发环境(IDE)，其亮点：出色的界面设计及编程辅助工具。跨平台运行，包括windows，Mac，Linux以及网页版。免费RStudio的创始人和CEO/81.9M检索帮助文档查看文件获取扩展环境中的变量历史记录文本编辑命令运行类和对象数值

Numeric如

100,0,-4.335字符

Character如“China”逻辑值Logical如TRUE,FALSE因子型Factor表示不同类别复数型Complex如：2+3i基本元素类和对象向量vector

一系列元素的组合如c(1,2,3);c("a","a","b","b","c")矩阵matrix二维的数据表

是数组的一个特例x<-1:12;dim(x)<-c(3,4)[,1][,2][,3][,4][1,]14710[2,]25811[3,]36912数组array矩阵的自然推广,数据也只能有一种模式。数据框dataframe由于不同的列可以包含不同模式（数值型、字符型等）的数据，数据框的概念较矩阵来说更为普遍。类和对象列表list列表就是一些对象的有序集合。列表允许你整合若干（可能无关的）对象到单个对象名下。简单会话>1+1[1]2>print("Helloworld!")[1]"Helloworld!“>x<-c(1,2,5)>x[1]125>x[3][1]5>mean(x)#求平均值[1]2.666667R是一种解释性语言，输入后可直接给出结果。功能靠函数实现。函数：print(),c(),mean()赋值符：<-，有时可用=,->代替注释符：#下标：从1开始函数名(变量,参数=)数学计算函数是R中最基础的部分，除了以下列举了几个常用的基本计算函数，还有三角函数计算、复数计算、方程计算等等。加减乘除+;-;*；/交集intersect(x,y)绝对值abs(x)并集union(x,y)判断正负sign(x)差集setdiff(x,y)指数函数exp(x)取唯一unique(x)对数函数log(x,base=a)比较计算==,>,<,!=,<=,>=求和sum(x)逻辑计算&,|,&&,||,xor均值mean(x)中位数median(x)最大值max(x)最小值min(x)向量赋值与计算练习R的安装简单数值运算10名婴儿的体重分布；体重与年龄的关系赋值计算练习向量创建向量（创建向量1,2,3,5,6,7，-3）访问向量（访问向量中的第3个元素；第3、第5元素；第2-5元素）矩阵R函数第二章常用函数数学相关计算数学相关计算PART1统计分析

R作为一门为统计计算和作图而生的语言，其统计分析可以称之致胜强项。这部分内容我们介绍常用的描述性统计量、假设检验、回归等的R包和函数。

一、描述性统计我们以R语言中的鸢尾花(iris)数据为例，进行描述性统计分析。iris是数据挖掘常用到的一个数据集，包含150种鸢尾花的信息，每50种取自三个鸢尾花种之一（setosa,versicolour或virginica)。每个花的特征用下面的5种属性描述萼片长度(Sepal.Length)、萼片宽度(Sepal.Width)、花瓣长度(Petal.Length)、花瓣宽度(Petal.Width)、类(Species)。

数据结构如下：描述性统计最常用的基础函数是summary()，如以下结果所示，该函数可提供数值型变量的最小值（Min）、最大值（Max）、四分位数（1stQu.，Median,3rdQu.）和均值（Mean），以及因子向量和逻辑型向量的频数统计。>summary(iris)Sepal.Length Sepal.Width Petal.Length Petal.Width SpeciesMin.:4.300 Min.:2.000 Min.:1.000 Min.:0.100 setosa:501stQu.:5.100 1stQu.:2.800 1stQu.:1.600 1stQu.:0.300 versicolor:50Median:5.800 Median:3.000 Median:4.350 Median:1.300 virginica:50Mean:5.843 Mean:3.057 Mean:3.758 Mean:1.199 3rdQu.:6.400 3rdQu.:3.300 3rdQu.:5.100 3rdQu.:1.800 Max.:7.900 Max.:4.400 Max.:6.900 Max.:2.500 二、假设检验分析

假设检验(HypothesisTesting)是数理统计学中根据一定假设条件由样本推断总体的一种方法。首先要对总体分布参数设定一个假设（零假设H0），然后从总体分布中抽样，通过样本计算所得的统计量来对总体参数进行推断。假定零假设为真，如果计算获得观测样本的统计量的概率非常小，小于预先给定的显著性水平，便可以拒绝原假设，接受它的对立面（研究假设H1）。

常用的假设检验方法有u检验法、t检验法、χ2检验法(卡方检验)、F检验法，秩和检验等。1）t检验法

适合于变量为连续型的两组间比较，且并假设两组数据均呈正态分布，可用t.test()函数进行分析，格式为t.test(x~y,data)（x为连续型变量，y为二分类变量，data为包含了变量的数据框），或者t.test(x,y)（x,y分别为两组的连续型变量）

。t.test(x,y,var.equal=TRUE)#T检验（双尾，非配对，等方差）t.test(x=(1:10),y=c(7:20))结论：由于p<0.05，拒绝H0，接受H1，即setosa和versicolour的萼片长度Sepal.Length不同。

以上为两独立样本t检验的例子，配对样本则添加参数paired=TRUE，此外还有参数var.equal=TRUE以假定方差相等，参数alternative="less"或alternative="greater"来进行有方向的检验。

更多参数的详细介绍请参考帮助文档。

函数作用的对象是什么样的？有哪些参数？如何改动？……help(t.test)?t.testDescription

#函数描述Usage

#默认选项Arguments

#参数Details

#详情Author(s)

#作者References

#参考文献Examples

#举例帮助文件中的内容：2）方差分析当变量多于两组时，t检验不再适用，则单因素方差分析，函数为aov()。如下，检验三种Species的iris萼片长度Sepal.Length是否相同：>fit<-aov(Sepal.Length~Species,iris)>summary(fit)结论：三种Species的iris萼片长度Sepal.Length存在显著差异（p<0.05）

3)卡方检验卡方检验常常用来评价分类资料中类别型变量的关系，包括：两个率或两个构成比比较的卡方检验；多个率或多个构成比比较的卡方检验以及分类资料的相关分析等。典型的H0假设是变量之间独立，H1假设是不独立。如：A、B组病人，对某药物的反应如下表所示，比较两组病人对药物的反应是否相同。>data<-matrix(c(180,90,133,200),nrow=2)>chisq.test(data,correct=TRUE)Pearson'sChi-squaredtestwithYates'continuitycorrectiondata:dataX-squared=41.601,df=1,p-value=1.119e-10结论：拒绝H0，接受H1，两组病人对药物的反应不同。敏感不敏感A组180133B组902004)Fisher精确检验Fisher精确检验适用于分析2x2列联表，并检验行变量和列变量是否独立。它基于超几何分布，作用于离散变量。相比于卡方检验，p值对于所有样本数量都是准确的，而不是近似的结果。

如：检验肺癌与吸烟是否有关的数据：>data<-matrix(c(9,4,6,12),nrow=2)>fisher.test(data)Fisher'sExactTestforCountDatadata:datap-value=0.07317alternativehypothesis:trueoddsratioisnotequalto195percentconfidenceinterval:0.784656327.9487259sampleestimates:oddsratio4.268185结论：不拒绝H0，该数据不能说明肺癌与吸烟有关。有无吸烟96不吸烟412三、相关性分析相关性分析是指对两个或多个具备相关性的变量元素进行分析，从而衡量两个变量因素的相关密切程度。相关系数可以用来描述定量变量之间的关系。相关系数的符号（+/-）表明关系的方向（正相关或负相关），其值的大小表示关系的强弱程度（完全不相关时为0，完全相关时为1）。R可计算多种相关系数，如最基础函数cor()可计算的Pearson相关系数、Spearman相关系数、kendall相关系数。>cor(x,y,method="")#其中x,y为拟进行相关分析的变量，method用来指定相关的类型相关性检验

相关系数进行显著性检验，常用的H0假设为变量间不相关（即总体的相关系数为0）。

>cor.test(x,y,method="")Pearson'sproduct-momentcorrelationdata:xandyt=16.763,df=18,p-value=1.98e-12alternativehypothesis:truecorrelationisnotequalto095percentconfidenceinterval:0.92277910.9880752sampleestimates:cor0.969433四、回归分析回归分析是确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法，可以由给出的自变量估计因变量的条件期望。回归分析按照涉及的变量的多少，分为一元回归和多元回归分析；按照自变量和因变量之间的关系类型，可分为线性回归分析和非线性回归分析。1）简单线性回归我们依然利用R语言中内置的iris数据，希望以花瓣长度(Petal.Length)为因变量，以花瓣宽度(Petal.Width)为自变量进行线性回归，得到以花瓣宽度预测花瓣长度的公式。#拟合>fit<-lm(Petal.Length~Petal.Width,data=iris)>summary(fit)#作图>plot(iris$Petal.Width,iris$Petal.Length)>abline(fit)Call:lm(formula=Petal.Length~Petal.Width,data=iris)Residuals:Min1QMedian3QMax-1.33542-0.30347-0.029550.257761.39453Coefficients:EstimateStd.ErrortvaluePr(>|t|)(Intercept)1.083560.0729714.85<2e-16***Petal.Width2.229940.0514043.39<2e-16***---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1Residualstandarderror:0.4782on148degreesoffreedomMultipleR-squared:0.9271,AdjustedR-squared:0.9266F-statistic:1882on1and148DF,p-value:<2.2e-16结论：Petal.Length=1.08+2.23×Petal.Width回归系数（2.03）显著不为0（p<0.05），花瓣宽度每增加1，花瓣长度将预期增加2.23；此模型可解释花瓣长度92.7%的方差。

2）多元线性回归以萼片长度（Sepal.Length）和

萼片宽度（Sepal.Width）为自变量，花瓣长度(Petal.Length)为因变量，进行线性回归。参考代码如下，结果解释与上一部分相同。>fit<-lm(Petal.Length~Sepal.Length+Sepal.Width,data=iris)>summary(fit)3）logistic回归

logistic回归又称logistic回归分析，是一种广义的线性回归分析模型。可通过一系列连续型和/或类别型预测变量来预测二分类结果变量。

结论：除了outcome2之外，其他变量的贡献都不显著。outcome2回归系数为-0.4543。QQ图卡方检验chisq.test()Fisher精确检验fisher.test(

)方差分析aov()回归分析lm()逐步回归step()相关性分析cor()qqnorm()qqline()生存分析Surv()survdiff()coxph()K-M生存曲线统计学中的常用函数路径选择修改路径：setwd()、或手动选择目前路径：getwd()不管是读取数据还是写入，R都是在工作路径中完成的。所以首先我们要知道我们的R所在的工作路径是在哪里。使用getwd()函数来获取我们的工作路径。查看目前路径下文件：dir()下面查看工作路径里面有哪些文件，使用dir()函数文本输入/输出文本读取：read.table()文本输出：write.table()Pdf输出：pdf(file)data<-read.table(file=“test.txt",sep="\t",header=TRUE,s=1)write.table(dif,"dif.txt",s=TRUE,sep="\t")练习进入路径，并读入文件test1.txt切换工作路径学习几个函数的example输出文件（txt文件，pdf文件）到制定路径共表达分析、fisher检验分析作图第三章常用绘图函数两大绘图系统：传统绘图系统直观地实时反映绘图和分析数据的逻辑两步=图+修饰/添加=执行一系列函数grid绘图系统（ggplot2和latttice）绘图=使用一次函数调用（一次成图）

R具备卓越的绘图功能，通过参数设置对图形进行精确控制。可以输出Jpg、tiff、eps、emf、pdf、png等各种格式。

Sourse:https:///x<-runif(50,0,2)#随机产生y<-runif(50,0,2)plot(x,y)abline(h=.6,v=.6)text(0.6,0.6,"Hi!",col="red")传统绘图系统图+修饰/添加分步绘图此处记得先下载R包，否则会报错箱式图boxplot（）test<-read.table(“expr.matrix.txt”)pdf("boxplot.pdf",width=4,height=4)par(mar=c(2+round(max(nchar(sampleNames(RawData)))/2),4,2,1))title<-paste("GSExxxxx","afternormalization",sep='')boxplot(test,boxwex=0.6,notch=T,main=title,outline=FALSE,las=3)dev.off()生信魅力第四章帮助预测疾病相关的潜在基因，以及该基因潜在的作用靶点、上游调控转录因子等，从而指导实验方向、缩小试验范围、简化试验流程。为基金申请提供支持，通过强大的信息数据的收集整理，减少投入增强研究目的性；且通过整合技术优势，指导提高临床诊断水平。丰富的免费的数据资源和生物软件工具：数据库：NCBI,EBI,DDBJ,PDB,SWISS-PROT等；软件工具：Cytoscape，BLAST,R语言，String,David等生物信息学魅力59一篇好的生信文章应该包含先进的生信分析思路（流程）较好的原始数据漂亮的图片、表格创新的结论原始数据来源基于客户样本的高通量数据测序数据转录组、基因组、表观组、蛋白组…芯片数据基于公共数据库的二次挖掘NCBIGEO、TCGA、EBI等基于临床数据（客户/公共）医学统计学分析生存曲线分析样本量大质量好分组好芯片简介生物芯片也称基因芯片(GeneChip)、生物芯片(Biochip)或微阵列(Microarray)。DNA芯片原理DNA芯片的基本原理与生物学中Southern杂交等实验技术相似，都是利用DNA双螺旋序列的互补性，即两条寡聚核苷酸链以碱基之间形成氢键配对(A与T配对，形成两个氢键；G与C配对，形成三个氢键)。DNA芯片通常以尼龙膜、玻璃、塑料、硅片等为基质材料，固着特定序列DNA单链探针(Oligo)，并与被检测序列单链cDNA序列互补结合(通常称杂交)。被检测序列用生物素或荧光染料标记，通过荧光染料信号强度，可推算每个探针对应的样品量。一张DNA芯片，可固着成千上万个探针，具体数目则取决于芯片设计和制备方法。63芯片类型用于基因组研究的SNP和CNV芯片用于mRNA表达研究的基因表达谱芯片用于转录调控研究的microRNA芯片和LncRNA芯片以及用于表观遗传研究的DNA甲基化芯片64lncRNA芯片SNP、CNV芯片microRNA芯片mRNA芯片甲基化芯片主流的芯片制造商affymetrix公司美国agilent公司美国illumina公司65二代测序简介测序（Sequencing）–使用专业仪器精确测定DNA或RNA序列碱基顺序的过程SangersequencingIn2006,NextGenerationSequencing(NGS)——UpdateofSangersequencingNGS的三大特点并行测序“massivelyparallelsequencing”Hundredsofthousandsto

hundredsofmillionsofreactionsimaged

perinstrumentrun测序耗费少Sharplydecreasingcost序列读长短Shorterreadlengthsthancapillarysequencers66测序技术发展历程67不同厂商提供的NGS技术特征目前，Roche454测序仪业务已经关闭68NGS技术vs

芯片技术芯片技术的内在缺陷闭合杂交系统Themicroarraycouldbedescribedasa“closedsystem”becauseinformationaboutRNAsorDNAsislimitedbythetargetsavailableforhybridization.探针杂交信号饱和问题Maximumexpressionvalueisknowninadvancelackfidelity探针杂交信号检测灵敏度下限问题69测序芯片检测未知基因可以检测未知基因只能检测已知基因分析内容分析内容多（可变剪切，融合基因）可分析内容少检测范围比较能反应真实的表达情况（最小1个read，最大10E+6）限于检测仪的荧光强度和灵敏极限（10E+4）案例解读第五章转录组数据解读Html版本Pdf版本转录组数据解读原始数据转录组数据解读原始数据一般是fq格式文件右图是1个样本：双端测序，所以有两个fq文件质控：FastQC软件测序数据单碱基质量分布图Basecontent分布图GCcontentSequencebasequality分布图转录组数据解读过滤测序下机数据包含一些带接头、低质量的reads，这些序列会对后续的信息分析造成很大的干扰，为了保证后续信息分析质量，需要对下机数据进行进一步过滤1)接头污染去除，采用cutadapt（version1.2.1）(MartinM，2011）去除3'端的接头污染，至少10bpoverlap（AGATCGGAAG），允许20%的碱基错误率；2)质量过滤，采用5bp滑动窗口法，窗口的平均质量分数>=Q20，允许3’截短，最短不低于50bp，序列中无不确定碱基（”N”）。转录组数据解读比对将测序获得的序列与基因组进行比对，将其还原到基因组上，再根据已有的基因组注释，可得知是何基因表达产物。对于转录组来说，还需要考虑到测序的实际序列为经过转录剪切后序列，需要进一步识别可变剪切位点。专用于转录组比对和分析的软件有基于Bowtie/Bowtie2的Tophat/Tophat2等转录组数据解读定量统计比对到每一个基因上readcount值，作为基因的原始表达量，获得样本的Reads值矩阵。分析软件：HTSeq0.6.1p2转录组数据解读横行代表基因，可能有几万行；列代表样本。差异分析基因表达差异分析用于检测不同样品之间基因表达模式的变化，通过基因表达量差异来研究基因的功能基因表达差异分析分为两样品差异分析和多样品差异分析。两样品差异分析常见于处理/对照、野生/突变等类型的比较分析。实验设计简单，结果可靠。多样品比较差异表达基因，包括单因素多水平处理：浓度梯度处理样本、时间梯度处理等。差异阀值（经验值）表达量倍数差异阈值|logFC|>1，显著性阈值P-value<0.05转录组数据解读差异分析聚类分析：聚类分析用于判断不同样本基因（或基因模块）表达模式，并将相同或相近表达模式的基因按照距离进行排序分析软件：heatmap.2（hclust）转录组数据解读功能富集分析GO（GeneOntology）富集分析分子功能(molecularfunction)细胞位置(cellularcomponent)生物学过程(biologicalprocess)KEGGpathway富集转录组数据解读蛋白互作网络分析（PPI）或调控网络（eg：miRNA-mRNA）转录组数据解读案例解读美国科学院报（ProceedingsoftheNationalAcademyofSciences

of

theUnitedStatesofAmerica，简称PANS）曾于2011年发表非小细胞肺癌成纤维细胞的基因芯片研究。PANS是世界公认的四大名刊之一（Cell，Nature，Science，PNAS），2014年的影响因子为9.674，该期刊发表的文章可代表某方向研究的一流水平。分析流程第一步：原始数据及预处理文章从15个非小细胞肺癌患者中提取出肺癌成纤维细胞和正常成纤维细胞样本，并做了Affymetrix平台下的芯片数据。这个芯片数据如何产生就不在这里描述了。以affymetrix芯片为例：读取数据

从NCBIGEO（或Arrayexpress）下载原始数据可整个数据包下载，比如GSE22863_RAW.tar，解压到GSE22863_RAW文件夹(放在工作目录里)，会发现里面包含一系列.cel.gz压缩包，不用再解压了。选取需要的GSM文件(normal和cancer)，和研究无关的删除。或者下载时直接下载需要的样本文件，点击download里面的custom，勾选要下载的文件。通过R中的Affy包进行预处理。案例第二步：差异基因筛选得到处理过的芯片数据后，文章分析肺癌成纤维细胞和正常成纤维细胞的基因表达水平差异，筛选出肺癌成纤维细胞中高表达或者低表达基因（通常称之为差异表达基因,DEGs）。Limma是常被用于计算基因表达水平变化值。右图展示输出结果，从表格可以很明显看出哪些基因在癌症样本中的表达水平发生变化。案例第三步：差异基因热度图文章展示的差异基因分析热度图Differentiallyexpressedgenesin15pairedCAFvs.NFandincorrespondingNSCLCstromavs.normallung.Atotalof46differentiallyexpressedgenes[22up-regulated(bluesidebar)and24down-regulated(yellowsidebar)]wereidentifiedbyusingpairedSAManalysis(absolutefoldchange>2;q<0.1).案例第四步：富集分析差异表达基因GO功能和KEGG通路GOanalysisstillshowedthatthesegeneswereenrichedingenesoftheextracellularregion(9/14,P=0.006)andKEGGpathways,includingECM–receptorinteraction(3/14,P=0.005)andfocaladhesion(3/14,P=0.01),在得到高表达和低表达基因后，原文作者随后对这些基因作了Geneontology和KEGG通路富集分析。云生信提供了齐全的基因富集分析工具，见模块“富集分析工具”。案例第五步：生存曲线分析作者根据收集非小细胞肺癌病人的临床数据做出了生存曲线图，相信广大医生对此很感兴趣。云生信也可以对临床数据做生存曲线分析，那就是模块“作图工具”中的“生存曲线分析”小工具。案例第6步：蛋白质网络分析作者基于差异表达基因做了蛋白质-蛋白质相互作用网络分析。案例实验验证免疫组化染色RT-qPCR第六章芯片数据分析中的常用R包R程序包是什么？R程序包是多个函数的集合，具有详细的说明和示例。Window下的R程序包是经过编译的zip包。每个程序包包含R函数、数据、帮助文件、描述文件等。为什么要安装R包？R程序包是R功能扩展，特定的分析功能，需要用相应的程序包实现。例如：系统发育分析，常用到ape程序包，群落生态学vegan包等。R

packages站在巨人的肩膀上七个核心包R包的类型二进制代码包（Binarypackage）即得即用型，但是依赖平台，比如Windows和Linux平台下不同。源代码包（Sourcepackage）

此类包可以跨平台使用，但用之前需要处理或者编译。同时，源代码包可以查看到程序源代码，便于查找、修改和引用。数学计算，统计计算，日期函数，包加载函数，数据处理函数，函数操作函数，图形设备函数R

packages9694个R

package部分展示以ggplot2包为例：一、R包的安装1、用函数install.packages()：如果已经连接到互联网，在括号中输入要安装的程序包名称，选择镜像后，程序将自动下载并安装程序包。如：install.packages(“ggplot2")2、Bioconductor平台R包的安装：##tryhttp://ifhttps://URLsarenotsupportedsource("/biocLite.R")biocLite("illuminaHumanv2.db")

3

安装本地zip包Packages>installpackagesfromlocalfiles>Browse>InstallRstudioR

packages如何使用：library(ggplot2)二、芯片预处理（affy包）

affy包可以读取Affymetix公司的基因表达芯片CEL格式数据，处理成表达矩阵，并进行预处理，生成可以直接进行后续差异分析的文件。

流程主要包括以下几步：以affy为例（附件额外代码提供Illumina、agilent等测序平台芯片预处理方法代码）##加载所需R包，需安装library(affy)##将CEL文件所在的文件夹设置为默认路径RawData<-ReadAffy()##背景校正和标准化RMAdata<-rma(RawData)##输出表达矩阵，行为探针，列为样本rma.data<-exprs(RMAdata)##Exportingdata

write.table(rma.data,file="norm.txt")二、芯片预处理（affy包）

#####对探针进行注释######。library（annotate）probesets<-s(rma.data)#获得探针的名称library(paste(annotation(RMAdata),".db",sep=""),character.only=TRUE)#确定注释包，首次使用需要下载annotation(RMAdata)#确认该芯片的注释包名称source("/biocLite.R")biocLite("hgu133plus2.db")#如果该包已经安装了，那么就不用再次安装。pro.symbols<-getSYMBOL(probesets,annotation(RMAdata))symbols.matrix<-data.frame(Symbols=pro.symbols,rma.data)expr.matrix<-aggregate(symbols.matrix[,-1],by=list(symbols.matrix[,1]),mean)#多个探针对应到同一个基因，取均值s(expr.matrix)<-expr.matrix[,1]#行名赋值为symbolexpr.matrix<-expr.matrix[,-1]#删除第一列write.table(expr.matrix,file='uniquenorm.txt',sep='\t',quote=F,s=T)三、差异分析（limma包）

limma包的功能十分强大，既可以读取、预处理芯片数据，同时也有差异化基因分析的算法（limma:LinearModelsforMicroarrayData），支持单通道、双通道芯片等所有芯片，甚至定量PCR和RNA-Seq，并且每一步处理都提供了多种可选方法，在芯片数据处理流程中占有举足轻重的地位。

我们演示以limma包经验贝叶斯方法进行芯片数据差异分析的流程：library(limma)#加载limma包#产生一个design矩阵，用于描述样品分组data<-expr.matrix#数据读入group<-c(rep(1,3),rep(2,3))#分组design<-model.matrix(~-1+factor(group))#设计矩阵colnames(design)<-c("case","con")#命名#差异分析信息contrast.matrix<-makeContrasts(case-con,levels=design)#差异分析的关键，case，con是前面命名的组名。#建立线性模型，并进行分组差异分析和p值计算、校正fit<-lmFit(data,design)fit1<-contrasts.fit(fit,contrast.matrix)fit2<-eBayes(fit1)#结果文件预览及输出results<-decideTests(fit2,method="global",adjust.method="none",p.value=0.05,lfc=0.585)summary(results)dif<-topTable(fit2,coef=1,n=nrow(fit2),lfc=0)write.table(dif,"allgene.dif.txt",s=TRUE,sep="\t")#所有基因的差异表达情况输出。#limma也支持多组的两两比较只需修改参数：colnames(design)<-c(“case”,“case2”，"control")result2<-topTable(fit2,coef=2,number=nrow(data),adjust.method="BH",sort.by="B",resort.by="M")limma包的分析结果展示:Log2(FoldChange)平均表达值t检验效应量显著性校正后pvalue差异表达的对数比四、主成分分析主成分分析（PCA）是一种数据降维技巧，它能将大量相关变量转化为一组很少的不相关变量，这些无关变量称为主成分。在基因芯片分析中，可获得样本在实验组和对照组之间的直观分布情况，从而便于对异常样本进行检测和去除，否则异常样本的存在将会对差异基因的鉴定等后续分析造成不利影响。install.packages('ggfortify')library(ggfortify)dataPCA<-t(rma.data)#转置行为样本、列为gene（没有rma.data的读取课件norm.txt）dataPCA<-as.data.frame(dataPCA)#转为表格dataPCA$group<-c(rep("Case",3),rep("Ctrl",3))#dataPCA$group<factor(dataPCA$group,levels=c("Case","Ctrl"))pdf("PCA.pdf",height=6,width=6)autoplot(prcomp(dataPCA[,-ncol(dataPCA)]),data=dataPCA,colour='group')dev.off()

从图中看出，两组样本被很好地分开，若有其他颜色的点混入某个组分，则后续分析可考虑剔除。五、聚类热图（pheatmap包）聚类热图使用颜色直观呈现多样本多个基因的全局表达量变化，可呈现多样本或多基因表达量的聚类关系（样本聚类和基因聚类）。

常用函数有gplots包：heatmap.2()、ggplot2包:ggplot()和pheatmap包：pheatmap()。#一般步骤##导入数据（从外部文件导入，一般是txt，csv格式，行为基因，列为样本）：>data<-read.csv(“F:/filepath/heatmap.csv”,s=1)

#你自己的工作目录>pheatmap(data)#先尝试直接生成热图，再逐步调整*data，你自己的表达谱数据>pheatmap(data,scale="row",#根据row标化数据color=colorRampPalette(c("navy","white","firebrick3"))(50),#更改颜色cutree_rows=2)

#将热图分为高低表达两部分##附件代码中

pheatmap.r对刚刚做完预处理及差异分析得到的差异基因提取表达值作聚类热图。各类参数设置自己可以尝试修改，看看每次修改，热图的哪些部分会产生变化。六、GO/pathway富集分析（clusterProfiler包）GO是Geneontology的缩写，分别从功能、参与的生物途径及细胞中的定位对基因产物进行了标准化描述，通过GO富集分析可以粗略了解差异基因富集在哪些生物学功能、途径或者细胞定位。Pathway指代谢通路，对差异基因进行pathway分析，可以了解实验条件下显著改变的代谢通路，在机制研究中显得尤为重要。clusterProfiler包可以对基因集合进行GO和pathway富集的统计分析和可视化。目前clusterProfiler包在最新版本的R中容易报错，如果报错，调整所使用的R版本。1.GO富集分析和可视化>library(clusterProfiler)>dif_enrich<-dif[which(abs(dif$logFC)>0.585&dif$P.Value<0.05),]#提取差异基因列表gene<-rownames(dif_enrich)#转为向量##将genesymbol转换为Entrez

ID>gene.id<-bitr(gene,fromType=“SYMBOL”,toType=“ENTREZID”,OrgDb=“org.Hs.eg.db”)

#gene你的gene

symbol列表>egoBP<-enrichGO(gene=gene.id$ENTREZID,#需要富集的基因集OrgDb='org.Hs.eg.db',#种属ont="BP",#BP、CC或者MFpvalueCutoff=0.05,#筛选阈值readable=TRUE)>summary(egoBP)#结果展示

结果可视化有多种类型，条形图和气泡图较为常用：>barplot(egoBP,showCategory=8)2.pathway富集分析和可视化Pathway富集分析和可视化与GO富集类似：>egoPathway<-enrichKEGG(gene=gene.id$ENTREZID,organism="hsa",pvalueCutoff=0.01)>dotplot(egoBP)七、Venn图（venndiagram包）用圆圈显示元素集合重叠区域的图示，用来表示两组（或多组）数据之间的关系。常用函数有gplots包：venn()、Venneuler包:venneuler()和VennDiagram包:venn.diagram()数据：a,b,c三个向量，分别是肿瘤分期stageI~stageIII与正常对照组比较的差异基因集合library(VennDiagram)#加载>data<-read.table("venn.txt",