二讲基因表达数据分析

第二讲基因表达数据分析

PartI基因测定平台及数据库介绍(GEO、TCGA、SEER数据库的下载）PartII基因芯片数据的预处理分析PartIII基因功能分析及DAVID网络工具应用PartIV综合案例：肿瘤基因表达数据分析PartI基因表达测定平台及数据库介绍近20年来三种不同高通量基因表达测定技术的应用趋势基因表达测定平台与数据库

1.cDNA芯片2.Affymetrix芯片

3.RNA-Seq技术对没有已知参考基因组信息的非模式生物，也可测定转录信息；RNA-Seq技术可以测定转录边界的精度达到一个碱基；RNA-Seq可以用来研究复杂的转录关系；RNA-Seq可以同时测定序列的变异；RNA-Seq背景信号很小，测定的动态范围很大。常用基因表达数据库名称数据库内容GeneExpressionOmnibus（GEO）目前最常用的基因表达数据（NCBI）ExpressionAtlas欧洲生物信息学中心的基因表达数据库SMDStanford基因表达数据库RNA-SeqAtlas正常组织的基因表达谱数据GEPdb基因型、表型和基因表达关系GXD老鼠发育基因表达信息EMAGE老鼠胚胎的时空表达信息AGEMAP老鼠老化的基因表达数据常用基因表达数据库

下载膀胱癌数据TCGA数据库的下载推荐UCSCXena下载TCGA数据

下载膀胱癌数据用RPKM值矩阵数据注释数据基因表达矩阵数据基因注释数据表型数据总生存时间无复发生存时间SEER数据库TheSurveillance,Epidemiology,andEndResultsSEER数据库（）美国国立癌症研究所“监测、流行病学和结果数据库”，是北美最具代表性的大型肿瘤登记注册数据库之一，收集了大量循证医学的相关数据，为临床医师的循证实践及临床医学研究提供了系统的证据支持和资料。SEER数据库及软件下载SEER数据库不像其他数据库，SEER需要注册，而且需要签份保证书，得到官方认可才可以得到一个用户名和账号，用于数据和软件的下载安装。

注册入口：填写信息提交后将迅速收到邮件，进行确认签署协议，发邮件后等半天~1天就会收到SEER发过来的用户名和密码。SEER数据库SEERStat软件双击SEERstat图标链接下载的数据库data是选择数据来源的范围，SEER数据库最近更新到2015年的随访数据，所以尽量选最新的数据集。selection是选择病例筛选的条件。点击edit，然后会有很多选项，选择你需要的病例范围、肿瘤类型等等。SEER数据库有许多自己定义的简写和缩写，比如metsatdx代表远处转移，转移部位和病理类型都用代号表示。建议官网下载“CScodinginstructions”的文件，里面详细列出了该数据库使用的每个coding的解释table是选择需要的临床信息。左边的Row和Column是需要选择的变量。Execute后就可以下载病例数据了，最后出现的是表格形式，可以将表格的数据直接复制粘贴到excel里此数据是汇总数据，也可以按照个体数据下载，但数据量大，下载时间会比较长。通过SEER数据库（优点是数据量大）的可获得性，结合统计学分析和数据挖掘等就可以撰写和发表论文了。PartII基因表达数据预处理分析

（一）基因芯片数据cDNA微阵列芯片荧光信号定性信息提取：P/A/M（Present/Absent/Marginal）定量信息提取：基于探针集汇总后的基因水平的荧光信号强度值对芯片数据做对数化转换后，数据可近似正态分布。对数转换前对数转换后数据转换数据过滤数据过滤的目的是去除表达水平是负值或很小的数据或者明显的噪声数据。过闪耀现象物理因素导致的信号污染杂交效能低点样问题其他缺失值填补1.数据缺失类型非随机缺失基因表达丰度过高或过低。随机缺失与基因表达丰度无关，数据补缺主要针对随机缺失情况。2.数据补缺方法（1）简单补缺法missingvalues=0expressionmissingvalues=1expressionmissingvalues=row（gene）averagemissingvalues=column（array）average（2）k近邻法选择与具有缺失值基因的k个邻居基因用邻居基因的加权平均估计缺失值数据标准化片内标化（within-slidenormalization）方法（局部加权回归）平行试验数据的标准化（数据点排序，重复试验的平均值）片间标化（multiple-slidenormalization）（平均数，中位数标准化）1）倍数法实验条件下的表达值对照条件下的表达值通常以2倍差异为阈值，判断基因是否差异表达差异表达基因的判断Fold_change2）t检验法

运用t检验法可以判断基因在两不同条件下的表达差异是否具有显著性

3）SAM法（significanceanalysisofmicroarrays）

确定差异表达基因阈值运用该阈值，统计在值中超过该阈值的假阳性基因个数，估计假阳性发现率FDR值。调整FDR值的大小得到差异表达基因。（二）RNA-seq测序技术及数据分析1）RNA-seq数据的主要数据分析流程RNA-seq数据的原始读数（rawreads）的质控包括序列质量分析、GC含量分析、测序误差分析或污染分析等。质控主要针对重复序列、错误比对以及碱基质量的校正，R软件包NOISeq及EDASeq等都可以绘制出相应的质控图。这样经过处理后的比对结果能够有效地减少了测序和比对过程中产生的假阳性和假阴性。2）RNA-seq数据的质控数据过滤和质控后，首先基本的分析是将reads映射到参考基因组或全基因组进行比对，进而推断和发现新的转录本，包括识别可变剪切、变异、非编码RNA和基因融合等。如果没有参考注释，则测序数据还需要进行拼接和重组。3）RNA-seq数据的比对RPKM（ReadsPerKilobasesperMillionreads）和FPKM（FragmentsPerKilobasesperMillionreads）常用来计算RNA-seq技术的基因表达水平。由于RNA-seq是基于reads量化获得的基因表达水平，因此差异表达基因的分析方法是基于泊松或负二项分布的离散型分布。常用的R分析软件包有DESeq、DESeq2、edgeR和baySeq等。筛选出来的差异表达基因可以进一步作基因的功能分析，如GO（GeneOntology）或KEGG（KyotoEncyclopediaofGenesandGenomes）富集分析。4）将Reads转换为基因表达水平5）基于RNA-seq数据的差异表达基因分析统计基因或转录本对应的读段（reads）计数，然后对读段计数进行标准化，使样本间和样本内的表达水平能够进行精确比较。最常见的一个指标是RPKM，每百万读段中来自某一基因每千碱基长度的读段数目：当前RPKM方法是最为常用的基因表达水平估算方法。进一步对标准化后读段计数分布进行统计学模型拟合，利用统计学检验评估基因的差异表达，得到相应的P值和差异倍数（foldchange），并完成多重检验校正，最后根据特定阈值（例如FDR<0.05）提取显著差异表达的基因。常用的R软件包有DEseq、DEGseq和baySeq等。

RNA-seq差异表达基因的筛选其中，第一列表示转录本名称；第二列baseMean表示归一化后该转录本在所有样本中的平均表达值；第三列log2FoldChange表示以2为底的对数化的倍数变化（倍数变化为两个条件下的均值之比）；第四列表示标准误；第五列表示负二项分布检验的Wald统计量；第六列表示的是负二项分布检验获得的p值；第七列表示的是调整后的p值，一般“padj”的阈值常设为0.05或0.01。adjustedP-value当前，已开发了很多R软件包,可根据数据类型（基因芯片microarray数据或RNA-seq数据）选择差异表达基因的分析方法，我们会在后面进行介绍。PartIII基因集功能富集分析

进行基因集功能富集分析的原因富集分析方法通常是分析一组基因在某个功能结点上是否过出现（over-presentation）。这个原理可以由单个基因的注释分析发展到大基因集合的成组分析。由于分析的结论是基于一组相关的基因，而不是根据单个基因，所以富集分析方法增加了研究的可靠性，同时也能够识别出与生物现象最相关的生物过程。GO富集分析中常用的统计方法有超几何分布、Fisher精确检验等。超几何分布：m表示人类所有的基因数，n表示筛选出的疾病相关的风险基因数（用户提供的基因集合中的基因数），t表示注释到GO结点的所有基因数，r为n和t的交集，表示注释到GO结点中的疾病风险基因数（一）GO富集分析Fisher精确检验：基因与GO结点的关系疾病风险基因数人类基因数映射到GO结点ab未映射到GO结点cdFisher精确检验公式为：这里以目前应用较为广泛的DAVID为例对基因集进行具体分析。DAVID是一个综合工具，不但提供基因富集分析，还提供基因间ID的转换、基因功能的分类等。三、富集应用分析实例

DAVID应用工具首页（）点击“FunctionalAnnotation”后，第一步为提交基因集，选择基因标识名和基因集类型；第二步得到注释结果摘要，包括多种注释数据；然后选择感兴趣的注释内容得到富集分析结果。点击“FunctionalAnnotation”，首先需要上传基因列表（UploadGeneList）。在Step1中可以通过粘贴基因列表（Pastealist）或者选择保存好的基因列表文件（ChooseFromaFile）上传。在Step2的基因类别识别（SelectIdentifier）中可以根据提供的基因名称进行选择。基因识别类别中包括Affymeterix、entrez、ensemble、uniprot、genesymbol等近40种基因识别ID，这里我们选择官方的基因名称识别（OFFICIAL_GENE_SYMBOL）。在Step3的列表类型中（ListType），有两种选择：一种是将上传的基因作为要分析的基因集合（GeneList），另一种是将上传的基因作为背景基因（Background）。这里我们勾选“GeneList”。Step4是提交基因列表（SubmitList）。点击提交后就可以看到富集分析结果了。点击已勾选的“GOTERM_BP_FAT”对应的灰色按钮“Chart”，可以得到如下图的输出结果：GO富集分析输出结果界面中包括GO条目的描述，富集分析的P值和Benjamini-Hochberg校正后的p值等。点击“DownloadFile”可将输出的结果保存。通路分析是现在经常被使用的芯片数据基因功能分析法。通路分析法利用的资源是许多已经研究清楚的基因之间的相互作用，即生物学通路。研究者可以把表达发生变化的基因集导入通路分析软件中，进而得到变化的基因都存在于哪些已知通路中，并通过统计学方法计算哪些通路与基因表达的变化最为相关。（二）通路富集分析1)KEGG通路数据库KEGG（KyotoEncyclopediaofGenesandGenomes）（）是系统分析基因功能、基因组信息的数据库，它整合了基因组学、生物化学以及系统功能组学的信息，有助于研究者把基因及表达信息作为一个整体进行研究。

KEGG不仅提供了所有可能的代谢通路，还对催化各步反应的酶进行了全面的注解，包含其氨基酸序列以及到PDB数据库的链接等。此外，KEGG还提供基于Java的图形工具访问基因组图谱、比较基因组图谱和操作表达图谱以及其他序列比较、图形比较和通路计算的工具。因此，KEGG数据库是进行生物体内代谢分析、代谢网络分析等研究的强有力工具之一。KEGG中的葡萄糖代谢通路KEGG中的非代谢通路（信号传导通路）KEGG通路中，结点间的各种关系查看乳腺癌基因BRCA1注释到哪些通路，在空白框里输入672red（672是BRCA1的EntrezID号，red表示在通路中该基因用红色表示）。点击Exec按钮执行搜索过程。搜索基因注释到哪些通路：搜索出了相关通路点击第三条通路（hsadd03460），得到该通路的通路图。2)其他通路数据库简介（1）Biocarta：Biocarta通路数据库从分子的关系角度描绘了一个网络图模型，通过不断整合蛋白质组信息迅速发展壮大起来。该数据库与KEGG数据库齐名，尽管其数据量目前可能没有KEGG的多，但在一些通路相关的研究中，这个数据库也和KEGG一道用来研究基因的功能，注释信息及进行富集分析等。（2）Reactome：Reactome是信号通路数据库，是一个汇集了由专家撰写，经同行评阅的有关人体内各项反应及生物学路径的文章的数据库。数据库目前覆盖了UniProt数据库中两万个经人工注释过的人类蛋白质中大约70%以上的蛋白质，对主要人类生物学研究领域，如细胞凋亡、DNA复制、转录、碳水化合物代谢途径等进行了注释。相对于KEGG数据库，Reactome是一个改良的搜索及数据挖掘工具，可以简化与生物学途径相关的数据搜索。（3）Biocyc：Biocyc数据库提供了上千种生物体基因组和代谢通路的参考信息。截止到2016年12月，Biocyc已融合了9300个数据库。Biocyc也包含了很多如搜索、可视化、比较和分析基因组和通路信息的软件工具。（4）GeneDB：GeneDB数据库是关于原核和真核病原体及相关生物体的基因组数据库。该数据库提供了基因组序列和注释数据，同时也整合了其他的数据库资源并发展了数据库驱动的注释工与GO功能富集分析一样，疾病通路的富集分析常采用的方法也是超几何检验或Fisher精确检验。超几何检验的公式如下：m表示人类所有基因数；t表示通路中所有的基因数；n表示疾病风险基因数（用户提供的基因集合中的基因数）以及r表示注释到通路中的疾病风险基因数（疾病风险基因和通路基因的交集）。3)KEGG通路富集分析如果采用Fisher精确检验，则表格可设计为四格表的形式：基于Fisher精确检验的通路富集分析基因与通路的关系疾病风险基因数人类基因数在通路中(Inpathway)ab不在通路中(Notinpathway)cdFisher精确检验公式为：DAVID软件也可用于通路的富集分析。方法与步骤与GO功能富集分析一致，提交基因列表：点击“KEGG_PATHWAY”对应的“Chart”按钮，可得到输出结果：KEGG富集分析输出结果界面中包括KEGG通路名称，富集分析的P值和Benjamini-Hochberg校正后的p值等。点击“DownloadFile”可将输出的结果保存。PartIV案例：肿瘤基因表达数据分析

R的CancerSubtypes软件包应用Bioconductor的安装语句source("")biocLite("CancerSubtypes")数据导入及数据预处理（数据正态化，数据的缺失值填补）聚类热图的绘制提取差异表达基因绘制火山图应用非负矩阵分解聚类法进行肿瘤亚型分析（肿瘤亚型比较，绘制亚型聚类热图和Silhouette图）（一）数据格式read.table("d:\\genematrix.csv",header=TRUE,sep=",")->aread.table("d:\\genematrixnohead.csv",header=FALSE,sep=",")->datadata<-as.matrix(data)rownames(data)<-a[,1]colnames(data)<-c(rep("D",60),c(rep("N",15)))data[1:10,1:10]（二）将整理好的excel数据导入library(CancerSubtypes)（三）检查数据的分布data.checkDistribution(data)(四）缺失值填补data=data.imputation(data,fun="median")（五）数据正态化data=data.normalization(data,type="feature_Median",log2=FALSE)type类型：feature_Median中位正态化feature_Mean均值正态化feature_zscore按基因标化sample_zscore按样本标化log2如果是TRUE,做以2为底的对数转换（常用于RNAseq数据）（六）绘制双向聚类图library(gplots)heatmap.2(data,col=redgreen)（七）求差异表达基因T_data<-data[,1:60]N_data<-data[,61:75]result=DiffExp.limma(Tumor_Data=T_data,Normal_Data=N_data,topk=NULL,RNAseq=FALSE)RNAseq=TRUE表示用RNAseq表达，RNAseq=FALSE表示microarrray数据。write.csv(result,file="d:\\result.csv")差异表达基因结果：（八）用输出结果绘制火山图

数据格式如下：library(ggplot2)read.table("d:\\volcano.csv",header=TRUE,sep=",")->datathreshold<-as.factor(abs(data$logFC)>=1.5&data$Pvalue<0.05)ggplot(data,aes(x=logFC,y=-log10(Pvalue),colour=threshold))+geom_point(alpha=0.4,size=1.75)+geom_vline(xintercept=c(-1,1),lty=4,col="grey",