|
这是一份大三学生的课程实践作业~~~ 找出胶质细胞瘤特异性甲基化区域，为临床诊断提供理论依据
可靠性 一、数据下载 首先进入TCGA下载数据GBM的RNA-seq和甲基化数据，从下表可见GBM共有172套RNA-seq数据以及437套DNA甲基化数据，由于TCGA提供Infinium HumanMethylation27 BeadChip和Infinium HumanMethylation450 BeadChip两种芯片平台的数据，为了避免后续不同芯片平台间数据合并的困难，仅下载HumanMethylation450的芯片数据，共计154套。 图表 1TCGA数据汇总 二、初步整理数据 使用TCGA-Assembler.2.0.5进行GBM数据批量下载与初步整理，并且绘制RNA-seq 基因表达量盒型 图 以及甲基化芯片数据盒型图 ，由于数据量较大，此处不贴图。 三、整体可视化 首先对于甲基化数据，选取ID为TCGA.06.AABW.11A.31D.A368.05的数据，查看总体甲基化程度。由于每个位点真实情况只存在：甲基化/非甲基化两种，所以对全部位点甲基化程度进行统计，也应该是大部分位点处于“完全甲基化”(Methylation state=1)和“完全非甲基化”(Methylation state=0)两种状态，下图绘制了数据的频数柱状图，可以明显看出形状处于“两头高，中间低”，反向说明芯片数据质量较好。 图表 2单个样本CpG甲基化程度统计 接下来，对多个样本绘制CpG甲基化程度小提琴图，同一行是同一个病人，左边样本来源于Primary Solid Tumor，右边样本来源于Recurrent Solid Tumor，除了甲基化程度大部分分布于0和1附近外，还能看出来源于同一病人肿瘤的甲基化程度依旧会有略微差异。 TCGA barcode：https://wiki.nci./display/TCGA/TCGA+barcode 图表 3小提琴图 同样的，对于RNA-seq数据也可以进行一些初步可视化，除了数据下载后绘制的盒型图，亦可以进行PCA初步查看数据分布，下图左为PCA陡坡图，反映了第一主成分、第二主成分…等等所拥有信息量的比例，下图右为使用PCA1和PCA2绘制的散点图，可以发现5个正常样本距离较近，从侧面反映数据可信度较好。 最后，对于RNA-seq表达谱数据，使用系统聚类方法，绘制树状图，可以发现5个正常样本距离也是很近，数据质量还行。 四、差异甲基化区域筛选 为了更加科学高效地筛选差异甲基化位点，参考bioconductor中甲基化芯片的分析流程，使用minfi包进行差异甲基化分析，得到差异甲基化位点。 http://www./help/workflows/methylationArrayAnalysis/ 在检测的526733个CpG位点中，共有4927个CpG位点P值<0.01，且在肿瘤样本中保持着甲基化程度高于0.7，对应2054个基因。 五、差异表达基因筛选 由于数据源自RNA-seq，最主流的分析方法当然是基于负二项分布模型的DESeq2包。 先用MA-plot查看差异表达基因大致分布。意外的是，图形左侧有大概七条线状条纹，最初我怀疑这是sample之间有batch effect导致，需要用其他更好normalize的方法，后来用identify方法挨个找出每条线上的基因名及其对应的表达量，发现这些基因在172套样本中表达量几乎全为0，仅有一两个样本有一点点表达，这种数据的存在导致这些线状条纹的产生。 图表 4 MA-plot 然后， 选取p值最小的差异表达基因，绘制其在不同组间表达量，确实差异很显著。 图表 5表达量散点图 接着，绘制差异表达基因在不同组间的表达量热图，正常样本是图片最左边的五列，当然如果需要解释具体的生物学问题，需要将聚类出来的每一类，将差异表达基因进行GO以及KEGG注释，结合有关的生物学表型，探讨其分子机制及意义 图表 6表达量热图 最后，选取筛选条件为p值小于0.01且log2FoldChange<-2的差异表达基因，在肿瘤样本低表达的基因，共计1657个 六、抑癌基因的获取 在pubmed中查询研究人员整理的tumor suppressor genes，果然在Nucleic Acids Res发表过TSGene数据库，存储了抑癌基因的列表。https://bioinfo./TSGene/download.cgi 下载全部1217个人类抑癌基因的列表。 All the 1217 human tumor suppressor genes with basic annotations
图表 7 TSGene database 七、数据整合 对于甲基化数据中，肿瘤样本高甲基化CpG附近的基因，RNA-seq中肿瘤样本低表达的基因，以及TSGene数据库中下载的抑癌基因列表，三者做overlap，找出特异性的候选靶标，为后续分析做准备，下图为三者overlap的韦恩图。
图表 8数据整合韦恩图 共计找出12个候选靶标基因。 八、靶标筛选 之前筛选选择的单个CpG的差异甲基化，而实际临床检测应用时候，可能需要多个CpG作为对照，因此统计了12个候选靶标基因TSS前1.5kb内所有CpG的甲基化程度，然后绘制热图，可以明显发现，虽然当初用CpG的差异甲基化位点筛出来的基因都是肿瘤样本高甲基化的，可是统计TSS前1.5kb内所有CpG的甲基化程度，这些基因却有很多在所有样本中都是低甲基化状态，而看上去很靠谱的是NUAK1基因，其正常样本在TSS前1.5kb内低甲基化，肿瘤样本中对应区域高甲基化。
图表 9 methylation heatmap NUAK1基因TSS前1.5kb内共检测了7个CpG，这7个CpG在154个样本中检测出来的甲基化程度如下图，可以明显看出来这7个CpG在Tumor组织中甲基化程度都相对高，而在Normal组织中甲基化程度相对较低。
图表 10 NUAK1的TSS区CpG甲基化程度 使用Cpgplot预测CpG island位置，这7个CpG在promoter5’到3’序列图上相对位置如下 1035 1094 1408 1413 1444 1448 1471
图表 11 CpG island预测 参数使用： 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数据库搜索NUAK1相关内容，可以发现基因全称NUAK family kinase 1，还是个激酶，激酶的话就对调控会有很大作用了，而在HPA RNA-seq normal tissues项目中，又看出来这个激酶在脑中表达量明显高于其他组织，这又与发生在脑部的GBM不谋而合。
图表 12 NUAK1相关讨论 十、分子机制探讨 对于肿瘤组织中高甲基化CpG附近的，并且在肿瘤样本中低表达的intersect共计274个基因，使用Gene Ontology进行富集分析，可以明显发现在GO biological process生物学过程中的“神经系统发育”、“化学性突触传递”和“细胞膜的组织”等部分里面有着富集，特别是“中枢神经系统的髓鞘形成”，富集程度达到26.95倍，这又与研究的多发生于脑补的GBM有着密切的联系，反向验证实验结果的正确性。
图表 13 GO富集分析 十一、FurtherMore 根据生物学知识可以得到，CpG的甲基化会调控基因的转录，因此，Transcript Start Site(TSS)附近的甲基化程度值得进行一番深入研究，选用人类基因组hg19版本，对23056基因共计46489个转录起始位点，进行转录起始位点富集甲基化程度统计。 统计TSS前后5000bp内CpG甲基化程度，并且使用曲线进行拟合，可以发现TSS处的CpG Methylation水平明显降低，这也与科学常识相吻合。
图表 14 TSS附近甲基化程度 |
|
|
