第二講基因表達數據分析

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）定量信息提?。夯谔结樇瘏R總后的基因水平的熒光信號強度值對芯片數據做對數化轉換后，數據可近似正態分布。對數轉換前對數轉換后數據轉換數據過濾數據過濾的目的是去除表達水平是負值或很小的數據或者明顯的噪聲數據。過閃耀現象物理因素導致的信號污染雜交效能低點樣問題其他缺失值填補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）中可以根據提供的基因名稱進行選擇?；蜃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”可將輸出的結果保存。通路分析是現在經常被使用的芯片數據基因功能分析法。通路分析法利用的資源是許多已經研究清楚的基因之間的相互作用，即生物學通路。研究者可以把表達發生變化的基因集導入通路分析軟件中，進而得到變化的基因都存在于哪些已知通路中，并通過統計學方法計算哪些通路與基因表達的變化最為相關。（二）通路富集分析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",