|
這是一份大三學生的課程實踐作業(yè)~~~ 找出膠質(zhì)細胞瘤特異性甲基化區(qū)域，為臨床診斷提供理論依據(jù)
可靠性 一、數(shù)據(jù)下載 首先進入TCGA下載數(shù)據(jù)GBM的RNA-seq和甲基化數(shù)據(jù)，從下表可見GBM共有172套RNA-seq數(shù)據(jù)以及437套DNA甲基化數(shù)據(jù)，由于TCGA提供Infinium HumanMethylation27 BeadChip和Infinium HumanMethylation450 BeadChip兩種芯片平臺的數(shù)據(jù)，為了避免后續(xù)不同芯片平臺間數(shù)據(jù)合并的困難，僅下載HumanMethylation450的芯片數(shù)據(jù)，共計154套。 圖表 1TCGA數(shù)據(jù)匯總 二、初步整理數(shù)據(jù) 使用TCGA-Assembler.2.0.5進行GBM數(shù)據(jù)批量下載與初步整理，并且繪制RNA-seq 基因表達量盒型 圖 以及甲基化芯片數(shù)據(jù)盒型圖 ，由于數(shù)據(jù)量較大，此處不貼圖。 三、整體可視化 首先對于甲基化數(shù)據(jù)，選取ID為TCGA.06.AABW.11A.31D.A368.05的數(shù)據(jù)，查看總體甲基化程度。由于每個位點真實情況只存在：甲基化/非甲基化兩種，所以對全部位點甲基化程度進行統(tǒng)計，也應該是大部分位點處于“完全甲基化”(Methylation state=1)和“完全非甲基化”(Methylation state=0)兩種狀態(tài)，下圖繪制了數(shù)據(jù)的頻數(shù)柱狀圖，可以明顯看出形狀處于“兩頭高，中間低”，反向說明芯片數(shù)據(jù)質(zhì)量較好。 圖表 2單個樣本CpG甲基化程度統(tǒng)計 接下來，對多個樣本繪制CpG甲基化程度小提琴圖，同一行是同一個病人，左邊樣本來源于Primary Solid Tumor，右邊樣本來源于Recurrent Solid Tumor，除了甲基化程度大部分分布于0和1附近外，還能看出來源于同一病人腫瘤的甲基化程度依舊會有略微差異。 TCGA barcode：https://wiki.nci./display/TCGA/TCGA+barcode 圖表 3小提琴圖 同樣的，對于RNA-seq數(shù)據(jù)也可以進行一些初步可視化，除了數(shù)據(jù)下載后繪制的盒型圖，亦可以進行PCA初步查看數(shù)據(jù)分布，下圖左為PCA陡坡圖，反映了第一主成分、第二主成分…等等所擁有信息量的比例，下圖右為使用PCA1和PCA2繪制的散點圖，可以發(fā)現(xiàn)5個正常樣本距離較近，從側(cè)面反映數(shù)據(jù)可信度較好。 最后，對于RNA-seq表達譜數(shù)據(jù)，使用系統(tǒng)聚類方法，繪制樹狀圖，可以發(fā)現(xiàn)5個正常樣本距離也是很近，數(shù)據(jù)質(zhì)量還行。 四、差異甲基化區(qū)域篩選 為了更加科學高效地篩選差異甲基化位點，參考bioconductor中甲基化芯片的分析流程，使用minfi包進行差異甲基化分析，得到差異甲基化位點。 http://www./help/workflows/methylationArrayAnalysis/ 在檢測的526733個CpG位點中，共有4927個CpG位點P值<0.01，且在腫瘤樣本中保持著甲基化程度高于0.7，對應2054個基因。 五、差異表達基因篩選 由于數(shù)據(jù)源自RNA-seq，最主流的分析方法當然是基于負二項分布模型的DESeq2包。 先用MA-plot查看差異表達基因大致分布。意外的是，圖形左側(cè)有大概七條線狀條紋，最初我懷疑這是sample之間有batch effect導致，需要用其他更好normalize的方法，后來用identify方法挨個找出每條線上的基因名及其對應的表達量，發(fā)現(xiàn)這些基因在172套樣本中表達量幾乎全為0，僅有一兩個樣本有一點點表達，這種數(shù)據(jù)的存在導致這些線狀條紋的產(chǎn)生。 圖表 4 MA-plot 然后， 選取p值最小的差異表達基因，繪制其在不同組間表達量，確實差異很顯著。 圖表 5表達量散點圖 接著，繪制差異表達基因在不同組間的表達量熱圖，正常樣本是圖片最左邊的五列，當然如果需要解釋具體的生物學問題，需要將聚類出來的每一類，將差異表達基因進行GO以及KEGG注釋，結(jié)合有關的生物學表型，探討其分子機制及意義 圖表 6表達量熱圖 最后，選取篩選條件為p值小于0.01且log2FoldChange<-2的差異表達基因，在腫瘤樣本低表達的基因，共計1657個 六、抑癌基因的獲取 在pubmed中查詢研究人員整理的tumor suppressor genes，果然在Nucleic Acids Res發(fā)表過TSGene數(shù)據(jù)庫，存儲了抑癌基因的列表。https://bioinfo./TSGene/download.cgi 下載全部1217個人類抑癌基因的列表。 All the 1217 human tumor suppressor genes with basic annotations
圖表 7 TSGene database 七、數(shù)據(jù)整合 對于甲基化數(shù)據(jù)中，腫瘤樣本高甲基化CpG附近的基因，RNA-seq中腫瘤樣本低表達的基因，以及TSGene數(shù)據(jù)庫中下載的抑癌基因列表，三者做overlap，找出特異性的候選靶標，為后續(xù)分析做準備，下圖為三者overlap的韋恩圖。
圖表 8數(shù)據(jù)整合韋恩圖 共計找出12個候選靶標基因。 八、靶標篩選 之前篩選選擇的單個CpG的差異甲基化，而實際臨床檢測應用時候，可能需要多個CpG作為對照，因此統(tǒng)計了12個候選靶標基因TSS前1.5kb內(nèi)所有CpG的甲基化程度，然后繪制熱圖，可以明顯發(fā)現(xiàn)，雖然當初用CpG的差異甲基化位點篩出來的基因都是腫瘤樣本高甲基化的，可是統(tǒng)計TSS前1.5kb內(nèi)所有CpG的甲基化程度，這些基因卻有很多在所有樣本中都是低甲基化狀態(tài)，而看上去很靠譜的是NUAK1基因，其正常樣本在TSS前1.5kb內(nèi)低甲基化，腫瘤樣本中對應區(qū)域高甲基化。
圖表 9 methylation heatmap NUAK1基因TSS前1.5kb內(nèi)共檢測了7個CpG，這7個CpG在154個樣本中檢測出來的甲基化程度如下圖，可以明顯看出來這7個CpG在Tumor組織中甲基化程度都相對高，而在Normal組織中甲基化程度相對較低。
圖表 10 NUAK1的TSS區(qū)CpG甲基化程度 使用Cpgplot預測CpG island位置，這7個CpG在promoter5’到3’序列圖上相對位置如下 1035 1094 1408 1413 1444 1448 1471
圖表 11 CpG island預測 參數(shù)使用： NUAK1promoter from 1 to 1500 Observed/Expected ratio > 0.40 Percent C + Percent G > 40.00 Length > 100 CpG island詳細信息: Length 101 (1086..1186) Length 105 (1366..1470) 這七個CpG基本都在CpGisland中，具體序列見附錄 九、靶標基因相關討論 進入Gene數(shù)據(jù)庫搜索NUAK1相關內(nèi)容，可以發(fā)現(xiàn)基因全稱NUAK family kinase 1，還是個激酶，激酶的話就對調(diào)控會有很大作用了，而在HPA RNA-seq normal tissues項目中，又看出來這個激酶在腦中表達量明顯高于其他組織，這又與發(fā)生在腦部的GBM不謀而合。
圖表 12 NUAK1相關討論 十、分子機制探討 對于腫瘤組織中高甲基化CpG附近的，并且在腫瘤樣本中低表達的intersect共計274個基因，使用Gene Ontology進行富集分析，可以明顯發(fā)現(xiàn)在GO biological process生物學過程中的“神經(jīng)系統(tǒng)發(fā)育”、“化學性突觸傳遞”和“細胞膜的組織”等部分里面有著富集，特別是“中樞神經(jīng)系統(tǒng)的髓鞘形成”，富集程度達到26.95倍，這又與研究的多發(fā)生于腦補的GBM有著密切的聯(lián)系，反向驗證實驗結(jié)果的正確性。
圖表 13 GO富集分析 十一、FurtherMore 根據(jù)生物學知識可以得到，CpG的甲基化會調(diào)控基因的轉(zhuǎn)錄，因此，Transcript Start Site(TSS)附近的甲基化程度值得進行一番深入研究，選用人類基因組hg19版本，對23056基因共計46489個轉(zhuǎn)錄起始位點，進行轉(zhuǎn)錄起始位點富集甲基化程度統(tǒng)計。 統(tǒng)計TSS前后5000bp內(nèi)CpG甲基化程度，并且使用曲線進行擬合，可以發(fā)現(xiàn)TSS處的CpG Methylation水平明顯降低，這也與科學常識相吻合。
圖表 14 TSS附近甲基化程度 |
|
|
