第一章：从Raw Data到发表级图表——甲基化分析全景概览
DNA甲基化是表观遗传调控的核心机制之一，广泛参与基因表达沉默、X染色体失活及肿瘤发生等生物学过程。高通量测序技术的发展使得全基因组甲基化分析成为可能，典型流程涵盖原始数据获取、质量控制、比对、甲基化位点识别至可视化呈现。
数据预处理与质控
原始测序数据（如Bisulfite-Seq或WGBS）通常以FASTQ格式存储。首先需进行适配子剪切和低质量过滤：
# 使用Trimmomatic进行数据清洗 java -jar trimmomatic.jar PE -phred33 \ sample_R1.fq.gz sample_R2.fq.gz \ cleaned_R1.fq cleaned_R2.fq \ ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 \ SLIDINGWINDOW:4:20 MINLEN:50
该命令去除接头序列，并采用滑动窗口策略过滤质量低于Q20的片段。
比对与甲基化 Calling
经质控后的数据需比对至参考基因组（如hg38），常用工具包括Bismark或BS-Seeker2。Bismark基于Bowtie2实现亚硫酸盐转化序列比对：
bismark --genome_folder hg38_index/ cleaned_R1.fq
随后运行甲基化提取模块：
bismark_methylation_extractor *.bam
输出每个CpG位点的甲基化水平（如chr1:10000, 5/10 reads methylated）。
结果可视化
高质量图表是成果发表的关键。可使用R语言ggplot2绘制甲基化密度图或热图：
- 加载甲基化β值矩阵
- 聚类样本并构建层次聚类热图
- 标注临床分组信息以揭示模式差异
| 样本 | 平均甲基化水平 (%) | 检测CpG数 |
|---|
| Tumor_01 | 78.3 | 2,845,102 |
| Normal_01 | 62.1 | 2,790,441 |
graph LR A[Raw FASTQ] --> B[Quality Control] B --> C[Alignment] C --> D[Methylation Calling] D --> E[Data Visualization]
第二章：GEO数据库甲基化数据挖掘与获取
2.1 甲基化芯片技术原理与GEO数据结构解析
甲基化芯片技术原理
甲基化芯片通过探针特异性结合基因组CpG位点，检测DNA甲基化水平。以Illumina Infinium平台为例，利用亚硫酸盐处理DNA，将未甲基化的胞嘧啶转化为尿嘧啶，而甲基化的胞嘧啶保持不变，随后通过荧光信号强度计算β值（β = M / (M + U + α)），反映甲基化程度。
# 示例：计算甲基化β值 beta_value <- methylated_intensity / (methylated_intensity + unmethylated_intensity + 100)
该公式中，分子为甲基化探针信号，分母包含非甲基化信号及常数α用于稳定低强度噪声。
GEO数据结构解析
GEO数据库中甲基化数据通常以Matrix格式存储，包含GPL（平台）、GSM（样本）、GSE（系列）三类记录。常见文件如
series_matrix.txt提供元数据，而原始信号值则分布于RAW文件中。
| 字段 | 含义 |
|---|
| GPL | 芯片平台信息 |
| GSM | 单个样本表达谱 |
| GSE | 实验系列集合 |
2.2 基于R的GEO数据批量下载与元信息提取
数据获取流程
使用
GEOquery包可实现GEO数据库中高通量表达数据的批量下载。通过指定多个GSE编号，自动化获取原始表达矩阵与样本注释信息。
library(GEOquery) gse_ids <- c("GSE12345", "GSE67890") gse_list <- lapply(gse_ids, getGEO)
上述代码遍历GSE编号列表，调用
getGEO()函数从NCBI服务器下载对应数据集，返回ExpressionSet或ESETList对象，包含表达值、探针注释和样本元数据。
元信息解析
提取样本表型信息是后续分析的关键步骤。可通过
pData()函数获取每个数据集的表型数据框。
- 样本分组（如疾病 vs 正常）
- 临床特征（年龄、性别等）
- 实验平台编号（GPLXX）
结合
meta()函数可进一步提取数据集级元信息，如研究标题、作者与DOI链接，便于构建可追溯的数据目录。
2.3 IDAT文件与表达矩阵的获取策略
数据来源与IDAT文件结构
Illumina甲基化芯片生成的IDAT文件存储了每个探针的荧光信号强度，分为红绿双通道。这些二进制文件需通过背景校正与归一化处理，转化为可用的甲基化β值。
表达矩阵构建流程
使用
minfi包读取IDAT文件并构建
RGSet对象：
library(minfi) targets <- read.metharray.sheet("sample_sheet.csv") rgSet <- read.metharray.exp(targets = targets)
该代码段解析样本表并批量导入IDAT数据。
read.metharray.exp自动匹配文件名与元数据，生成包含M和U信号的原始信号集。
质量控制与矩阵转换
经背景校正后，采用SWAN或Functional Normalization方法进行批次效应校正，最终通过
betaEstimate()函数输出标准化后的表达矩阵，供下游差异分析使用。
2.4 数据质量控制与样本分组设计
在构建可靠的数据分析模型前，必须确保输入数据的准确性与一致性。数据质量控制涵盖缺失值处理、异常检测和重复记录清洗等关键步骤。
常见数据清洗策略
- 使用均值或中位数填补数值型字段的空缺
- 通过IQR方法识别并剔除离群点
- 基于唯一键去重，避免样本偏差
样本分组设计原则
为保证实验结果的可解释性，常采用随机分层抽样。以下为Python示例代码：
from sklearn.model_selection import train_test_split X_train, X_test, y_train, y_test = train_test_split( X, y, stratify=y, test_size=0.2, random_state=42 )
该代码实现按标签分布分层划分训练集与测试集，参数`stratify=y`确保各类别比例一致，`test_size=0.2`表示测试集占比20%，`random_state`保障结果可复现。
2.5 表型数据整理与临床信息关联
数据清洗与标准化
表型数据常来源于电子病历、问卷调查和临床试验，存在缺失值、单位不统一等问题。需进行去重、归一化和编码转换。例如，将性别字段“男/女”映射为“M/F”，使用以下代码实现：
import pandas as pd # 加载原始表型数据 phenotype_df = pd.read_csv("phenotype_raw.csv") # 标准化性别字段 phenotype_df['gender'] = phenotype_df['gender'].replace({'男': 'M', '女': 'F'}) # 填充年龄缺失值为中位数 phenotype_df['age'].fillna(phenotype_df['age'].median(), inplace=True)
上述代码首先读取原始数据，随后对分类变量进行一致性编码，并对数值型缺失字段采用中位数填充策略，提升数据完整性。
临床信息整合
通过唯一患者ID，将表型数据与基因组数据、影像学报告等临床信息进行横向关联，构建统一分析数据集。常用方式如下：
| PatientID | Age | Gender | Diagnosis | Genetic_Marker |
|---|
| PT001 | 45 | M | Diabetes | SNP_X_12345 |
| PT002 | 52 | F | Hypertension | SNP_Y_67890 |
该表格展示了整合后的多维数据结构，支持后续的关联分析与机器学习建模。
第三章：甲基化数据预处理与标准化
3.1 使用minfi进行原始信号读取与背景校正
在Illumina甲基化芯片数据分析流程中，`minfi`是R语言中广泛使用的工具包，专门用于处理IDAT格式的原始信号文件。它能够高效读取探针级强度数据，并提供多种背景校正方法以提升数据质量。
加载与读取IDAT文件
library(minfi) baseDir <- "path/to/idat/files" targets <- read.metharray.sheet(baseDir) rgSet <- read.metharray.exp(baseDir = baseDir)
该代码段首先读取样本信息表（targets），然后通过
read.metharray.exp()函数批量导入所有IDAT文件，生成包含红绿通道信号的
RGChannelSet对象。
背景校正策略
minfi支持多种校正算法，如"noob"（normal-exponential out-of-band）：
bg_corrected <- preprocessNoob(rgSet, fixOutliers = TRUE)
preprocessNoob()函数对每个探针执行基于负对照通道的背景扣除，有效降低技术噪声，尤其适用于低信号强度探针。设置
fixOutliers = TRUE可进一步修正异常值，提升数据稳定性。
3.2 探针过滤与性别相关位点去除
在基因型数据预处理中，探针过滤是确保数据质量的关键步骤。需首先排除检测失败率高或变异类型异常的探针，避免对后续分析造成偏差。
常见过滤标准
- 检出率低：移除在超过5%样本中未能检出的探针
- 非多态性位点：过滤固定为单一等位基因的SNP
- 性别染色体相关位点：常引入性别偏差，需特殊处理
性别相关位点去除策略
为避免性别染色体对常染色体分析的干扰，通常从参考数据库（如dbSNP）中提取位于X、Y染色体的SNP列表，并从主数据集中剔除。
# 使用R语言进行探针过滤示例 filtered_probes <- probes[probes$chromosome != "X" & probes$chromosome != "Y" & probes$detection_pval < 0.05, ]
上述代码保留染色体非X/Y且检测P值小于0.05的探针，有效提升数据一致性。参数
detection_pval反映探针信号可靠性，阈值通常设为0.05。
3.3 β值计算与M值转换：生物学意义与数学基础
在DNA甲基化研究中，β值和M值是量化CpG位点甲基化水平的核心指标。β值表示甲基化信号占总信号的比例，其取值范围为[0,1]，直观反映生物学意义上的甲基化程度。
β值的数学定义
β = \frac{M}{M + U + \alpha}，其中M为甲基化荧光信号强度，U为非甲基化信号强度，α为平滑常数（通常设为100），用于避免低信号时的噪声干扰。
M值的转换逻辑
M值即对数比值：M = log₂\left(\frac{M + \alpha}{U + \alpha}\right)，具备更好的统计分布特性，适用于差异分析。
beta <- function(methylated, unmethylated, alpha = 100) { (methylated + alpha) / (methylated + unmethylated + 2 * alpha) }
该函数实现β值计算，输入为M和U信号向量，输出为对应β值。加入α可稳定低表达CpG位点的估计。
- β值便于解释：接近1表示高度甲基化
- M值更适合建模：近似正态分布，利于回归分析
第四章：差异甲基化分析与功能注释
4.1 DMP与DMR检测：统计模型选择与R实现
在DNA甲基化分析中，差异甲基化位点（DMP）和差异甲基化区域（DMR）的识别依赖于合适的统计模型。常用方法包括基于β值的t检验、Wilcoxon秩和检验，以及更复杂的广义线性模型（GLM）。
常用R包与函数选择
limma：适用于标准化后的β值矩阵，支持线性模型拟合；missMethyl：专为Illumina甲基化芯片设计，可校正CpG位点间的相关性；DMRcate：基于limma输出结果，合并邻近DMP形成DMR。
# 使用limma进行DMP检测 library(limma) design <- model.matrix(~0 + group) fit <- lmFit(beta_matrix, design) fit <- eBayes(fit) dmp_results <- topTable(fit, coef = 2, number = Inf)
上述代码构建分组设计矩阵，通过经验贝叶斯平滑估计差异显著性。其中
coef = 2指定比较第二组别，
eBayes提升小样本下的稳定性。后续可将
dmp_results输入至
DMRcate进行区域聚合。
4.2 差异结果可视化：热图、火山图与箱线图绘制
热图展示差异表达模式
热图适用于呈现多个样本间基因或蛋白的表达趋势。使用R语言
pheatmap包可快速生成高质量热图：
library(pheatmap) pheatmap(log2(expr_matrix + 1), scale = "row", clustering_distance_rows = "correlation", show_rownames = FALSE)
其中，
scale = "row"对每行进行标准化，突出表达模式；
clustering_distance_rows调整聚类距离算法，提升分组逻辑性。
火山图识别显著差异分子
火山图结合统计显著性与变化倍数，直观筛选关键分子。常用参数包括log2FoldChange和-log10(padj)。
- 横轴表示变化倍数（log2FC）
- 纵轴表示显著性（-log10(padj)）
- 显著上调/下调点以不同颜色标注
4.3 基因组注释与CpG岛富集分析
基因组注释基础
基因组注释是识别DNA序列中功能元件的过程，包括基因、启动子及非编码RNA等。常用工具如GENCODE和Ensembl提供高质量的参考注释文件（GTF格式），可用于后续分析。
CpG岛检测与富集分析
CpG岛是富含CG二核苷酸的区域，常位于基因启动子区，与基因表达调控密切相关。通过bedtools可进行CpG岛与基因组功能区域的重叠分析：
# 计算CpG岛与启动子区域的交集 bedtools intersect -a cpg_islands.bed -b promoters.bed -wa -wb > cpg_in_promoters.bed
该命令输出位于启动子区内的CpG岛记录，用于评估其在转录调控中的潜在作用。
富集结果可视化
使用表格汇总富集结果，便于比较不同区域的CpG分布：
| 基因组区域 | CpG岛数量 | 占比(%) |
|---|
| 启动子区 | 1250 | 45.2 |
| 基因内区 | 980 | 35.5 |
| 基因间区 | 530 | 19.3 |
4.4 功能富集分析与通路关联（GO/KEGG/GSVA）
功能富集分析概述
功能富集分析用于识别差异表达基因在生物学过程、分子功能和细胞组分中的显著性聚集。常用方法包括GO（Gene Ontology）和KEGG（Kyoto Encyclopedia of Genes and Genomes）通路分析。
- GO分析分为三个维度：BP（生物过程）、MF（分子功能）、CC（细胞组分）
- KEGG则聚焦于基因参与的代谢与信号通路
- GSVA（Gene Set Variation Analysis）实现样本级通路活性评分
gsva_result <- gsva(expr_matrix, gene_sets, method="ssgsea")
该代码使用GSVA对表达矩阵
expr_matrix进行通路活性评分，
gene_sets为基因集列表，
method="ssgsea"指定采用单样本GSEA算法，适用于无分组信息的连续样本分析。
结果可视化示例
| 通路名称 | p值 | FDR |
|---|
| Apoptosis | 0.001 | 0.008 |
| Cell Cycle | 0.0003 | 0.002 |
第五章：迈向高质量科研图表与论文投稿
选择合适的图表类型提升数据表达力
科研图表应根据数据特性选择恰当的可视化形式。例如，时间序列数据适合使用折线图，分类比较推荐柱状图，而相关性分析可采用散点图。在 Python 中，Matplotlib 与 Seaborn 提供了高度可定制的绘图接口：
import seaborn as sns import matplotlib.pyplot as plt # 设置美观主题 sns.set_theme(style="whitegrid") plt.figure(figsize=(8, 6)) # 绘制带置信区间的回归图 sns.regplot(data=df, x='variable_x', y='variable_y', scatter_kws={'alpha':0.6}) plt.xlabel("自变量（单位：mm）") plt.ylabel("因变量（单位：kPa）") plt.title("变量间线性关系及95%置信区间") plt.savefig("regression_plot.pdf", bbox_inches='tight')
遵循期刊图表规范确保顺利投稿
不同期刊对图像分辨率、字体大小和格式有明确要求。常见要求包括：
- 分辨率不低于 300 dpi
- 字体统一使用 Arial 或 Times New Roman
- 图像格式优先选择 PDF、EPS 或 TIFF
- 图注需独立于图像文件提供
LaTeX 环境下的图表嵌入实践
使用 LaTeX 撰写论文时，推荐通过
graphicx宏包插入矢量图，以保证印刷质量：
\begin{figure}[htbp] \centering \includegraphics[width=0.8\textwidth]{regression_plot.pdf} \caption{实验组与对照组的响应趋势对比} \label{fig:response} \end{figure}
投稿前的关键检查清单
| 检查项 | 完成状态 |
|---|
| 图表分辨率达标 | ✅ |
| 坐标轴标签完整 | ✅ |
| 图例清晰无遮挡 | ✅ |
