二講基因表達數(shù)據(jù)分析

第二講基因表達數(shù)據(jù)分析

PartI基因測定平臺及數(shù)據(jù)庫介紹(GEO、TCGA、SEER數(shù)據(jù)庫的下載）PartII基因芯片數(shù)據(jù)的預處理分析PartIII基因功能分析及DAVID網(wǎng)絡工具應用PartIV綜合案例：腫瘤基因表達數(shù)據(jù)分析PartI基因表達測定平臺及數(shù)據(jù)庫介紹近20年來三種不同高通量基因表達測定技術的應用趨勢基因表達測定平臺與數(shù)據(jù)庫

1.cDNA芯片2.Affymetrix芯片

3.RNA-Seq技術對沒有已知參考基因組信息的非模式生物，也可測定轉(zhuǎn)錄信息；RNA-Seq技術可以測定轉(zhuǎn)錄邊界的精度達到一個堿基；RNA-Seq可以用來研究復雜的轉(zhuǎn)錄關系；RNA-Seq可以同時測定序列的變異；RNA-Seq背景信號很小，測定的動態(tài)范圍很大。常用基因表達數(shù)據(jù)庫名稱數(shù)據(jù)庫內(nèi)容GeneExpressionOmnibus（GEO）目前最常用的基因表達數(shù)據(jù)（NCBI）ExpressionAtlas歐洲生物信息學中心的基因表達數(shù)據(jù)庫SMDStanford基因表達數(shù)據(jù)庫RNA-SeqAtlas正常組織的基因表達譜數(shù)據(jù)GEPdb基因型、表型和基因表達關系GXD老鼠發(fā)育基因表達信息EMAGE老鼠胚胎的時空表達信息AGEMAP老鼠老化的基因表達數(shù)據(jù)常用基因表達數(shù)據(jù)庫

下載膀胱癌數(shù)據(jù)TCGA數(shù)據(jù)庫的下載推薦UCSCXena下載TCGA數(shù)據(jù)

下載膀胱癌數(shù)據(jù)用RPKM值矩陣數(shù)據(jù)注釋數(shù)據(jù)基因表達矩陣數(shù)據(jù)基因注釋數(shù)據(jù)表型數(shù)據(jù)總生存時間無復發(fā)生存時間SEER數(shù)據(jù)庫TheSurveillance,Epidemiology,andEndResultsSEER數(shù)據(jù)庫（）美國國立癌癥研究所“監(jiān)測、流行病學和結果數(shù)據(jù)庫”，是北美最具代表性的大型腫瘤登記注冊數(shù)據(jù)庫之一，收集了大量循證醫(yī)學的相關數(shù)據(jù)，為臨床醫(yī)師的循證實踐及臨床醫(yī)學研究提供了系統(tǒng)的證據(jù)支持和資料。SEER數(shù)據(jù)庫及軟件下載SEER數(shù)據(jù)庫不像其他數(shù)據(jù)庫，SEER需要注冊，而且需要簽份保證書，得到官方認可才可以得到一個用戶名和賬號，用于數(shù)據(jù)和軟件的下載安裝。

注冊入口：填寫信息提交后將迅速收到郵件，進行確認簽署協(xié)議，發(fā)郵件后等半天~1天就會收到SEER發(fā)過來的用戶名和密碼。SEER數(shù)據(jù)庫SEERStat軟件雙擊SEERstat圖標鏈接下載的數(shù)據(jù)庫data是選擇數(shù)據(jù)來源的范圍，SEER數(shù)據(jù)庫最近更新到2015年的隨訪數(shù)據(jù)，所以盡量選最新的數(shù)據(jù)集。selection是選擇病例篩選的條件。點擊edit，然后會有很多選項，選擇你需要的病例范圍、腫瘤類型等等。SEER數(shù)據(jù)庫有許多自己定義的簡寫和縮寫，比如metsatdx代表遠處轉(zhuǎn)移，轉(zhuǎn)移部位和病理類型都用代號表示。建議官網(wǎng)下載“CScodinginstructions”的文件，里面詳細列出了該數(shù)據(jù)庫使用的每個coding的解釋table是選擇需要的臨床信息。左邊的Row和Column是需要選擇的變量。Execute后就可以下載病例數(shù)據(jù)了，最后出現(xiàn)的是表格形式，可以將表格的數(shù)據(jù)直接復制粘貼到excel里此數(shù)據(jù)是匯總數(shù)據(jù)，也可以按照個體數(shù)據(jù)下載，但數(shù)據(jù)量大，下載時間會比較長。通過SEER數(shù)據(jù)庫（優(yōu)點是數(shù)據(jù)量大）的可獲得性，結合統(tǒng)計學分析和數(shù)據(jù)挖掘等就可以撰寫和發(fā)表論文了。PartII基因表達數(shù)據(jù)預處理分析

（一）基因芯片數(shù)據(jù)cDNA微陣列芯片熒光信號定性信息提取：P/A/M（Present/Absent/Marginal）定量信息提?。夯谔结樇瘏R總后的基因水平的熒光信號強度值對芯片數(shù)據(jù)做對數(shù)化轉(zhuǎn)換后，數(shù)據(jù)可近似正態(tài)分布。對數(shù)轉(zhuǎn)換前對數(shù)轉(zhuǎn)換后數(shù)據(jù)轉(zhuǎn)換數(shù)據(jù)過濾數(shù)據(jù)過濾的目的是去除表達水平是負值或很小的數(shù)據(jù)或者明顯的噪聲數(shù)據(jù)。過閃耀現(xiàn)象物理因素導致的信號污染雜交效能低點樣問題其他缺失值填補1.數(shù)據(jù)缺失類型非隨機缺失基因表達豐度過高或過低。隨機缺失與基因表達豐度無關，數(shù)據(jù)補缺主要針對隨機缺失情況。2.數(shù)據(jù)補缺方法（1）簡單補缺法missingvalues=0expressionmissingvalues=1expressionmissingvalues=row（gene）averagemissingvalues=column（array）average（2）k近鄰法選擇與具有缺失值基因的k個鄰居基因用鄰居基因的加權平均估計缺失值數(shù)據(jù)標準化片內(nèi)標化（within-slidenormalization）方法（局部加權回歸）平行試驗數(shù)據(jù)的標準化（數(shù)據(jù)點排序，重復試驗的平均值）片間標化（multiple-slidenormalization）（平均數(shù)，中位數(shù)標準化）1）倍數(shù)法實驗條件下的表達值對照條件下的表達值通常以2倍差異為閾值，判斷基因是否差異表達差異表達基因的判斷Fold_change2）t檢驗法

運用t檢驗法可以判斷基因在兩不同條件下的表達差異是否具有顯著性

3）SAM法（significanceanalysisofmicroarrays）

確定差異表達基因閾值運用該閾值，統(tǒng)計在值中超過該閾值的假陽性基因個數(shù)，估計假陽性發(fā)現(xiàn)率FDR值。調(diào)整FDR值的大小得到差異表達基因。（二）RNA-seq測序技術及數(shù)據(jù)分析1）RNA-seq數(shù)據(jù)的主要數(shù)據(jù)分析流程RNA-seq數(shù)據(jù)的原始讀數(shù)（rawreads）的質(zhì)控包括序列質(zhì)量分析、GC含量分析、測序誤差分析或污染分析等。質(zhì)控主要針對重復序列、錯誤比對以及堿基質(zhì)量的校正，R軟件包NOISeq及EDASeq等都可以繪制出相應的質(zhì)控圖。這樣經(jīng)過處理后的比對結果能夠有效地減少了測序和比對過程中產(chǎn)生的假陽性和假陰性。2）RNA-seq數(shù)據(jù)的質(zhì)控數(shù)據(jù)過濾和質(zhì)控后，首先基本的分析是將reads映射到參考基因組或全基因組進行比對，進而推斷和發(fā)現(xiàn)新的轉(zhuǎn)錄本，包括識別可變剪切、變異、非編碼RNA和基因融合等。如果沒有參考注釋，則測序數(shù)據(jù)還需要進行拼接和重組。3）RNA-seq數(shù)據(jù)的比對RPKM（ReadsPerKilobasesperMillionreads）和FPKM（FragmentsPerKilobasesperMillionreads）常用來計算RNA-seq技術的基因表達水平。由于RNA-seq是基于reads量化獲得的基因表達水平，因此差異表達基因的分析方法是基于泊松或負二項分布的離散型分布。常用的R分析軟件包有DESeq、DESeq2、edgeR和baySeq等。篩選出來的差異表達基因可以進一步作基因的功能分析，如GO（GeneOntology）或KEGG（KyotoEncyclopediaofGenesandGenomes）富集分析。4）將Reads轉(zhuǎn)換為基因表達水平5）基于RNA-seq數(shù)據(jù)的差異表達基因分析統(tǒng)計基因或轉(zhuǎn)錄本對應的讀段（reads）計數(shù)，然后對讀段計數(shù)進行標準化，使樣本間和樣本內(nèi)的表達水平能夠進行精確比較。最常見的一個指標是RPKM，每百萬讀段中來自某一基因每千堿基長度的讀段數(shù)目：當前RPKM方法是最為常用的基因表達水平估算方法。進一步對標準化后讀段計數(shù)分布進行統(tǒng)計學模型擬合，利用統(tǒng)計學檢驗評估基因的差異表達，得到相應的P值和差異倍數(shù)（foldchange），并完成多重檢驗校正，最后根據(jù)特定閾值（例如FDR<0.05）提取顯著差異表達的基因。常用的R軟件包有DEseq、DEGseq和baySeq等。

RNA-seq差異表達基因的篩選其中，第一列表示轉(zhuǎn)錄本名稱；第二列baseMean表示歸一化后該轉(zhuǎn)錄本在所有樣本中的平均表達值；第三列l(wèi)og2FoldChange表示以2為底的對數(shù)化的倍數(shù)變化（倍數(shù)變化為兩個條件下的均值之比）；第四列表示標準誤；第五列表示負二項分布檢驗的Wald統(tǒng)計量；第六列表示的是負二項分布檢驗獲得的p值；第七列表示的是調(diào)整后的p值，一般“padj”的閾值常設為0.05或0.01。adjustedP-value當前，已開發(fā)了很多R軟件包,可根據(jù)數(shù)據(jù)類型（基因芯片microarray數(shù)據(jù)或RNA-seq數(shù)據(jù)）選擇差異表達基因的分析方法，我們會在后面進行介紹。PartIII基因集功能富集分析

進行基因集功能富集分析的原因富集分析方法通常是分析一組基因在某個功能結點上是否過出現(xiàn)（over-presentation）。這個原理可以由單個基因的注釋分析發(fā)展到大基因集合的成組分析。由于分析的結論是基于一組相關的基因，而不是根據(jù)單個基因，所以富集分析方法增加了研究的可靠性，同時也能夠識別出與生物現(xiàn)象最相關的生物過程。GO富集分析中常用的統(tǒng)計方法有超幾何分布、Fisher精確檢驗等。超幾何分布：m表示人類所有的基因數(shù)，n表示篩選出的疾病相關的風險基因數(shù)（用戶提供的基因集合中的基因數(shù)），t表示注釋到GO結點的所有基因數(shù)，r為n和t的交集，表示注釋到GO結點中的疾病風險基因數(shù)（一）GO富集分析Fisher精確檢驗：基因與GO結點的關系疾病風險基因數(shù)人類基因數(shù)映射到GO結點ab未映射到GO結點cdFisher精確檢驗公式為：這里以目前應用較為廣泛的DAVID為例對基因集進行具體分析。DAVID是一個綜合工具，不但提供基因富集分析，還提供基因間ID的轉(zhuǎn)換、基因功能的分類等。三、富集應用分析實例

DAVID應用工具首頁（）點擊“FunctionalAnnotation”后，第一步為提交基因集，選擇基因標識名和基因集類型；第二步得到注釋結果摘要，包括多種注釋數(shù)據(jù)；然后選擇感興趣的注釋內(nèi)容得到富集分析結果。點擊“FunctionalAnnotation”，首先需要上傳基因列表（UploadGeneList）。在Step1中可以通過粘貼基因列表（Pastealist）或者選擇保存好的基因列表文件（ChooseFromaFile）上傳。在Step2的基因類別識別（SelectIdentifier）中可以根據(jù)提供的基因名稱進行選擇?；蜃R別類別中包括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”可將輸出的結果保存。通路分析是現(xiàn)在經(jīng)常被使用的芯片數(shù)據(jù)基因功能分析法。通路分析法利用的資源是許多已經(jīng)研究清楚的基因之間的相互作用，即生物學通路。研究者可以把表達發(fā)生變化的基因集導入通路分析軟件中，進而得到變化的基因都存在于哪些已知通路中，并通過統(tǒng)計學方法計算哪些通路與基因表達的變化最為相關。（二）通路富集分析1)KEGG通路數(shù)據(jù)庫KEGG（KyotoEncyclopediaofGenesandGenomes）（）是系統(tǒng)分析基因功能、基因組信息的數(shù)據(jù)庫，它整合了基因組學、生物化學以及系統(tǒng)功能組學的信息，有助于研究者把基因及表達信息作為一個整體進行研究。

KEGG不僅提供了所有可能的代謝通路，還對催化各步反應的酶進行了全面的注解，包含其氨基酸序列以及到PDB數(shù)據(jù)庫的鏈接等。此外，KEGG還提供基于Java的圖形工具訪問基因組圖譜、比較基因組圖譜和操作表達圖譜以及其他序列比較、圖形比較和通路計算的工具。因此，KEGG數(shù)據(jù)庫是進行生物體內(nèi)代謝分析、代謝網(wǎng)絡分析等研究的強有力工具之一。KEGG中的葡萄糖代謝通路KEGG中的非代謝通路（信號傳導通路）KEGG通路中，結點間的各種關系查看乳腺癌基因BRCA1注釋到哪些通路，在空白框里輸入672red（672是BRCA1的EntrezID號，red表示在通路中該基因用紅色表示）。點擊Exec按鈕執(zhí)行搜索過程。搜索基因注釋到哪些通路：搜索出了相關通路點擊第三條通路（hsadd03460），得到該通路的通路圖。2)其他通路數(shù)據(jù)庫簡介（1）Biocarta：Biocarta通路數(shù)據(jù)庫從分子的關系角度描繪了一個網(wǎng)絡圖模型，通過不斷整合蛋白質(zhì)組信息迅速發(fā)展壯大起來。該數(shù)據(jù)庫與KEGG數(shù)據(jù)庫齊名，盡管其數(shù)據(jù)量目前可能沒有KEGG的多，但在一些通路相關的研究中，這個數(shù)據(jù)庫也和KEGG一道用來研究基因的功能，注釋信息及進行富集分析等。（2）Reactome：Reactome是信號通路數(shù)據(jù)庫，是一個匯集了由專家撰寫，經(jīng)同行評閱的有關人體內(nèi)各項反應及生物學路徑的文章的數(shù)據(jù)庫。數(shù)據(jù)庫目前覆蓋了UniProt數(shù)據(jù)庫中兩萬個經(jīng)人工注釋過的人類蛋白質(zhì)中大約70%以上的蛋白質(zhì)，對主要人類生物學研究領域，如細胞凋亡、DNA復制、轉(zhuǎn)錄、碳水化合物代謝途徑等進行了注釋。相對于KEGG數(shù)據(jù)庫，Reactome是一個改良的搜索及數(shù)據(jù)挖掘工具，可以簡化與生物學途徑相關的數(shù)據(jù)搜索。（3）Biocyc：Biocyc數(shù)據(jù)庫提供了上千種生物體基因組和代謝通路的參考信息。截止到2016年12月，Biocyc已融合了9300個數(shù)據(jù)庫。Biocyc也包含了很多如搜索、可視化、比較和分析基因組和通路信息的軟件工具。（4）GeneDB：GeneDB數(shù)據(jù)庫是關于原核和真核病原體及相關生物體的基因組數(shù)據(jù)庫。該數(shù)據(jù)庫提供了基因組序列和注釋數(shù)據(jù)，同時也整合了其他的數(shù)據(jù)庫資源并發(fā)展了數(shù)據(jù)庫驅(qū)動的注釋工與GO功能富集分析一樣，疾病通路的富集分析常采用的方法也是超幾何檢驗或Fisher精確檢驗。超幾何檢驗的公式如下：m表示人類所有基因數(shù)；t表示通路中所有的基因數(shù)；n表示疾病風險基因數(shù)（用戶提供的基因集合中的基因數(shù)）以及r表示注釋到通路中的疾病風險基因數(shù)（疾病風險基因和通路基因的交集）。3)KEGG通路富集分析如果采用Fisher精確檢驗，則表格可設計為四格表的形式：基于Fisher精確檢驗的通路富集分析基因與通路的關系疾病風險基因數(shù)人類基因數(shù)在通路中(Inpathway)ab不在通路中(Notinpathway)cdFisher精確檢驗公式為：DAVID軟件也可用于通路的富集分析。方法與步驟與GO功能富集分析一致，提交基因列表：點擊“KEGG_PATHWAY”對應的“Chart”按鈕，可得到輸出結果：KEGG富集分析輸出結果界面中包括KEGG通路名稱，富集分析的P值和Benjamini-Hochberg校正后的p值等。點擊“DownloadFile”可將輸出的結果保存。PartIV案例：腫瘤基因表達數(shù)據(jù)分析

R的CancerSubtypes軟件包應用Bioconductor的安裝語句source("")biocLite("CancerSubtypes")數(shù)據(jù)導入及數(shù)據(jù)預處理（數(shù)據(jù)正態(tài)化，數(shù)據(jù)的缺失值填補）聚類熱圖的繪制提取差異表達基因繪制火山圖應用非負矩陣分解聚類法進行腫瘤亞型分析（腫瘤亞型比較，繪制亞型聚類熱圖和Silhouette圖）（一）數(shù)據(jù)格式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數(shù)據(jù)導入library(CancerSubtypes)（三）檢查數(shù)據(jù)的分布data.checkDistribution(data)(四）缺失值填補data=data.imputation(data,fun="median")（五）數(shù)據(jù)正態(tài)化data=data.normalization(data,type="feature_Median",log2=FALSE)type類型：feature_Median中位正態(tài)化feature_Mean均值正態(tài)化feature_zscore按基因標化sample_zscore按樣本標化log2如果是TRUE,做以2為底的對數(shù)轉(zhuǎn)換（常用于RNAseq數(shù)據(jù)）（六）繪制雙向聚類圖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數(shù)據(jù)。write.csv(result,file="d:\\result.csv")差異表達基因結果：（八）用輸出結果繪制火山圖

數(shù)據(jù)格式如下：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",