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

因?yàn)橛袔讉€(gè)EBV(-)樣本跟EBV(+)樣本的距離是一樣的，所以最終選到了20個(gè)EBV(-)樣本，來(lái)跟那25個(gè)EBV(+)樣本配對(duì)。這樣每個(gè)配對(duì)之間的距離，就比25個(gè)EBV(+)對(duì)260個(gè)EBV(-)之間的平均距離小多了，控制了混雜因素，使組間具有可比性。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

這步得到22個(gè)基因模塊與19個(gè)EBV基因計(jì)數(shù)矩陣相關(guān)（P< 0.1，理由同上），如下表。其中核心EBV基因?yàn)長(zhǎng)MP-1，BALF1，BALF2，BARF1，BNRF1，LF1和BZLF1。

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

差異分析結(jié)果

差異基因功能分析

差異基因信號(hào)通路分析

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

生存分析

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

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

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

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

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

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

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