背景

章老师之前送了胚和胚乳的甲基化数据，现在数据返回来了，可以进行分析了。

数据位置

/data3/jiangchenhao/methylation

参考文献

既然是方法性质的工作，那么随便找几篇参考文献来看看。

香榧NC

首先是香榧的Nc，我手里有香榧Nc的甲基化结果，所以可以直接那来看。

香榧Nc文章对应的测序报告，其中的DMR是（differentially methylated regions），其中的

现在我手里的数据对应的步骤相当于是到Clean data

根据报告上说的，从Clean data开始，他使用的甲基化参考基因组比对工具是Bismark 。

得到比对结果之后，采用Bismark（Krueger, 2011）进行甲基化位点检测。

也就是说，是一个软件干了所有的事情，所以这个流程还是不太可靠的。

董老师的脚本，只包含了比对

Plos one 文章

A comparison of DNA methylation detection between HiFi sequencing and whole genome bisulfite sequencing in monozygotic twins with Down syndrome

具体的方法

WGBS data was analyzed to determine methylated CpGs (mCs) using the wg-blimp analysis pipeline v0.9.10 [15]. In brief, with the input of FASTQ files and the reference genome, sequence reads were aligned to the hg38 reference genome using Bwa-Meth [16] and deduplicated with picard [17]. After that, read quality was evaluated by FastQC [18] and Qualimap [19]. Then, methylation calling was performed by MethylDackel [20]. While wg-blimp was selected for its comprehensive and reliable WGBS workflow, Bismark v0.24.2 [21] was also used to validate the results and rule out tool-specific biases. In this approach, reads were aligned to a bisulfite-converted genome using Bismark with default settings, followed by deduplication and methylation extraction. Moreover, methylation levels in non-CpG contexts (CHG and CHH), as determined by Bismark [21] and MethylDackel [20], were analyzed to assess broader methylation patterns and bisulfite conversion efficiency. The estimated bisulfite conversion efficiency was calculated a

这个工作没有具体的参数啊

Bib调查性综述

A survey of the approaches for identifying differential methylation using bisulfite sequencing data

这个文章发在2018年

这个文章好，用来了解甲基化的来龙去脉，是很不错的。

摘要

DNA甲基化是一种重要的表观遗传机制，在细胞调控系统中发挥着关键作用。近年来测序技术的进步使我们能够生成高通量甲基化数据，并实现单碱基分辨率的甲基化测量。然而，海量数据也带来了挑战，其中DNA甲基化研究的关键挑战之一在于识别不同生物条件下碱基对甲基化水平的显著差异。目前已有多种计算方法被开发用于基于亚硫酸氢盐测序数据识别差异甲基化，但现有方法之间尚未达成明确共识。对这些方法进行全面梳理将极大帮助潜在用户和研究者全面了解现有资源。本文详细综述了22种此类方法，重点关注其基础统计模型、主要特征、核心优势和重要局限性。值得注意的是，本综述指出的这些方法固有缺陷可能为未来研究提供改进方向。

引言

表观遗传学是研究活细胞内基因如何、在何处以及何时被激活或抑制的科学领域。DNA甲基化作为该领域研究最深入的表观遗传机制之一，在多种生命进程中发挥着关键作用[1]。由于其在基因表达调控中的重要性，DNA甲基化参与细胞发育、分化等基本生物学过程。此外，人类癌症中已发现异常高甲基化模式，这为研究此类复杂疾病的发生发展提供了新视角[2]。特别是在癌症中，抑癌基因沉默的一个重要原因就与高甲基化有关。

目前研究最充分的DNA甲基化形式是5-甲基胞嘧啶(5-mC)，即在DNA链胞嘧啶(C)碱基的5号碳原子上添加甲基基团。虽然人类基因组中仅有约5%的胞嘧啶会发生甲基化，但CG二核苷酸(CpG位点)的甲基化频率高达70-80%[3,4]。非CpG环境(如CHG和CHH位点，其中H代表C、T或A)也可发生甲基化，这种现象在植物和干细胞中尤为常见[5,6]。最新研究发现TET蛋白能将5-mC氧化为5-羟甲基胞嘧啶(5-hmC)、5-甲酰基胞嘧啶(5-fC)和5-羧基胞嘧啶(5-caC)，但这些修饰变体的丰度远低于5-mC[7]。鉴于大多数分析方法都针对CpG环境下的5-mC设计，本次研究也聚焦于此类型表观修饰。

当启动子区域的CpG位点发生甲基化时，通常通过阻碍特定转录因子(TF)的结合来抑制该区域的转录活性；反之，非甲基化的启动子CpG位点则允许转录因子结合[8-10]。鉴于DNA甲基化在细胞活动中的重要调控作用，鉴定不同生物条件下DNA甲基化的变化具有重大研究价值。

随着参考基因组和先进测序技术的普及，现已开发出能在全基因组范围内提供高分辨率甲基化谱的分析方法。根据甲基化水平的检测分辨率，当前基于测序的技术可分为两类：（1）基于富集的方法；（2）基于亚硫酸氢盐测序的方法[11,12]。前者可实现100-200碱基分辨率的甲基化水平检测，而后者则能达到单碱基分辨率。

全基因组甲基化检测面临的主要挑战之一是对生物样本量的需求，这一要求直到最近才达到适合临床样本的水平[13]。其他挑战涉及处理新兴技术产生的数据，并将其与其他数据类型（如甲基化与基因表达）进行有意义的整合以获取生物学洞见。本综述将重点关注基于亚硫酸氢盐测序的分析方法。

亚硫酸氢盐测序

检测胞嘧啶甲基化的金标准是亚硫酸氢盐测序技术，其优势在于可实现单碱基分辨率的甲基化检测。该技术通过亚硫酸钠处理DNA，使未甲基化胞嘧啶(C)脱氨基转变为尿嘧啶(U)，而甲基化胞嘧啶保持不变。在测序过程中，尿嘧啶会被识别为胸腺嘧啶(T)。通过计算C/(C+T)比值即可估算每个CpG位点的甲基化水平，从而实现甲基化与非甲基化CpG位点的序列特异性鉴别[18]。

基于亚硫酸盐转化的甲基化检测技术已发展出多种方案。其中全基因组亚硫酸氢盐测序(WGBS)可提供最全面的全基因组DNA甲基化谱，但对大基因组生物进行全基因组测序成本较高。限制性酶切位点亚硫酸氢盐测序(RRBS)及其增强版等更具成本效益的方案，通过靶向富集符合特定长度要求的CpG富集区域，显著降低了测序需求[19]，因此更适用于需要多重重复的研究。

图2展示了亚硫酸氢盐测序数据分析的整体流程，包含六大核心环节：(1)输入数据(FASTA/FASTQ格式的甲基化数据和参考基因组)；(2)数据处理与质量控制；(3)短读段与参考基因组比对；(4)比对后分析；(5)差异甲基化分析；(6)输出结果(差异甲基化位点DMCs、差异甲基化区域DMRs及甲基化模式)。各环节的具体内容将在后续章节详述。

预-分析

数据预处理

到我们手里都是处理之后的clean data

读段比对

经过质量控制后，亚硫酸氢盐测序读段可以与参考基因组进行比对以估计甲基化水平。直接使用标准比对工具会导致比对效率低下，因为亚硫酸氢盐处理将未甲基化的胞嘧啶转化为胸腺嘧啶，导致测序读段与参考基因组之间存在额外的差异。因此，研究人员提出了新的亚硫酸氢盐测序读段比对策略。

现有的亚硫酸氢盐测序比对方法可分为两类：三字母比对器和通配符比对器。三字母比对器（如Bismark[26]、BS Seeker[27]、MethylCoder[28]、BRAT[29]和GNUMAP-bs[30]）会将参考基因组正链的所有C转换为T，反链的所有G转换为A。然后使用标准基因组比对工具（如Bowtie[31]和Bowtie2[32]）将经过相同转换的读段与这些预转换的参考基因组进行比对。

我们Nc的文章使用的就是Bismark这样的工具

相比之下，通配符比对器（如BSMAP[33]、RRBSMAP[34]、GSNAP[35]和RMAP[36]）会将参考基因组的C替换为通配符Y，使其能与测序读段中的C和T都匹配。比对结果通常以SAM/BAM文件格式存储。

比对后分析

在完成读段比对后，可在差异甲基化分析前进行可选的比对后处理步骤，以从比对结果中提取具有生物学意义的信息。目前已开发出多款比对后分析工具，包括：BiQ Analyzer [37] QUMA [38] BRAT [29] MethyQA [39] BSPAT [40] MethGo [41]。大多数工具可提供甲基化数据的统计摘要、质量评估和可视化功能。部分工具还具备额外特性：读段比对（BSPAT和BRAT）

DNA甲基化共现模式识别（BSPAT），单核苷酸多态性与拷贝数变异检测（MethGo），等位基因特异性甲基化模式分析（BSPAT）

这两个标黄的东西，一个可以用甲基化数据来分析胚和胚乳的2n是否完全相同。一个可以用来看甲基化有没有单倍型特异性。

我们自己使用的流程

数据预处理

虽然上面已经说了不需要数据预处理这一步，但是还是有必要把两组测序reads给合并一下，他当时给了我们两组测序reads，据说是因为数据量不够。

奇怪，我明明用nohup挂后台了，但是这个命令还是挂掉了，那么重新跑一遍。

nohup zcat B20903_1.clean.paired.R2.fq.gz  B20903_2.clean.paired.R2.fq.gz > B20903_merged_R2.fq.gz &
nohup zcat B20903_1.clean.paired.R1.fq.gz B20903_2.clean.paired.R1.fq.gz > B20903_merged_R1.fq.gz

读段比对

Bismark安装

直接使用mamba从bioconda安装就可以了，很简单

conda install -c bioconda bismark

尽管这个程序有minimap2等依赖，但是bismark都会提前配置好的。

Bismark软件使用说明

https://felixkrueger.github.io/Bismark/quick_reference/

支持的BS-Seq文件类型

1. 序列格式：FastQ 或 FastA
2. 单端（single-end）或双端（paired-end）读段
3. 输入文件可为未压缩或 gzip 压缩（扩展名 .gz）

我们的数据

1. 支持可变读长
2. 定向（directional）或非定向（non-directional）BS-Seq 文库

比对模式参数

1. 可调整种子长度（adjustable seed length）
2. 允许错配数（number of mismatches）
3. 插入片段大小（insert size）

甲基化统计报告

1. 分析的序列总数（Number of sequences analysed）
2. 唯一最佳比对的序列数（unique best alignment）
3. 甲基化胞嘧啶数（methylated cytosines）
4. 非甲基化胞嘧啶数（unmethylated cytosines）
5. 甲基化百分比计算：
   • CpG、CHG、CHH 上下文（context）
   • 计算公式：

快速参考

Bismark需要安装Perl并在命令行运行。此外，您的计算机还需安装Bowtie 2或HISAT2。如需了解如何搭配Bowtie 2运行Bismark，请参阅本手册末尾部分。

从0.14.0或更高版本开始，Bismark支持并行化处理比对和甲基化提取步骤。具体参数请查阅下方--parallel/--multicore的说明。

准备工作

1. 下载参考基因组：从Ensembl或NCBI等网站获取（示例使用Homo sapiens基因组），保存至指定文件夹。
   • 格式要求：FastA格式（扩展名为.fa或.fasta），支持单条目或多条目文件。

示例数据

示例文件test_dataset.fastq可从Bismark项目或Github页面下载（包含10,000条读段）：

• 格式：FastQ（Phred33质量值）
• 读长：50 bp
• 文库类型：人源定向BS-Seq文库

基因组准备

用法：

bismark_genome_preparation [选项] <基因组文件夹路径>

典型基因组索引命令示例：

/bismark/bismark_genome_preparation --path_to_aligner /usr/bin/bowtie2/ --verbose /data/genomes/homo_sapiens/GRCh37/

明白了，这个东西也是和bwa一样需要构建索引的
我们构建好索引的脚本的路径在/data3/92chuand/methylation/ref_genome/

序列比对

用法：

bismark [选项] --genome <基因组文件夹> {-1 <配对文件1> -2 <配对文件2> | <单端文件>}

典型比对示例：

bismark --genome /data/genomes/homo_sapiens/GRCh37/ test_dataset.fastq

将生成两个输出文件：

• test_dataset_bismark_bt2.bam（包含所有比对结果及甲基化调用字符串）
• test_dataset_bismark_SE_report.txt（包含比对和甲基化摘要）

注意：

当前工作目录必须包含待分析的序列文件。

去重

deduplicate_bismark --bam [选项] <文件名>

此命令将对Bismark比对的BAM文件进行去重，仅保留相同位置和相同方向的一条读段。建议用于全基因组亚硫酸盐测序样本，但不适用于RRBS、扩增子或靶向富集文库。

甲基化提取

用法：

bismark_methylation_extractor [选项] <文件名>

典型上下文依赖（CpG/CHG/CHH）甲基化提取命令示例：

bismark_methylation_extractor --gzip --bedGraph test_dataset_bismark_bt2.bam

将生成三个甲基化输出文件：

• CpG_context_test_dataset_bismark_bt2.txt.gz
• CHG_context_test_dataset_bismark_bt2.txt.gz
• CHH_context_test_dataset_bismark_bt2.txt.gz

以及bedGraph和Bismark覆盖文件。

样本报告

用法：

bismark2report [选项]

此命令尝试查找Bismark比对、去重和甲基化提取报告以及M-bias文件，为目录中的每个样本生成图形化HTML报告（示例见Bismark双端报告）。

汇总报告

用法：

bismark2summary [选项]

此命令扫描当前工作目录中的Bismark比对、去重和甲基化提取报告，为目录中的所有文件生成图形化汇总HTML报告和数据表格（示例见Bismark汇总报告）。适用于快速可视化查看大量样本（数十、数百或数千个）的比对统计信息，单个样本报告请使用bismark2report。