全基因组亚硫酸盐测序分析实战：从原始数据到甲基化可视化

本教程适用于课堂上的 Arabidopsis thaliana TAIR10 chr1 教学数据。

数据类型：single-end WGBS
参考文件目录：/public/home/bio2023/shiqc/dataforwgbs/ref
原始 reads 目录：/public/home/bio2023/shiqc/dataforwgbs/rawdata-0.1
环境名：wgbs

课前信息

一、课堂数据来源与参考

WGBS / BS-seq / methylC-seq 在拟南芥中的经典开创性工作，可参考 Lister et al. 2008, Cell：Highly Integrated Single-Base Resolution Maps of the Epigenome in Arabidopsis。
本课堂实际使用的两个样本 SRR534177_colWT 和 SRR534239_met1 来自 Stroud et al. 2013, Cell：Comprehensive analysis of silencing mutants reveals complex regulation of the Arabidopsis methylome。
这两个 SRA 运行在 NCBI 中都标记为 Layout: SINGLE，因此本教程按单端 single-end 数据处理。
本教程的课堂化流程和部分命令写法，参考了 SAGC 的 WGBS 教学教程。
课堂提供的原始 reads 使用的是原始公开数据的 10% 子集，用于缩短课堂运行时间。

mindmap
  root((WGBS 教程总览))
    服务器准备
      登录主服务器
      进入计算节点
      创建个人目录
      激活 wgbs 环境
    数据准备
      检查 fa 和 gff3
      复制 raw reads
      生成 genes.bed6
    原始数据处理
      FastQC
      Trim Galore
    比对与甲基化提取
      Bismark 建索引
      Bismark 比对
      去重复和排序
      甲基化提取
    可视化与结果整理
      bigWig 轨道
      IGV 查看
      deepTools profile
      gene region methylation profile
      "CG / CHG / CHH 扇形图"

本教程从登录服务器开始，依次完成原始 reads 质控、trimming、Bismark 比对、甲基化提取、bigWig 轨道构建、IGV 查看、gene region methylation profile 以及 CG / CHG / CHH 扇形图绘制。

二、先登录主服务器并进入计算节点

目的

先登录服务器，再进入一个计算节点，后续分析都在计算节点中完成。

命令

先登录主服务器（具体地址和口令请联系课程 TA 获取）。

登录成功后，再进入一个计算节点。下面以 node05 为例，实际使用时可根据可用节点调整。

ssh node05

快速检查

hostname
pwd

三、开始前要知道的文件

公共参考目录中应该有：

Arabidopsis_thaliana_TAIR10_chr1.fa
Arabidopsis_thaliana_TAIR10_chr1.gff3

公共原始数据目录中应该有两个单端样本，例如：

SRR534177.train10pct.fastq.gz
SRR534239.train10pct.fastq.gz

这两个文件分别对应：

SRR534177.train10pct.fastq.gz：colWT 的 10% 子集
SRR534239.train10pct.fastq.gz：met1 的 10% 子集

这里是单端数据，因为每个样本只有一个 fastq.gz，没有 R1 / R2 成对文件。

四、创建工作目录并激活环境

目的

建立自己的分析目录，并确认 wgbs 环境已经激活。

命令

mkdir -p ~/wgbs
cd ~/wgbs

mkdir -p 00_reference 01_rawdata 02_fastqc 03_trimmed 04_bismark_index 05_alignment 06_methylation 07_igv 08_plots 10_summary scripts

conda activate wgbs

pwd
which bismark
which trim_galore
which samtools
which MethylDackel

快速检查

MethylDackel 2>&1 | head

这里先确认当前目录正确，并且 bismark / trim_galore / samtools / MethylDackel 都能在 wgbs 环境中正常找到。

五、复制参考文件和原始数据

目的

把课堂用到的参考文件和 10% 原始 reads 复制到自己的工作目录中。

命令

cp /public/home/bio2023/shiqc/dataforwgbs/ref/Arabidopsis_thaliana_TAIR10_chr1.fa 00_reference/
cp /public/home/bio2023/shiqc/dataforwgbs/ref/Arabidopsis_thaliana_TAIR10_chr1.gff3 00_reference/

cp /public/home/bio2023/shiqc/dataforwgbs/rawdata-0.1/SRR534177.train10pct.fastq.gz 01_rawdata/SRR534177_colWT.classroom.fastq.gz
cp /public/home/bio2023/shiqc/dataforwgbs/rawdata-0.1/SRR534239.train10pct.fastq.gz 01_rawdata/SRR534239_met1.classroom.fastq.gz

快速检查

ls -lh 00_reference
ls -lh 01_rawdata
zcat 01_rawdata/SRR534177_colWT.classroom.fastq.gz | head

这里应该能看到：

00_reference/ 下有 fa 和 gff3
01_rawdata/ 下有两个约 240 MB 的单端 fastq.gz
zcat ... | head 能显示标准的 FASTQ 四行结构

六、从 gff3 生成基因区间 BED6

目的

从 gff3 中提取 gene 区间，生成后续 IGV、deepTools 和基因区域 profile 会用到的 BED6 文件。

命令

awk 'BEGIN{FS=OFS="\t"} $3=="gene" {id="."; split($9,a,";"); for(i in a){if(a[i] ~ /^ID=/){sub(/^ID=/,"",a[i]); id=a[i]}} print $1,$4-1,$5,id,".",$7}' \
00_reference/Arabidopsis_thaliana_TAIR10_chr1.gff3 > 00_reference/Arabidopsis_thaliana_TAIR10_chr1.genes.bed6

快速检查

head 00_reference/Arabidopsis_thaliana_TAIR10_chr1.genes.bed6

这里使用 BED6，主要是为了同时保留 gene 区间、gene 名称和链方向，方便后面做 5' / 3' 相关的区域统计与作图。

七、原始 reads 质控

目的

先查看原始数据的基本质量情况。

命令

fastqc -o 02_fastqc 01_rawdata/SRR534177_colWT.classroom.fastq.gz
fastqc -o 02_fastqc 01_rawdata/SRR534239_met1.classroom.fastq.gz

这里的 -o 02_fastqc 表示把 FastQC 结果统一输出到 02_fastqc/ 目录。

结果文件

02_fastqc/SRR534177_colWT.classroom_fastqc.html
02_fastqc/SRR534177_colWT.classroom_fastqc.zip
02_fastqc/SRR534239_met1.classroom_fastqc.html
02_fastqc/SRR534239_met1.classroom_fastqc.zip

快速检查

ls -lh 02_fastqc

八、trimming

目的

去除 adapter、低质量碱基，并对 read 两端各额外剪去 5 bp。

命令

trim_galore --cores 2 --fastqc --clip_R1 5 --three_prime_clip_R1 5 --output_dir 03_trimmed 01_rawdata/SRR534177_colWT.classroom.fastq.gz
trim_galore --cores 2 --fastqc --clip_R1 5 --three_prime_clip_R1 5 --output_dir 03_trimmed 01_rawdata/SRR534239_met1.classroom.fastq.gz

结果文件

03_trimmed/SRR534177_colWT.classroom_trimmed.fq.gz
03_trimmed/SRR534239_met1.classroom_trimmed.fq.gz
03_trimmed/*.trimming_report.txt
trimming 后的 FastQC 文件

快速检查

ls -lh 03_trimmed
head 03_trimmed/SRR534177_colWT.classroom.fastq.gz_trimming_report.txt

从 Trim Galore 日志里，最值得注意的是：

Trimming mode: single-end
Adapter sequence: 'AGATCGGAAGAGC'
All Read 1 sequences will be trimmed by 5 bp from their 5' end
All Read 1 sequences will be trimmed by 5 bp from their 3' end
Writing final adapter and quality trimmed output to ...trimmed.fq.gz

Maximum trimming error rate: 0.1 (default) 是接头匹配时的默认容错率，不是数据本身出错。它表示在识别 adapter 时允许少量碱基不完全匹配，这在真实测序数据中是正常的。

九、Bismark 原理简介

为什么 WGBS 不能直接按普通 DNA-seq 来比对

亚硫酸氢盐处理后：

未甲基化 C 会变成 T
甲基化 C 保留为 C

所以 WGBS read 和参考基因组之间会出现大量 C/T 变化，普通 DNA-seq 的比对方式不够合适。

Bismark 在做什么

先准备 bisulfite-aware 参考索引
再调用 Bowtie2 做比对
最后直接衔接甲基化提取

本教程为什么不加 `--non_directional`

本教程按最常见的 directional 文库处理，因此 Bismark 比对时不额外加 --non_directional。

可以简单记住：

如果提取结果只有 OT / OB，通常就是 directional 文库
只有在 --non_directional 下比对率明显提高，并且后续出现 OT / CTOT / OB / CTOB 四组文件时，才考虑是 non-directional 文库

本教程当前数据只出现 OT / OB，因此按 directional 处理即可。

十、准备 Bismark 索引

目的

为完整 chr1 参考建立 Bismark 使用的 bisulfite 索引。

命令

mkdir -p 04_bismark_index/Arabidopsis_thaliana_TAIR10_chr1

ln -s ../../00_reference/Arabidopsis_thaliana_TAIR10_chr1.fa \
04_bismark_index/Arabidopsis_thaliana_TAIR10_chr1/Arabidopsis_thaliana_TAIR10_chr1.fa

bismark_genome_preparation --bowtie2 04_bismark_index/Arabidopsis_thaliana_TAIR10_chr1

这里的 --bowtie2 表示后面使用 Bowtie2 作为 Bismark 的比对后端，因此建索引时也按 Bowtie2 的格式准备参考。

快速检查

ls -lh 04_bismark_index/Arabidopsis_thaliana_TAIR10_chr1

这里使用软链接，是为了节省空间，但不要删除 00_reference 里的原始 fa。

十一、Bismark 比对

目的

把 trimming 后的 reads 比对到 chr1 参考上。

命令

可以顺序运行，也可以并行后台运行。下面给出一个常见的并行写法：

nohup bismark -p 2 --bam --bowtie2 \
  --genome 04_bismark_index/Arabidopsis_thaliana_TAIR10_chr1 \
  -o 05_alignment \
  03_trimmed/SRR534177_colWT.classroom_trimmed.fq.gz \
  > 05_alignment/SRR534177.bismark.log 2>&1 &

nohup bismark -p 2 --bam --bowtie2 \
  --genome 04_bismark_index/Arabidopsis_thaliana_TAIR10_chr1 \
  -o 05_alignment \
  03_trimmed/SRR534239_met1.classroom_trimmed.fq.gz \
  > 05_alignment/SRR534239.bismark.log 2>&1 &

这里几个参数可以简单理解为：

nohup：让任务在终端断开后继续运行
-p 2：使用 2 个线程
--bam：直接输出 BAM 格式结果
--bowtie2：使用 Bowtie2 进行比对
--genome：指定 Bismark 索引目录
-o 05_alignment：把比对结果输出到 05_alignment/
> ... 2>&1 &：把日志写入文件，并让任务在后台运行

快速检查

tail -f 05_alignment/SRR534177.bismark.log
tail -f 05_alignment/SRR534239.bismark.log

运行时如果看到 Failed to launch x86-64-v3 version, staying with default，通常只是程序回退到默认二进制版本，不影响比对结果。

十二、去重复、排序、建索引

目的

整理 BAM，为后续甲基化提取和 IGV 查看做准备。

命令

为了让 deduplicate_bismark 的输出直接落在 05_alignment/，这里先进入该目录再去重：

cd 05_alignment

deduplicate_bismark --bam SRR534177_colWT.classroom_trimmed_bismark_bt2.bam
deduplicate_bismark --bam SRR534239_met1.classroom_trimmed_bismark_bt2.bam

cd ..

samtools sort -o 05_alignment/SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.bam 05_alignment/SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.bam
samtools sort -o 05_alignment/SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted.bam 05_alignment/SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.bam

samtools index 05_alignment/SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.bam
samtools index 05_alignment/SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted.bam

这里的 deduplicate_bismark --bam 表示输入文件是单端 BAM，程序会按比对位置去除可能的 PCR duplicate。

这里的 samtools sort 默认按染色体坐标排序，而不是按 read 名字排序。n 表示按数值排序，避免把 1、10、100、2 这种坐标按字符顺序排错。

结果文件

05_alignment/*.deduplicated.sorted.bam
05_alignment/*.deduplicated.sorted.bam.bai

快速检查

samtools flagstat 05_alignment/SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.bam
samtools flagstat 05_alignment/SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted.bam

samtools coverage 05_alignment/SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.bam
samtools coverage 05_alignment/SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted.bam

这里可以简单这样看：

in total：当前 BAM 中的总 reads 数
mapped：当前保留下来的 reads 中成功比对到参考的数量
单端数据里 paired in sequencing / read1 / read2 通常都是 0
coverage：chr1 上被覆盖到的碱基比例
meandepth：chr1 全长平均下来的覆盖深度

这里看到的 mapped 100%，指的是当前这个 deduplicated sorted BAM 相对于 chr1 参考的结果，不是完整基因组层面的总体比对率。

可以直接拿 SRR534177_colWT 这个样本来理解：

63820 in total：这个 BAM 里一共有 63820 条保留下来的 reads
63820 mapped：这些 reads 都成功比对到了当前的 chr1 参考
coverage = 11.7251：完整 chr1 上大约有 11.7% 的碱基至少被一条 read 覆盖到
meandepth = 0.144453：如果把 chr1 全长平均来看，平均覆盖深度约为 0.14x

这说明这套 10% 教学数据已经能在 chr1 上提供一定覆盖，但仍然不是高深度全覆盖数据，因此后面的轨道和 profile 图会有一定稀疏性。

十三、甲基化提取

目的

提取甲基化信息，并为后面的轨道和统计结果做准备。

这里分成两部分：

先用 bismark_methylation_extractor 保留 Bismark context 文件，供后面的 gene_region_profile.py 和扇形图统计使用
再用 MethylDackel 直接生成 CG / CHG / CHH 三种 bedGraph，供后面的 IGV 和 deepTools 使用

命令

先运行 Bismark 的 context 提取。这里把 genome_folder 写成绝对路径：

GENOME_DIR="$(pwd)/04_bismark_index/Arabidopsis_thaliana_TAIR10_chr1"

bismark_methylation_extractor \
  --single-end \
  --bedGraph \
  --counts \
  --gzip \
  --genome_folder "$GENOME_DIR" \
  -o 06_methylation \
  05_alignment/SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.bam

bismark_methylation_extractor \
  --single-end \
  --bedGraph \
  --counts \
  --gzip \
  --genome_folder "$GENOME_DIR" \
  -o 06_methylation \
  05_alignment/SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted.bam

这里先不额外设置 coverage cutoff，先把 Bismark 的 context 结果完整提出来，后面如果需要更严格过滤，再单独加阈值。

如果后续想做更严格的过滤，可以考虑这些可选参数：

bismark_methylation_extractor --cutoff N：只在 bedGraph / coverage 输出中保留覆盖度至少为 N 的位点，例如 --cutoff 3 或 --cutoff 5
MethylDackel extract -d N：设置最小覆盖度，例如 -d 4
MethylDackel extract -q N：设置最小 mapping quality，例如 -q 10
MethylDackel extract -p N：设置最小 base quality，例如 -p 5

本教程当前数据量不算高，因此主线先不加这些过滤参数，避免过早丢掉信号。

接着使用 MethylDackel 直接生成 CG / CHG / CHH 三种 bedGraph：

MethylDackel extract --CHG --CHH --fraction \
  -o 06_methylation/SRR534177_colWT \
  00_reference/Arabidopsis_thaliana_TAIR10_chr1.fa \
  05_alignment/SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.bam

MethylDackel extract --CHG --CHH --fraction \
  -o 06_methylation/SRR534239_met1 \
  00_reference/Arabidopsis_thaliana_TAIR10_chr1.fa \
  05_alignment/SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted.bam

这里把 CG / CHG / CHH 都放在同一条 MethylDackel extract 命令里，目的是让三种 context 的提取标准保持一致，后面更便于直接比较。

这里的参数可以简单理解为：

extract：从排序后的 BAM 中提取甲基化信息
--CHG --CHH：除了 CpG 之外，也一起输出 CHG 和 CHH
--fraction：输出按比例表示的甲基化信号
-o：设置输出文件前缀

后面的 IGV 和 deepTools 轨道也都统一使用 MethylDackel 生成的 CG / CHG / CHH bedGraph，不再把 Bismark 的 CpG bedGraph 和 MethylDackel 的 CHG / CHH 混在一起。

结果文件

06_methylation/CpG_OB_*.txt.gz
06_methylation/CpG_OT_*.txt.gz
06_methylation/CHG_OB_*.txt.gz
06_methylation/CHG_OT_*.txt.gz
06_methylation/CHH_OB_*.txt.gz
06_methylation/CHH_OT_*.txt.gz
06_methylation/*.bedGraph.gz
06_methylation/*.bismark.cov.gz
06_methylation/*_CpG.meth.bedGraph
06_methylation/*_CHG.meth.bedGraph
06_methylation/*_CHH.meth.bedGraph
06_methylation/*_splitting_report.txt
06_methylation/*.M-bias.txt

这里的 OT 和 OB 可以简单理解为链特异性的甲基化调用：

OT：Original Top，对应参考基因组原始正链上的信息
OB：Original Bottom，对应参考基因组原始负链上的信息

如果是少数 non-directional 文库，还会出现 CTOT 和 CTOB；它们分别是 Original Top 和 Original Bottom 的互补链信息。

之所以分开，是因为双硫酸盐处理后正负链的转换方式不同，Bismark 需要先按链分别记录，后面再合并成位点层面的甲基化比例。因此：

看链特异性结果时，可以分别看 OT 和 OB
看后面的 bedGraph / cov / bigWig 时，通常已经是合并后的结果

*_splitting_report.txt 是甲基化提取的汇总报告，适合快速看 CpG / CHG / CHH 的总体统计。*.M-bias.txt 则是在看 read 不同位置上的甲基化比例是否存在系统性偏差，也就是常说的 M-bias。这和前面保留 --clip_R1 5 --three_prime_clip_R1 5 的 trimming 设置是有关联的。

快速检查

ls -lh 06_methylation
head 06_methylation/SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.M-bias.txt

十四、生成 IGV 轨道

目的

把结果整理成 BAM + BAI + bigWig，用于 IGV。

这一节直接使用上一节 MethylDackel 生成的 CpG / CHG / CHH bedGraph 转成 bigWig。

命令

先准备染色体长度文件：

mkdir -p 07_igv

samtools faidx 00_reference/Arabidopsis_thaliana_TAIR10_chr1.fa
cut -f1,2 00_reference/Arabidopsis_thaliana_TAIR10_chr1.fa.fai > 07_igv/Arabidopsis_thaliana_TAIR10_chr1.chrom.sizes

再把 MethylDackel 生成的三个 bedGraph 整理成干净的 bedGraph，然后转成 bigWig：

grep -v '^track' 06_methylation/SRR534177_colWT_CpG.meth.bedGraph > 07_igv/SRR534177_colWT.CpG.clean.bedGraph
grep -v '^track' 06_methylation/SRR534177_colWT_CHG.meth.bedGraph > 07_igv/SRR534177_colWT.CHG.clean.bedGraph
grep -v '^track' 06_methylation/SRR534177_colWT_CHH.meth.bedGraph > 07_igv/SRR534177_colWT.CHH.clean.bedGraph

grep -v '^track' 06_methylation/SRR534239_met1_CpG.meth.bedGraph > 07_igv/SRR534239_met1.CpG.clean.bedGraph
grep -v '^track' 06_methylation/SRR534239_met1_CHG.meth.bedGraph > 07_igv/SRR534239_met1.CHG.clean.bedGraph
grep -v '^track' 06_methylation/SRR534239_met1_CHH.meth.bedGraph > 07_igv/SRR534239_met1.CHH.clean.bedGraph

bedGraphToBigWig 07_igv/SRR534177_colWT.CpG.clean.bedGraph 07_igv/Arabidopsis_thaliana_TAIR10_chr1.chrom.sizes 07_igv/SRR534177_colWT.CpG.bw
bedGraphToBigWig 07_igv/SRR534177_colWT.CHG.clean.bedGraph 07_igv/Arabidopsis_thaliana_TAIR10_chr1.chrom.sizes 07_igv/SRR534177_colWT.CHG.bw
bedGraphToBigWig 07_igv/SRR534177_colWT.CHH.clean.bedGraph 07_igv/Arabidopsis_thaliana_TAIR10_chr1.chrom.sizes 07_igv/SRR534177_colWT.CHH.bw

bedGraphToBigWig 07_igv/SRR534239_met1.CpG.clean.bedGraph 07_igv/Arabidopsis_thaliana_TAIR10_chr1.chrom.sizes 07_igv/SRR534239_met1.CpG.bw
bedGraphToBigWig 07_igv/SRR534239_met1.CHG.clean.bedGraph 07_igv/Arabidopsis_thaliana_TAIR10_chr1.chrom.sizes 07_igv/SRR534239_met1.CHG.bw
bedGraphToBigWig 07_igv/SRR534239_met1.CHH.clean.bedGraph 07_igv/Arabidopsis_thaliana_TAIR10_chr1.chrom.sizes 07_igv/SRR534239_met1.CHH.bw

这里统一使用 MethylDackel 生成的 CG / CHG / CHH bedGraph，这样三种 context 的轨道来自同一套提取标准，后面在 IGV 和 deepTools 里更容易直接比较。其中 grep -v '^track' 是去掉 bedGraph 文件开头的轨道说明行，bedGraphToBigWig 则是把文本格式的 bedGraph 转成更适合 IGV 和后续作图使用的 bigWig。

结果文件

07_igv/Arabidopsis_thaliana_TAIR10_chr1.chrom.sizes
07_igv/SRR534177_colWT.CpG.bw
07_igv/SRR534177_colWT.CHG.bw
07_igv/SRR534177_colWT.CHH.bw
07_igv/SRR534239_met1.CpG.bw
07_igv/SRR534239_met1.CHG.bw
07_igv/SRR534239_met1.CHH.bw

快速检查

ls -lh 07_igv

十五、IGV 查看

目的

把比对结果和甲基化轨道加载到 IGV 中。

载入文件

05_alignment/SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.bam
05_alignment/SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.bam.bai
07_igv/SRR534177_colWT.CpG.bw
07_igv/SRR534177_colWT.CHG.bw
07_igv/SRR534177_colWT.CHH.bw

另一个样本同理。

十七、deepTools 版本

目的

这一部分是 deepTools 版本，作为补充放在最后。主线结果以上面的 gene_region_profile.py 为主。

命令

这里的参数含义可以简单理解为：

每个样本只画一张图，图中三条线分别表示 CG、CHG 和 CHH
-S：输入的 bigWig 轨道，这里依次放入 CG / CHG / CHH
-R：输入的基因区域文件，这里使用全部 chr1 gene
-b 1000：每个基因上游取 1000 bp
-a 1000：每个基因下游取 1000 bp
--regionBodyLength 1000：把 gene body 统一缩放成 1000 bp
--binSize 100：每 100 bp 统计一个 bin，因此上游、gene body、下游各得到 10 个 bin
--skipZeros：跳过全为 0 的区域，减少空信号对图形的影响
--perGroup：把同一样本的 CG / CHG / CHH 三条线画在同一张图里

这里先运行 computeMatrix，是因为 plotProfile 不能直接拿 bigWig + BED 画图。computeMatrix scale-regions 会先把所有基因统一按 5' -> gene body -> 3' 的结构整理成矩阵：

每一行代表一个基因
每一列代表一个 bin
数值是对应 bin 中的平均甲基化信号

这里的 scale-regions 表示把不同长度的 gene body 统一缩放到同样长度，再与上下游一起比较。按当前参数：

上游 1000 bp 分成 10 个 bin
gene body 缩放成 1000 bp 后分成 10 个 bin
下游 1000 bp 分成 10 个 bin

因此每个 context 都有 30 个 bin，CG / CHG / CHH 三种 context 合起来就是一份可供后面 plotProfile 直接作图的矩阵文件。

这里图上的 y 轴是 deepTools 对 bigWig 信号按区域分箱后再做平均得到的结果，不再是单个位点的原始甲基化百分比。由于这套 chr1 教学数据覆盖较低，平均值通常会比较小。

computeMatrix scale-regions \
  -S 07_igv/SRR534177_colWT.CpG.bw 07_igv/SRR534177_colWT.CHG.bw 07_igv/SRR534177_colWT.CHH.bw \
  -R 00_reference/Arabidopsis_thaliana_TAIR10_chr1.genes.bed6 \
  -b 1000 -a 1000 \
  --regionBodyLength 1000 \
  --binSize 100 \
  --skipZeros \
  -o 08_plots/SRR534177_colWT.deeptools.matrix.gz

plotProfile \
  -m 08_plots/SRR534177_colWT.deeptools.matrix.gz \
  --perGroup \
  --samplesLabel CG CHG CHH \
  --yAxisLabel "Mean methylation signal" \
  --plotTitle "Col-0 deepTools profile" \
  -out 08_plots/SRR534177_colWT.deeptools_profile.png

computeMatrix scale-regions \
  -S 07_igv/SRR534239_met1.CpG.bw 07_igv/SRR534239_met1.CHG.bw 07_igv/SRR534239_met1.CHH.bw \
  -R 00_reference/Arabidopsis_thaliana_TAIR10_chr1.genes.bed6 \
  -b 1000 -a 1000 \
  --regionBodyLength 1000 \
  --binSize 100 \
  --skipZeros \
  -o 08_plots/SRR534239_met1.deeptools.matrix.gz

plotProfile \
  -m 08_plots/SRR534239_met1.deeptools.matrix.gz \
  --perGroup \
  --samplesLabel CG CHG CHH \
  --yAxisLabel "Mean methylation signal" \
  --plotTitle "met1 deepTools profile" \
  -out 08_plots/SRR534239_met1.deeptools_profile.png

如果想直接看看 matrix.gz 里面装的是什么，可以这样做：

zcat 08_plots/SRR534177_colWT.deeptools.matrix.gz | head -n 3

输出通常可以这样理解：

第 1 行是一个 JSON 头信息，记录了这份矩阵的参数，例如：
- sample_labels：这里对应 CG / CHG / CHH
- sample_boundaries:[0,30,60,90]：表示 90 列信号被分成了 3 段，每段 30 列
- group_boundaries:[0,2014]：表示这里总共统计了 2014 个基因
从第 2 行开始，每一行对应一个基因
前 6 列是原始 BED6 信息：chr / start / end / gene_id / score / strand
后面的 90 列就是这个基因在 CG / CHG / CHH 三种 context 下的 bin 信号：
- 第 1-30 列：CG
- 第 31-60 列：CHG
- 第 61-90 列：CHH

例如第 2 行开头的：

1  23145  31227  AT1G01040  .  +

说明这一行对应的是 chr1 上的基因 AT1G01040，链方向是 +。

后面出现的很多 nan，表示该 bin 在当前数据里没有足够信号可用于统计；像 0.000000、1.000000 这样的值则表示该 bin 的平均甲基化信号。由于这套 chr1 教学数据覆盖有限，所以矩阵里出现较多 nan 是正常的。

结果文件

08_plots/SRR534177_colWT.deeptools.matrix.gz
08_plots/SRR534177_colWT.deeptools_profile.png
08_plots/SRR534239_met1.deeptools.matrix.gz
08_plots/SRR534239_met1.deeptools_profile.png

十八、统计 CG / CHG / CHH 并画扇形图

目的

查看每个样本内部三种 context 的组成比例。这里直接读取 Bismark 的 splitting_report.txt 汇总结果来作图。

命令

python - <<'PY'
from pathlib import Path
import matplotlib.pyplot as plt

report_files = {
    "SRR534177_colWT": Path("06_methylation/SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted_splitting_report.txt"),
    "SRR534239_met1": Path("06_methylation/SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted_splitting_report.txt"),
}

for sample, report_path in report_files.items():
    counts = {"CG": 0, "CHG": 0, "CHH": 0}
    with report_path.open("r", encoding="utf-8") as fh:
        for line in fh:
            # splitting_report 中已经汇总了三种 context 的甲基化 C 数量
            if line.startswith("Total methylated C's in CpG context:"):
                counts["CG"] = int(line.rstrip().split("\t")[-1])
            elif line.startswith("Total methylated C's in CHG context:"):
                counts["CHG"] = int(line.rstrip().split("\t")[-1])
            elif line.startswith("Total methylated C's in CHH context:"):
                counts["CHH"] = int(line.rstrip().split("\t")[-1])

    # 按固定顺序整理成 matplotlib 画图需要的输入
    labels = ["CG", "CHG", "CHH"]
    values = [counts[x] for x in labels]

    plt.figure(figsize=(4, 4))
    plt.pie(values, labels=labels, autopct="%1.1f%%", startangle=90)
    plt.title(sample + " methylated cytosine context distribution")
    plt.tight_layout()
    plt.savefig(f"08_plots/{sample}.context_pie.png", dpi=200)
    plt.close()
PY

结果文件

08_plots/SRR534177_colWT.context_pie.png
08_plots/SRR534239_met1.context_pie.png

快速检查

ls -lh 08_plots

这里直接读取 splitting_report.txt 中已经汇总好的 Total methylated C's in CpG/CHG/CHH context 三行结果，因此比逐个扫描 CpG_OB/OT、CHG_OB/OT、CHH_OB/OT 文件更简单。这种写法适合画样本整体的扇形图；如果要做位点级或区域级统计，仍然需要前面的 context 原始文件。

十九、最后应该得到哪些结果

这门课的主线完成后，你至少应该得到：

00_reference/Arabidopsis_thaliana_TAIR10_chr1.genes.bed6
03_trimmed/*_trimmed.fq.gz
05_alignment/*.deduplicated.sorted.bam
05_alignment/*.deduplicated.sorted.bam.bai
06_methylation/CpG_OB_*.txt.gz
06_methylation/CpG_OT_*.txt.gz
06_methylation/CHG_OB_*.txt.gz
06_methylation/CHG_OT_*.txt.gz
06_methylation/CHH_OB_*.txt.gz
06_methylation/CHH_OT_*.txt.gz
06_methylation/*_CpG.meth.bedGraph
06_methylation/*_CHG.meth.bedGraph
06_methylation/*_CHH.meth.bedGraph
07_igv/*.bw
08_plots/*official_like_profile.norm.png
08_plots/*.context_pie.png

如果继续补充做 deepTools，还会额外得到：

08_plots/*deeptools_profile.png

附录、gene_region_profile.py 脚本说明

这里真正调用的是 scripts/gene_region_profile.py。它是一个普通的 Python 脚本，不依赖 deepTools，而是直接读取 Bismark 的 context 提取结果，自己完成分箱统计和作图。

统计口径固定为：

每个样本只画一张图，图中三条线分别表示 CG、CHG 和 CHH
全部 chr1 gene
上游 1000 bp
gene body 1000 bp
下游 1000 bp
10 / 10 / 10 个 bin
图上只标 5' 和 3'

这个脚本内部主要做了这几件事：

先读取 genes.bed6，拿到每个 gene 的起止位置和链方向
再读取 CpG_OB/OT、CHG_OB/OT、CHH_OB/OT，把同一位点的甲基化和未甲基化计数合并起来
以每个 gene 为单位，按 5' -> gene body -> 3' 的统一方向统计
上游 1000 bp 分成 10 个 bin，gene body 缩放成 1000 bp 后再分成 10 个 bin，下游 1000 bp 也分成 10 个 bin
对所有 gene 的对应 bin 累加 meth / unmeth 计数，再计算每个 bin 的甲基化比例
最后输出一份 tsv 数值表，并画出每个样本一张、包含 CG / CHG / CHH 三条线的 profile 图

其中：

--cpg：输入 CpG 的 context 文件
--noncpg：输入 CHG / CHH 的 context 文件
--normalize-max 表示把每条曲线按自身最大值归一化，便于在一张图里比较形状
--smooth-window 3 表示对相邻 bin 做一个轻微平滑，让曲线更容易观察总体趋势

命令

python scripts/gene_region_profile.py \
  --genes 00_reference/Arabidopsis_thaliana_TAIR10_chr1.genes.bed6 \
  --cpg 06_methylation/CpG_OB_SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz 06_methylation/CpG_OT_SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz \
  --noncpg 06_methylation/CHG_OB_SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz 06_methylation/CHG_OT_SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz 06_methylation/CHH_OB_SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz 06_methylation/CHH_OT_SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz \
  --sample SRR534177_colWT \
  --title "Col-0 (normalized to max)" \
  --out-tsv 08_plots/SRR534177_colWT.official_like_profile.tsv \
  --out-pdf 08_plots/SRR534177_colWT.official_like_profile.norm.png \
  --normalize-max \
  --smooth-window 3

python scripts/gene_region_profile.py \
  --genes 00_reference/Arabidopsis_thaliana_TAIR10_chr1.genes.bed6 \
  --cpg 06_methylation/CpG_OB_SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz 06_methylation/CpG_OT_SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz \
  --noncpg 06_methylation/CHG_OB_SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz 06_methylation/CHG_OT_SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz 06_methylation/CHH_OB_SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz 06_methylation/CHH_OT_SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz \
  --sample SRR534239_met1 \
  --title "met1 (normalized to max)" \
  --out-tsv 08_plots/SRR534239_met1.official_like_profile.tsv \
  --out-pdf 08_plots/SRR534239_met1.official_like_profile.norm.png \
  --normalize-max \
  --smooth-window 3

附录、常见命令行参数速查

下面这些参数不是 WGBS 特有参数，但在本教程里会反复出现：

ssh -p <port>：-p 表示指定远程登录端口，课堂服务器端口以 TA 通知为准
mkdir -p：-p 表示目录不存在时自动创建，已存在也不会报错
ln -s：-s 表示创建软链接，而不是复制原文件
tail -f：-f 表示持续跟踪文件末尾内容，常用来看日志是否还在更新
samtools faidx：给 FASTA 建立索引，并生成对应的 .fai 文件
cut -f1,2：取文本的第 1 和第 2 列，这里用来从 .fai 中提取染色体名和长度

课前信息#

联系方式（Teacher / TA）#

参考链接#

远程登录工具#

一、课堂数据来源与参考#

二、先登录主服务器并进入计算节点#

目的#

命令#

快速检查#

三、开始前要知道的文件#

四、创建工作目录并激活环境#

目的#

命令#

快速检查#

五、复制参考文件和原始数据#

目的#

命令#

快速检查#

六、从 gff3 生成基因区间 BED6#

目的#

命令#

快速检查#

七、原始 reads 质控#

目的#

命令#

结果文件#

快速检查#

八、trimming#

目的#

命令#

结果文件#

快速检查#

九、Bismark 原理简介#

为什么 WGBS 不能直接按普通 DNA-seq 来比对#

Bismark 在做什么#

本教程为什么不加 --non_directional#

十、准备 Bismark 索引#

目的#

命令#

快速检查#

十一、Bismark 比对#

目的#

命令#

快速检查#

十二、去重复、排序、建索引#

目的#

命令#

结果文件#

快速检查#

十三、甲基化提取#

目的#

命令#

结果文件#

快速检查#

十四、生成 IGV 轨道#

目的#

命令#

结果文件#

快速检查#

十五、IGV 查看#

目的#

载入文件#

十七、deepTools 版本#

目的#

命令#

结果文件#

十八、统计 CG / CHG / CHH 并画扇形图#

目的#

命令#

结果文件#

快速检查#

十九、最后应该得到哪些结果#

附录、gene_region_profile.py 脚本说明#

命令#

附录、常见命令行参数速查#

一起交流

课前信息

联系方式（Teacher / TA）

参考链接

远程登录工具

一、课堂数据来源与参考

二、先登录主服务器并进入计算节点

目的

命令

快速检查

三、开始前要知道的文件

四、创建工作目录并激活环境

目的

命令

快速检查

五、复制参考文件和原始数据

目的

命令

快速检查

六、从 gff3 生成基因区间 BED6

目的

命令

快速检查

七、原始 reads 质控

目的

命令

结果文件

快速检查

八、trimming

目的

命令

结果文件

快速检查

九、Bismark 原理简介

为什么 WGBS 不能直接按普通 DNA-seq 来比对

Bismark 在做什么

本教程为什么不加 `--non_directional`

十、准备 Bismark 索引

目的

命令

快速检查

十一、Bismark 比对

目的

命令

快速检查

十二、去重复、排序、建索引

目的

命令

结果文件

快速检查

十三、甲基化提取

目的

命令

结果文件

快速检查

十四、生成 IGV 轨道

目的

命令

结果文件

快速检查

十五、IGV 查看

目的

载入文件

十七、deepTools 版本

目的

命令

结果文件

十八、统计 CG / CHG / CHH 并画扇形图

目的

命令

结果文件

快速检查

十九、最后应该得到哪些结果

附录、gene_region_profile.py 脚本说明

命令

附录、常见命令行参数速查