DNA 甲基化数据分析

这是一份 WGBS 课堂教程笔记，示例数据为 Arabidopsis thaliana TAIR10 chr1 教学数据。流程覆盖 single-end WGBS 从原始 reads 到甲基化提取、IGV 轨道和基因区域 profile 的主要步骤。

教程信息

数据类型：single-end WGBS
参考基因组：Arabidopsis thaliana TAIR10 chr1
示例样本：SRR534177_colWT 和 SRR534239_met1
软件环境示例：wgbs

外部工具：

IGV 下载：https://igv.org/download/
MobaXterm：https://mobaxterm.mobatek.net/
FileZilla：https://filezilla-project.org/

数据来源与参考

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 run 在 NCBI 中均标记为 Layout: SINGLE，因此本教程按单端 single-end 数据处理。

流程总览

登录服务器并进入计算节点。
创建个人工作目录并激活 WGBS 环境。
准备参考文件、原始 reads 和课堂脚本。
从 GFF3 生成 genes.bed6。
使用 FastQC 检查原始 reads。
使用 Trim Galore 去接头和低质量末端。
使用 Bismark 建索引、比对、去重复和甲基化提取。
生成 BAM / BAI / bigWig 轨道并在 IGV 中查看。
绘制 gene region methylation profile。
统计 CG / CHG / CHH context 并绘制扇形图。

登录服务器并进入计算节点

后续分析建议在计算节点完成。

ssh -p <PORT> <USER>@<HOST>
ssh <NODE_NAME>

快速检查：

hostname
pwd

开始前要知道的文件

公共参考目录中应包含：

Arabidopsis_thaliana_TAIR10_chr1.fa
Arabidopsis_thaliana_TAIR10_chr1.gff3

公共原始数据目录中应包含两个单端样本：

SRR534177_colWT.classroom.fastq.gz
SRR534239_met1.classroom.fastq.gz

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

创建工作目录

在自己的目录中完成整个分析，不直接在公共目录里运行。

mkdir -p ~/<YOUR_NAME>/wgbs
cd ~/<YOUR_NAME>/wgbs
conda activate wgbs

快速检查：

pwd
which bismark
which trim_galore
which samtools

创建分析目录

把每一步的结果分目录保存，方便后续检查。

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

快速检查：

ls

复制参考文件和原始 reads

把公共数据复制到自己的工作目录。下面的 <CLASS_DATA_DIR> 请替换为课堂共享数据目录。

cp <CLASS_DATA_DIR>/ref/Arabidopsis_thaliana_TAIR10_chr1.fa 00_reference/
cp <CLASS_DATA_DIR>/ref/Arabidopsis_thaliana_TAIR10_chr1.gff3 00_reference/

cp <CLASS_DATA_DIR>/rawdata-chr1only/SRR534177_colWT.classroom.fastq.gz 01_rawdata/
cp <CLASS_DATA_DIR>/rawdata-chr1only/SRR534239_met1.classroom.fastq.gz 01_rawdata/

结果文件：

00_reference/Arabidopsis_thaliana_TAIR10_chr1.fa
00_reference/Arabidopsis_thaliana_TAIR10_chr1.gff3
01_rawdata/SRR534177_colWT.classroom.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

复制课堂脚本

后面生成 gene region profile 图时会用到两个脚本。下面的 <CLASS_SCRIPT_DIR> 请替换为课堂共享脚本目录。

cp <CLASS_SCRIPT_DIR>/gene_region_profile.py scripts/
cp <CLASS_SCRIPT_DIR>/context_calls_to_bedgraph.py scripts/

快速检查：

ls scripts

从 GFF3 提取 chr1 gene BED

后续区域统计和作图更适合使用 BED，所以先从 GFF3 提取 gene 区间。

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

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

快速检查：

head 00_reference/Arabidopsis_thaliana_TAIR10_chr1.genes.bed6
wc -l 00_reference/Arabidopsis_thaliana_TAIR10_chr1.genes.bed6

如果提取正确，这里应该看到大约 7509 个 gene。

FastQC 质控

先查看原始 reads 的质量情况。

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

结果文件：

02_fastqc/*html
02_fastqc/*zip

快速检查：

ls -lh 02_fastqc

Trim Galore

去除接头和低质量末端，得到更适合比对的 reads。

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

这里的 --clip_R1 5 --three_prime_clip_R1 5 用于处理 WGBS 常见的两端偏差。--clip_R1 5 表示从 read 的 5’ 端额外剪掉 5 bp，--three_prime_clip_R1 5 表示从 3’ 端也额外剪掉 5 bp。

结果文件：

03_trimmed/*_trimmed.fq.gz
03_trimmed/*html

快速检查：

ls -lh 03_trimmed

Bismark 原理简介

亚硫酸氢盐处理后：

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

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

Bismark 的主要工作包括：

准备 bisulfite-aware 参考索引。
调用 Bowtie2 做比对。
衔接甲基化位点提取。

准备 Bismark 参考索引

为 chr1 参考基因组建立 Bismark 索引。

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

这里使用软链接而不是再次复制 fa，可以节省空间；但要注意，00_reference 里的原始 fa 不能删除。

快速检查：

ls 04_bismark_index/Arabidopsis_thaliana_TAIR10_chr1

Bismark 比对

把 trimmed reads 比对到 chr1 参考基因组。

bismark --bowtie2 --bam --genome 04_bismark_index/Arabidopsis_thaliana_TAIR10_chr1 -o 05_alignment 03_trimmed/SRR534177_colWT.classroom_trimmed.fq.gz
bismark --bowtie2 --bam --genome 04_bismark_index/Arabidopsis_thaliana_TAIR10_chr1 -o 05_alignment 03_trimmed/SRR534239_met1.classroom_trimmed.fq.gz

结果文件：

05_alignment/*_bismark_bt2.bam
05_alignment/*_report.txt

快速检查：

ls -lh 05_alignment

去重复、排序、建索引

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

deduplicate_bismark --bam 05_alignment/SRR534177_colWT.classroom_trimmed_bismark_bt2.bam
deduplicate_bismark --bam 05_alignment/SRR534239_met1.classroom_trimmed_bismark_bt2.bam

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

结果文件：

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

flagstat 用来快速看总体比对情况，coverage 可以帮助查看 chr1 上的大致覆盖范围和平均深度。

Bismark 甲基化提取

提取每个位点的甲基化信息。

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 3 或 --cutoff 5。

结果文件：

06_methylation/CpG_context_*.txt.gz
06_methylation/Non_CpG_context_*.txt.gz
06_methylation/*.bedGraph.gz
06_methylation/*.bismark.cov.gz

快速检查：

ls -lh 06_methylation

生成 IGV 轨道

把结果整理成 BAM + BAI + bigWig，用于 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

解压 CpG bedGraph 并转成 bigWig：

zcat 06_methylation/SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.bedGraph.gz > 07_igv/SRR534177_colWT.CpG.bedGraph
zcat 06_methylation/SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted.bedGraph.gz > 07_igv/SRR534239_met1.CpG.bedGraph

grep -v '^track' 07_igv/SRR534177_colWT.CpG.bedGraph > 07_igv/SRR534177_colWT.CpG.clean.bedGraph
grep -v '^track' 07_igv/SRR534239_met1.CpG.bedGraph > 07_igv/SRR534239_met1.CpG.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/SRR534239_met1.CpG.clean.bedGraph 07_igv/Arabidopsis_thaliana_TAIR10_chr1.chrom.sizes 07_igv/SRR534239_met1.CpG.bw

从 Non_CpG_context 中提取 CHG 和 CHH。第 5 列区分不同类型的非 CpG 位点：

x：CHG 位点，未甲基化
X：CHG 位点，甲基化
h：CHH 位点，未甲基化
H：CHH 位点，甲基化

zcat 06_methylation/Non_CpG_context_SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz | \
awk 'BEGIN{FS=OFS="\t"} /^Bismark methylation extractor version/{next} {key=$3 OFS $4; if($5=="X") meth[key]++; else if($5=="x") unmeth[key]++} END{for(k in meth) seen[k]=1; for(k in unmeth) seen[k]=1; for(k in seen){split(k,a,OFS); cov=meth[k]+unmeth[k]; if(cov>0) print a[1], a[2]-1, a[2], 100*meth[k]/cov}}' | sort -k1,1 -k2,2n > 07_igv/SRR534177_colWT.CHG.bedGraph

zcat 06_methylation/Non_CpG_context_SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz | \
awk 'BEGIN{FS=OFS="\t"} /^Bismark methylation extractor version/{next} {key=$3 OFS $4; if($5=="H") meth[key]++; else if($5=="h") unmeth[key]++} END{for(k in meth) seen[k]=1; for(k in unmeth) seen[k]=1; for(k in seen){split(k,a,OFS); cov=meth[k]+unmeth[k]; if(cov>0) print a[1], a[2]-1, a[2], 100*meth[k]/cov}}' | sort -k1,1 -k2,2n > 07_igv/SRR534177_colWT.CHH.bedGraph

zcat 06_methylation/Non_CpG_context_SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz | \
awk 'BEGIN{FS=OFS="\t"} /^Bismark methylation extractor version/{next} {key=$3 OFS $4; if($5=="X") meth[key]++; else if($5=="x") unmeth[key]++} END{for(k in meth) seen[k]=1; for(k in unmeth) seen[k]=1; for(k in seen){split(k,a,OFS); cov=meth[k]+unmeth[k]; if(cov>0) print a[1], a[2]-1, a[2], 100*meth[k]/cov}}' | sort -k1,1 -k2,2n > 07_igv/SRR534239_met1.CHG.bedGraph

zcat 06_methylation/Non_CpG_context_SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz | \
awk 'BEGIN{FS=OFS="\t"} /^Bismark methylation extractor version/{next} {key=$3 OFS $4; if($5=="H") meth[key]++; else if($5=="h") unmeth[key]++} END{for(k in meth) seen[k]=1; for(k in unmeth) seen[k]=1; for(k in seen){split(k,a,OFS); cov=meth[k]+unmeth[k]; if(cov>0) print a[1], a[2]-1, a[2], 100*meth[k]/cov}}' | sort -k1,1 -k2,2n > 07_igv/SRR534239_met1.CHH.bedGraph

转成 bigWig：

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

结果文件：

05_alignment/*.bam
05_alignment/*.bam.bai
07_igv/*.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 profile 可选流程

先用 bigWig + deepTools 生成一版基因区域 profile 图。

参数含义：

-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。
--skipZeros：跳过全为 0 的区域，减少空信号对图形的影响。

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 \
  --samplesLabel CG CHG CHH \
  --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 \
  --samplesLabel CG CHG CHH \
  --plotTitle "met1 deepTools profile" \
  -out 08_plots/SRR534239_met1.deeptools_profile.png

结果文件：

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

全部 chr1 gene 的基因区域 profile 图

这一部分使用 gene_region_profile.py 统计并作图。

统计口径：

全部 chr1 gene
上游 1000 bp
gene body 1000 bp
下游 1000 bp
10 / 10 / 10 个 bin
每个样本一张图
图中显示 CG / CHG / CHH
图上只标 5’ 和 3'

python scripts/gene_region_profile.py \
  --genes 00_reference/Arabidopsis_thaliana_TAIR10_chr1.genes.bed6 \
  --cpg 06_methylation/CpG_context_SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz \
  --noncpg 06_methylation/Non_CpG_context_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_context_SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz \
  --noncpg 06_methylation/Non_CpG_context_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

快速检查：

ls -lh 08_plots
head 08_plots/SRR534177_colWT.official_like_profile.tsv
head 08_plots/SRR534239_met1.official_like_profile.tsv

统计 CG / CHG / CHH 并画扇形图

查看每个样本内部三种 context 的组成比例。

import gzip
from collections import Counter
import matplotlib.pyplot as plt

samples = {
    "SRR534177_colWT": "06_methylation/Non_CpG_context_SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz",
    "SRR534239_met1": "06_methylation/Non_CpG_context_SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz",
}

cpg_files = {
    "SRR534177_colWT": "06_methylation/CpG_context_SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz",
    "SRR534239_met1": "06_methylation/CpG_context_SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz",
}

for sample in samples:
    counts = Counter()
    with gzip.open(cpg_files[sample], "rt") as fh:
        for line in fh:
            if line.startswith("Bismark methylation extractor version"):
                continue
            parts = line.rstrip("\n").split("\t")
            if len(parts) == 5 and parts[4] == "Z":
                counts["CG"] += 1

    with gzip.open(samples[sample], "rt") as fh:
        for line in fh:
            if line.startswith("Bismark methylation extractor version"):
                continue
            parts = line.rstrip("\n").split("\t")
            if len(parts) != 5:
                continue
            if parts[4] == "X":
                counts["CHG"] += 1
            elif parts[4] == "H":
                counts["CHH"] += 1

    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()

也可以直接在 shell 中运行：

python - <<'PY'
import gzip
from collections import Counter
import matplotlib.pyplot as plt

samples = {
    "SRR534177_colWT": "06_methylation/Non_CpG_context_SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz",
    "SRR534239_met1": "06_methylation/Non_CpG_context_SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz",
}
cpg_files = {
    "SRR534177_colWT": "06_methylation/CpG_context_SRR534177_colWT.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz",
    "SRR534239_met1": "06_methylation/CpG_context_SRR534239_met1.classroom_trimmed_bismark_bt2.deduplicated.sorted.txt.gz",
}
for sample in samples:
    counts = Counter()
    with gzip.open(cpg_files[sample], "rt") as fh:
        for line in fh:
            if line.startswith("Bismark methylation extractor version"):
                continue
            parts = line.rstrip("\n").split("\t")
            if len(parts) == 5 and parts[4] == "Z":
                counts["CG"] += 1
    with gzip.open(samples[sample], "rt") as fh:
        for line in fh:
            if line.startswith("Bismark methylation extractor version"):
                continue
            parts = line.rstrip("\n").split("\t")
            if len(parts) != 5:
                continue
            if parts[4] == "X":
                counts["CHG"] += 1
            elif parts[4] == "H":
                counts["CHH"] += 1
    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

最后应该得到哪些结果

主线完成后，至少应该得到：

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_context_*.txt.gz
06_methylation/Non_CpG_context_*.txt.gz
07_igv/*.bw
08_plots/*official_like_profile.norm.png
08_plots/*.context_pie.png

DNA 甲基化数据分析

教程信息

数据来源与参考

流程总览

登录服务器并进入计算节点

开始前要知道的文件

创建工作目录

创建分析目录

复制参考文件和原始 reads

复制课堂脚本

从 GFF3 提取 chr1 gene BED

FastQC 质控

Trim Galore

Bismark 原理简介

准备 Bismark 参考索引

Bismark 比对

去重复、排序、建索引

Bismark 甲基化提取

生成 IGV 轨道

IGV 查看

deepTools profile 可选流程

全部 chr1 gene 的基因区域 profile 图

统计 CG / CHG / CHH 并画扇形图

最后应该得到哪些结果

一起交流

评论与讨论

教程信息#

数据来源与参考#

流程总览#

登录服务器并进入计算节点#

开始前要知道的文件#

创建工作目录#

创建分析目录#

复制参考文件和原始 reads#

复制课堂脚本#

从 GFF3 提取 chr1 gene BED#

FastQC 质控#

Trim Galore#

Bismark 原理简介#

准备 Bismark 参考索引#

Bismark 比对#

去重复、排序、建索引#

Bismark 甲基化提取#

生成 IGV 轨道#

IGV 查看#

deepTools profile 可选流程#

全部 chr1 gene 的基因区域 profile 图#

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

最后应该得到哪些结果#

一起交流

评论与讨论

登录后参与讨论

教程信息

数据来源与参考

流程总览

登录服务器并进入计算节点

开始前要知道的文件

创建工作目录

创建分析目录

复制参考文件和原始 reads

复制课堂脚本

从 GFF3 提取 chr1 gene BED

FastQC 质控

Trim Galore

Bismark 原理简介

准备 Bismark 参考索引

Bismark 比对

去重复、排序、建索引

Bismark 甲基化提取

生成 IGV 轨道

IGV 查看

deepTools profile 可选流程

全部 chr1 gene 的基因区域 profile 图

统计 CG / CHG / CHH 并画扇形图

最后应该得到哪些结果