從TCGA數(shù)據(jù)庫下載EBV(+)和EBV(-)的原始測序fastq數(shù)據(jù)，有25個EBV(+)和260個EBV(-)樣本。然后構(gòu)建一個臨床特征表格，包括性別、年齡、腫瘤部位、組織學(xué)病理學(xué)診斷等等，接著用R語言的dist()函數(shù)計算兩組樣本間臨床特征的歐式距離，找出距離最近的樣本進(jìn)行匹配。

因為有幾個EBV(-)樣本跟EBV(+)樣本的距離是一樣的，所以最終選到了20個EBV(-)樣本，來跟那25個EBV(+)樣本配對。這樣每個配對之間的距離，就比25個EBV(+)對260個EBV(-)之間的平均距離小多了，控制了混雜因素，使組間具有可比性。

該文章附件中的表格展示了配對樣本的一部分，左邊4列分別為EBV(+)樣本編號、組織學(xué)診斷、EBV(-)樣本編號、組織學(xué)診斷。右邊2列數(shù)字就是歐式距離，左邊的是EBV(+)樣本與它所匹配的EBV(-)樣本的距離，右邊的是該EBV(+)樣本與原260個EBV(-)樣本的平均距離，可見前者是比后者要小多，消除了樣本間差異的偏倚因素。

接下來對EBV基因進(jìn)行篩選，從原始的88個基因中，將在超過5個樣本中原始計數(shù)（raw count）為0的基因去除，剩下19個基因，這就獲得了EBV相關(guān)的特征分子標(biāo)志。

部分EBV基因在EBV(+)樣本中的表達(dá)

然后提取人類基因組的特征分子標(biāo)志。用Bioconductor的DEseq2包找出兩組樣本間差異表達(dá)的基因（DEx gene），這里就要用到Bonferroni法對P值進(jìn)行校正，原始P值達(dá)到1E-6，即校正后的P < 0.05，才算有統(tǒng)計學(xué)意義。共找到939個DEx基因，其中189個上調(diào)，750個下調(diào)。

再用Bioconductor中的GAGE (Generally Applicable Gene-set Enrichment for Pathway Analysis)分析法，找到兩組樣本中差異表達(dá)的通路，共計29個，在EBV(+)樣本中全是上調(diào)的。

同時用R語言的MEGENA包構(gòu)建那些DEx基因的共表達(dá)網(wǎng)絡(luò)，找到它們的聚類基因模塊（Gene modules）和其中的關(guān)鍵節(jié)點(diǎn)基因（hub genes）。前者modules是指一組具有較強(qiáng)相關(guān)性的基因，而hub則是指modules中具有全局性影響的樞紐基因。這一步是利用人類基因組中本來就有的交互網(wǎng)絡(luò)數(shù)據(jù)，來減少數(shù)據(jù)維度，聚焦到關(guān)鍵信息。這一步找到了91個節(jié)點(diǎn)基因和27個基因模塊 (一個節(jié)點(diǎn)基因可出現(xiàn)在不同的模塊中)。

然而，這些數(shù)據(jù)維度還是高，怎么進(jìn)一步抓取核心要素呢？作者又用了主成分分析法（PCA）提取其中起主要作用的信息，他們只提取了第一個主成分（PC1）來代表該模塊或通路。這樣，網(wǎng)絡(luò)構(gòu)建和相關(guān)性分析的素材提取完畢，最終獲得19個EBV基因、91個人類胃癌相關(guān)的hub基因、27個基因模塊和29個通路（的PC1）。

對于生信有基礎(chǔ)的同學(xué)，可按照上述方法利用工具自行探索，零基礎(chǔ)的同學(xué)我們還有后續(xù)教學(xué)。

單因素相關(guān)分析就是要從19個EBV基因、91個人類胃癌相關(guān)的hub基因、27個基因模塊和29個通路中再進(jìn)一步獲取相關(guān)性最高的組合，采用配對的Pearson相關(guān)性檢驗對三組數(shù)據(jù)分別進(jìn)行分析：1）基因模塊的PC1和EBV基因表達(dá)數(shù)據(jù)；2）差異通路的PC1和EBV基因表達(dá)數(shù)據(jù)；3）hub基因和EBV基因表達(dá)數(shù)據(jù)。用置換檢驗法獲取相關(guān)系數(shù)（用C++）。

至于相關(guān)系數(shù)的顯著性，則沒有采用通用的p<0.05，因為它假設(shè)兩個變量都是正態(tài)分布，這對RNA-seq數(shù)據(jù)來說通常不現(xiàn)實。所以研究者采用p < 0.1 作為統(tǒng)計學(xué)顯著性的門檻。

這一步得到了7個EBV基因和12個模塊的PC1呈顯著相關(guān)，其中相關(guān)性最強(qiáng)的是BLAF1基因和12號模塊，相關(guān)系數(shù)達(dá)到0.662，如下表。

3個EBV基因和4個人類胃癌的hub基因相關(guān)，其中最強(qiáng)的是LMP-1和Clorf115，相關(guān)系數(shù)0.754。

又有2個EBV基因和4個通路的PC1相關(guān)，其中BLAF4和4個通路都相關(guān)，最強(qiáng)的是BLAF4和磷脂酰肌醇信號通路，相關(guān)系數(shù)達(dá)到0.709。

好了，單因素分析這步到這里做完。

這步要用到稀疏典型相關(guān)分析法（sparse Canonical Correlation Analysis, sCCA）。典型相關(guān)分析法（CCA）是分析兩個數(shù)據(jù)矩陣的經(jīng)典方法，不需要降維提取主成分。但普通CCA只適用于行比列多的矩陣，即樣本比基因多，這情況顯然不對，是樣本少基因多，所以要用sCCA。這步用R語言的PMA包來做。

其中用CCA()函數(shù)可得到每個元素（比如每個EBV基因）的典型系數(shù)（canonical coefficient），它表示該元素對整個矩陣典型相關(guān)系數(shù)的貢獻(xiàn)，如果某個元素的典型系數(shù)非零，則表示它對全局相關(guān)性有很大貢獻(xiàn)，被選為核心（essential）基因。

這步得到22個基因模塊與19個EBV基因計數(shù)矩陣相關(guān)（P< 0.1，理由同上），如下表。其中核心EBV基因為LMP-1，BALF1，BALF2，BARF1，BNRF1，LF1和BZLF1。

此外，人類hub基因計數(shù)矩陣和EBV基因計數(shù)矩陣的典型相關(guān)系數(shù)是0.806。對這個相關(guān)性貢獻(xiàn)最大的核心EBV基因和核心人類hub基因就不一一列出了。所有的差異表達(dá)的通路都和EBV基因計數(shù)矩陣達(dá)到了顯著相關(guān)的典型系數(shù)。最強(qiáng)幾個的如下表所示。

至此，數(shù)據(jù)分析的主要步驟就做完了。但文章沒完呀，還要做個總結(jié)，比較一下Pearson和sCCA兩種相關(guān)性分析的結(jié)果。咱們接著看。

先總結(jié)一下EBV基因的情況。如下表，各個重要EBV基因在Pearson和sCCA分析中，分別與人類胃癌基因模塊、hub基因、通路相關(guān)的情況。

人類重要基因模塊在兩種相關(guān)性分析中的情況則如下表：

比較有意思的是模塊5和模塊12。模塊5在DAVID數(shù)據(jù)庫的注釋中，有膜、跨膜、離子通道等關(guān)鍵詞；而模塊12則和饑餓激素通路相關(guān)，在前人的研究中，它可能會促進(jìn)胃腸道和胰腺疾病的惡化。由此，文章數(shù)據(jù)分析的統(tǒng)計學(xué)結(jié)果終于與生物學(xué)意義產(chǎn)生了交匯！

下圖展示了模塊5和EBV基因相互作用的關(guān)系圖，作為上述相關(guān)性的一個示例：

再來看看重要的人類hub基因，它們是C1orf115，CNTD2和VANGL2。它們在Pearson分析中，和EBV的BALF1和LMP-1相關(guān)，在sCCA分析中也被選為核心基因。下圖展示了人類hub基因和EBV基因的互作網(wǎng)絡(luò)。

人類重要的通路則是凋亡、溶酶體、Jak-STAT信號通路和磷脂酰肌醇信號通路，它們在Pearson和sCCA分析中都與EBV的BALF4和/或BALF5基因相關(guān)。值得注意的是，BALF4基因與Jak-STAT信號通路、磷脂酰肌醇信號通路中的大多數(shù)基因都呈現(xiàn)正相關(guān)。

總體而言，在第一步差異表達(dá)分析中所找到的DEx基因、模塊和通路，與篩選出來的EBV基因，大多都有相關(guān)性，這也驗證了前人的觀察和實驗的結(jié)果，即EBV和胃癌確實存在基因?qū)用娴南嗷プ饔谩?/p>

生信和統(tǒng)計分析完成后，得到了一堆具有統(tǒng)計學(xué)相關(guān)性的基因、通路，這時候需要找出它們的生物學(xué)相關(guān)性，才是拔高文章水平的點(diǎn)睛之筆。作者分析了他們找到的每個具有統(tǒng)計學(xué)相關(guān)性的元素在既往文獻(xiàn)報道中的相關(guān)性，文獻(xiàn)檢索的工作相當(dāng)細(xì)致。例如LMP-1編碼的潛伏膜蛋白1可導(dǎo)致細(xì)胞生長失調(diào)，是一個致癌蛋白；BALF1編碼抗凋亡蛋白Bcl2同源體，它的表達(dá)也和幾種癌癥的發(fā)生有相關(guān)性，等等等等。

當(dāng)然除了驗證前人的研究，BALF2可謂一個新發(fā)現(xiàn)，之前沒有報道過它和癌癥有什么關(guān)系，本文是第一次提出，如果能夠經(jīng)過實驗驗證這文章絕非五分這般平庸啊。

TCGA的挖掘潛力日益受到科研界的重視，搞腫瘤的這一優(yōu)勢得天獨(dú)厚，應(yīng)善加利用。生信的技能學(xué)習(xí)是一個由淺入深的漸進(jìn)過程，解螺旋將推出系列課程，今天放出第一波，由我們的合作伙伴GCBI制作的TCGA數(shù)據(jù)庫基礎(chǔ)應(yīng)用教學(xué)，課件免費(fèi)下載，在服務(wù)號完成分享任務(wù)后獲取16mins視頻。

懶癌晚期患者，請加入解螺旋鉆石會員，不但能夠?qū)W習(xí)臨床研究，基金申請和SCI寫作的深度課程，當(dāng)下還可享受“TCGA基因萬能篩”活動，1000元，給你用生信方式篩選一個靠譜的研究靶標(biāo)基因，實驗驗證不成功免費(fèi)重篩，把專業(yè)的事交給專業(yè)的人辦。

下面是“TCGA基因萬能篩”活動獲得的靶標(biāo)基因示例：

針對研究癌癥卻沒有方向的會員，整理TCGA數(shù)據(jù)庫中特定癌癥的癌/癌旁數(shù)據(jù)，包括miRNA或mRNA數(shù)據(jù)，篩選癌/癌旁差異基因，并結(jié)合多種生物信息學(xué)分析內(nèi)容（包括差異基因功能、信號通路、生存以及基因網(wǎng)絡(luò)分析等），為會員研究特定癌癥的發(fā)生發(fā)展機(jī)制提供有價值的靶基因及相應(yīng)的理論基礎(chǔ)。

差異分析結(jié)果

差異基因功能分析

差異基因信號通路分析

基于KEGG數(shù)據(jù)庫，對差異基因利用Fisher精確檢驗和卡方檢驗，把目標(biāo)基因參與的Pathway進(jìn)行顯著性分析，按照P value<0.05進(jìn)行篩選，得到顯著性的Pathway。

生存分析

利用單基因在所有樣本中的表達(dá)情況，按照中位數(shù)進(jìn)行分組，利用時序檢驗(log-rank)檢驗進(jìn)行比較兩組的總生存期的差異性，篩選條件為：p<0.05，研究基因表達(dá)對預(yù)后的影響。

調(diào)控關(guān)系

每個基因與基因之間都有非常多的調(diào)控關(guān)系，我們從基因相關(guān)的轉(zhuǎn)錄因子、miRNA、lncRNA和上下游相關(guān)基因的四個角度來展示基因與基因間的關(guān)系。

基因CFH中，已經(jīng)驗證的上游靶向結(jié)合的miRNA為has-miR-146a-5p。在基因網(wǎng)絡(luò)中，CFH能夠抑制基因CFB和C3的表達(dá)。

文獻(xiàn)報道中和CFH相關(guān)的lncRNA的71個。

轉(zhuǎn)錄因子預(yù)測

根據(jù)位置結(jié)合關(guān)系，預(yù)測結(jié)合CFH啟動子區(qū)域（轉(zhuǎn)錄起始位點(diǎn)上游2000bp，下游500bp）的轉(zhuǎn)錄因子。并依據(jù)Transfac數(shù)據(jù)庫的預(yù)測分值和COSMIC數(shù)據(jù)庫以及dbSNP數(shù)據(jù)庫是否存在對應(yīng)注釋信息整合出推薦度；以推薦度來表示預(yù)測出來的轉(zhuǎn)錄因子的科研價值，推薦度越高越好。