前面我们系统性的总结了circRNA的相关背景知识:
同样的策略,我们也可以应用到其它领域的知识背景快速学习,比如我们的lncRNA系列,miRNA系列,现在我们一起学习一下甲基化吧。
我们已经强调过主要是DNA的甲基化,而且前面我们介绍了一些背景知识,主要是理解什么是DNA甲基化,为什么要检测它,以及芯片和测序两个方向的DNA甲基化检测技术。具体介绍在:甲基化的一些基础知识
提到这些甲基化芯片背景知识,不得不强推中科院毕业的生信博士joshua
,目前在英国从事博士后研究,他是 ChAMP
这个分析(一站式)450K甲基化芯片数据的R包作者,我们生信技能树也多次推送过相关教程。比如850K甲基化芯片数据的分析
甲基化芯片主要是450K和850K,都是采用了两种探针Infinium Ⅰ 和Infinium Ⅱ对甲基化进行测定,Infinium I采用了两种bead(甲基化M和非甲基化U),而II只有一种bead(即甲基化和非甲基化在一起),这也导致了它们在后续荧光探测的不同,450K采用了两种荧光探测信号(红光和绿光)。数量分布如下:
illumina的软件(Illumina GenomeStudio软件)自动就把450k芯片的IDAT原始文件
转换成了甲基化信号文件,一般是3~8列值每个样本。通常是直接使用Illumina的Bead Studio软件和甲基化模块v3.2计算每个CpG靶标的β值。
β值=来自甲基化珠粒类型的强度值/(来自甲基化的强度值+来自未甲基化珠粒类型的强度值+ 100)。
然后使用p值对数据进行质量过滤。
p值大于0.01的β值被认为低于最小强度,阈值显示为“NA”。
有一个2010文章介绍了使用Beta value 和 M-value进行统计分析的利弊,大家可以参考一下。Comparison of Beta-value and M-value methods for quantifying methylation levels by microarray analysis
而minfi
包有getM
和getBeta
来分别计算M-values和Beta-values,包的作者认为
minfi包的一个函数read.450k.exp
也可以直接读取.IDAT文件
公式计算:平均β=信号B /(信号A +信号B + 100)
通过计算甲基化(信号A)和未甲基化(信号B)等位基因之间的强度比来确定DNA甲基化水平(β值)。
具体地,β值是由甲基化(M对应于信号A)和未甲基化(U对应于信号B)等位基因的强度计算的,荧光信号的比率β= Max(M,0)/ [Max( M,0)+ Max(U,0)+ 100]
。
因此,β值的范围从0(完全未甲基化)到1(完全甲基化)
一般来说,具体的β值
的意义是:
甲基化分析需要对应的注释文件,主流是EPIC和450K的,我先分析450K的Manifest
,首先原有的Manifest包含了从BeadChip到最终的文件的对应号,但是有一部分信息应该要提前过滤掉:一部分是开头的Header,另一部分是结尾的Control Probe,可以直接从illumina官网下载到对应的450K注释文件,如果把Header,Control Probe和SNP全部删掉,450K数据的行数正好就是:485512。
这就是450K甲基化注释的所有Probe数量,每一个Probe对应一个CpG位点。不是说人体的全基因组上只有这些位点,而是说illumina公司决定只将这些位点设计到芯片中完成测序。设计的时候,参考了很多公共数据库:
•Genomic Coordinates
•UCSC RefGene Name
•UCSC RefGene Accession
•UCSC RefGene Group
•UCSC CpG Islands Name
•Relation to UCSC CpG Island (Island, Shore, Shelf)
•Phantom
•DMR
•Enhancer
•HMM Island
•Regulatory Feature Name
•Regulatory Feature Group
而且位于基因组不同功能区域的甲基化探针信号值代表的生物学意义不一样,包括5’ UTR, first exon, gene body, 3’ UTR, CpG island, CpG shore, CpG shelf,这些注释,都是在厂商提供注释文件信息里面。
450K和850K芯片的甲基化探针数量庞大,而我们人类的蛋白编码基因也就两万个,所以每个基因评价可以有几十个探针,而我们比较关心的是该基因的启动子区域的甲基化探针的信号值水平。面临两个问题:
Using the median
of the methylation values of the CpG probes mapped at promoter regions
(e.g., 1500 nucleotides from the transcription start site TSS1500; 200 nucleotides from the transcription start site TSS200; 5’UTR; and 1st exon), a single methylation value was consolidated for each gene promoter represented in the HM450K platform.
参考文章:https://link.springer.com/article/10.1186/s41241-017-0041-9
关于GEOquery在GEO公共数据库挖掘的重要性我已经多次强调和介绍了,不赘述,不过它可不只是用来下载表达芯片矩阵的,也可以下载甲基化芯片信号矩阵,也就是贝塔值矩阵。
代码如下:
require(GEOquery)
require(Biobase)
GSE80559 <- getGEO("GSE80559")
beta.m <- exprs(GSE80559[[1]])
值得一提的是,甲基化信号矩阵有成熟的R包来区分细胞类型。
比如8种血液细胞的亚型 (B-cells, CD4+ T-cells, CD8+ T-cells, NK-cells, Monocytes, Neutrophils, Eosinophils, and Granulocytes. 你也可以看看单细胞水平的甲基化信号值,是否可以根据传统bulk的细胞marker来进行区分。
R包是:EpiDISH
大部分情况下的甲基化芯片数据分析,与mRNA表达芯片拿到的表达矩阵下游分析并没有本质区别。甲基化芯片数据分析的输入是每个探针在每个样本的甲基化信号值矩阵,而mRNA表达芯片拿到的表达矩阵是每个探针代表的基因在每个样本的表达量矩阵。
如果从workflow来看,主要是:
下载
,主要是GEO和TCGA数据库差异
甲基化,3个层次的差异注释分析
相关性分析
你会发现,拿到甲基化信号值矩阵后,仍然是可以走标准分析流程,火山图,热图,GO/KEGG数据库注释等等。这些流程的视频教程都在B站和GitHub了,目录如下:
感兴趣可以细读表达芯片的公共数据库挖掘系列推文 ;
所以实际上,还是看大家的R语言水平啦。以及大家的文献阅读功底啦。
比如表达矩阵可以进行分子分型,那么甲基化信号值矩阵也是如此。
比如表达量可以把病人分组做生存分析,那么甲基化信号值也是一样。这个属于TCGA数据挖掘系列文章了,大家搜索关键词, TCGA methylation signature ,会发现大量2019年附近的文章,具体原因我这里不方便细说。
比如基因型可以进行gwas研究,那么甲基化信号也可以。For example, a genome-wide DNA methylation association study in obesity that recruited 5,387 individuals, identified 278 CpGs associated with BMI (Wahl et al., 2017).
大部分数据分析方法都是借鉴成熟的mRNA表达芯片的表达矩阵来的,并无特殊之处,需要注意的仅仅是数据质控部分,以及后续的基于甲基化这个特性的分析。
这个步骤主要是靠多练习,数据集见得多了,就会有自己的判断力。当然了,bioconductor的几个集成好的workflow里面都是着重强调了甲基化数据质控、归一化。尤其是要多花一点时间在上面。
甲基化信号值矩阵有3个层次的差异分析:DMP,DMR,DMB:
简单来说,DMP是找出一个一个的差异甲基化CpG位点,DMR就是一个连续不断都比较长的差异片段,科学家们觉得,这样的连续差异片段,对于基因的影响会更加明显,只找这样的片段,可以使得计算生物学的打击精度更为准确,也可以让最终找出来的结论数据更少,便于实验人员筛选。另外一个类似的东西就是DMB,那个东西出现的原因是,有的科学家觉得,DMR这样的区域还不够显著,DNA上的甲基化出现变化,可能是绵延几千位点的!而且只会在基因以外的区域,但是这些基因以外的区域发生变化,却会导致基因的表达发生变化。你可以想象成,北京周边的河北在大炼钢铁,然后北京也跟着雾霾了,大概就是这意思。
差异甲基化位点的基因功能区域分类
前面我们提到了位于基因组不同功能区域的甲基化探针信号值代表的生物学意义不一样,包括5’ UTR, first exon, gene body, 3’ UTR, CpG island, CpG shore, CpG shelf,这些注释,都是在厂商提供注释文件信息里面。所以差异分析拿到的高甲基化或者低甲基化位点就可以分类,比如文章: Epigenetic regulation of diacylglycerol kinase alpha promotes radiation-induced fibrosis:
还可以更细致的展示某个基因附近的全部甲基化探针,如下:
(e) Detailed interrogation of DNA methylation using EpiTYPER in two distinct patient fibroblasts (fibrosis and fibrosis-free) at the DGKA 5′ UTR (upper panel) and the differentially methylated region (lower panel).
其实我们生信技能树前面的多个教程也系统性介绍过,比如850K甲基化芯片数据的分析 ,还有Bioconductor的DNA甲基化芯片分析流程
可以很清楚的看到,也是以差异分析为主:
甲基化特有的分析。
你可以下载GSE58651数据集,然后看看能不能重复出来同样的差异分析结果,并且绘制文章 里面的差异甲基化信号值位点热图:
Screening of the 201 significant CpGs for gender-specific probes, SNP probes, repeats and unidentified locations, resulted in 81 CpGs (54 hypermethylated and 27 hypomethylated) located in 64 genes.
再比如发表在 January 25, 2018 的文章,Identification of methylated genes and miRNA signatures in nasopharyngeal carcinoma by bioinformatics analysis ,集中在两个 (GSE52068 and GSE62336) 甲基化数据差异分析结果的overlap上面。也是值得大家借鉴!
另外一个简单数据挖掘例子,是GSE62336,文章在:https://www.spandidos-publications.com/10.3892/ol.2017.7083
大部分人感兴趣的是数据挖掘,通常我们说甲基化芯片数据挖掘,不仅仅是前面提到的单一数据的分析,往往是会把它跟mRNA表达芯片联合起来。因为它们可以有共同的基因做桥梁,甲基化探针设计在基因的不同功能区域,包括启动子,外显子,内含子,mRNA表达芯片探针通常是设计在基因的转录本上面。
比如发表在Front. Genet., 19 December 2018 | https://doi.org/10.3389/fgene.2018.00663 的文章,题目很长Using Integrative Analysis of DNA Methylation and Gene Expression Data in Multiple Tissue Types to Prioritize Candidate Genes for Drug Development in Obesit ,主要是关于Visceral adipose tissue (VAT) and subcutaneous adipose tissue (SAT)的分析。
发表在 September 27, 2017 文章 https://www.spandidos-publications.com/ol/14/6/7171 题目是Identification of potentially critical differentially methylated genes in nasopharyngeal carcinoma: A comprehensive analysis of methylation profiling and gene expression profiling
发表在 September 2019 :European Archives of Oto-Rhino-Laryngology 的 文章 ,Comprehensive analysis of gene expression and DNA methylation for human nasopharyngeal carcinoma 也做了挖掘:https://link.springer.com/article/10.1007/s00405-019-05525-2
实际上,这个时候你应该是意识到,之所以联合甲基化芯片和mRNA表达芯片,就是因为它们可以通过基因这个桥梁,所以理论上mRNA-seq当然也可以跟甲基化芯片或者测序联合起来分析的啦!发表在BMC Systems Biology December 2011的文章 An integrative analysis of DNA methylation and RNA-Seq data for human heart, kidney and liver 就是这样的一个例子。
Agilent 的SurePrint甲基化芯片:
Roche 的 NimbleGen甲基化芯片
感兴趣的朋友可以去GEO的GPL合辑里面搜索找一下,看看这些平台到底小众到什么程度。