WGBS上游分析

WGBS的分析全流程：

主要参考资料：
WGBS甲基化分析
Bismark软件使用入门
沉浸式体验WGBS(上游)
甲基化流程浅析
听说你不会处理WGBS数据？安排上
全基因组甲基化分析简述：使用BS-Seeker2
Bismark Bisulfite Mapper学习笔记（二）甲基化信息提取以及文件解读
DNA甲基化分析流程详解

查看用户quota：quota -uvs

数据清洗

Trim galore可以自动检测adapter对NGS数据进行质量过滤，还可以查看质量分布图。也就是等于fastqc+trimmomatic。trim_galore适用于所有高通量测序，包括RRBS(Reduced Representation Bisulfite-Seq ), Illumina、Nextera和smallRNA测序平台的双端和单端数据
Trim_galore的工作流程：
第一步 首先去除低质量碱基，然后去除3'末端的adapter(如果没有指定具体的adapter，程序会自动检测前1million的序列）
第二步 对比前12-13bp的序列是否符合以下类型的adapter:

• Illumina: AGATCGGAAGAGC # 如果输入参数 --Illumina,就会默认trimmed前13bp的adapter
• Small RNA: TGGAATTCTCGG # 同上 如果输入参数 --small_rna,就会默认trimmed前12bp的adapter
• Nextera: CTGTCTCTTATA # 同上 如果输入参数 --nextera,就会默认trimmed前12bp的adapter

在trim_galore运行中，会先生成trimmed文件，随后生成val文件，可以看官方讲解这二者的区别：Can you explain the difference between _val_1/_val_2 and _trimmed?

# 安装Trim galore
conda install -c bioconda trim-galore -y

fastqc --noextract -t 10 -f fastq -o ./fastqc ./rawdata/*_1.fastq.gz
fastqc --noextract -t 10 -f fastq -o ./fastqc ./rawdata/*_2.fastq.gz

# 双端测序数据
cleandata=./1_trim_galore
rawdata=./0_rawdata
cat sample.txt | while read id
do
echo "trim_galore --fastqc --paired -o ./1_trim_galore ./0_rawdata/${id}_1.fastq.gz ./0_rawdata/${id}_2.fastq.gz" 
done > trim_galore.sh

参数详解：

• --phred33：指定测序数据的Phred质量值，33或者64，默认为33。如何判断呢，可以见下面两个文章，但是最近的数据基本上都是33了：

fastq格式文件及phred33的判断
Fastq 格式说明 & (Phred33 or Phred64)

• -j 8：指定线程数，8为最大。它的使用基于python3版本，如果报错，可以尝试更新python环境
• -q 20：指定质量值，切除质量得分低于20的序列
• --length：设定输出reads长度阈值小于设定值会被抛弃，默认值20bp；即小于20bp的被去除。注意，在pe150下，不要设置太高，可以50或36
• --stringency 3：可以忍受的前后adapter重叠的碱基数，默认为1（非常苛刻）。可以适度放宽，因为后一个adapter几乎不可能被测序仪读到
• --fastqc：当分析结束后，使用默认选项对结果文件进行fastqc分析
• --paired：对于双端测序结果，一对reads中，如果有一个被剔除，那么另一个会被同样抛弃，而不管是否达到标准
• --max_n 3：read中允许包含的N的总数（整数），否则将完全删除该read。 在配对末端设置中，任一read超过此限制将导致整个配对从修剪后的输出文件中删除。
• -o：指定输出目录

数据质控

# 安装MultiQc
conda install -c bioconda multiqc -y
# 整合FastQC结果
multiqc *.zip

WGBS比对

准备参考基因组

# 下载NCBI上猪的参考基因组
wget https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/003/025/GCF_000003025.6_Sscrofa11.1/GCF_000003025.6_Sscrofa11.1_genomic.fna.gz
# 解压并重命名
gunzip GCF_000003025.6_Sscrofa11.1_genomic.fna.gz && mv GCF_000003025.6_Sscrofa11.1_genomic.fna pig.fa

# 下载Ensembl上猪的参考基因组
wget ftp://ftp.ensembl.org/pub/release-102/fasta/sus_scrofa/dna/Sus_scrofa.Sscrofa11.1.dna.toplevel.fa.gz

# 解压并重命名
gunzip Sus_scrofa.Sscrofa11.1.dna.toplevel.fa.gz && mv Sus_scrofa.Sscrofa11.1.dna.toplevel.fa pig.fa

探索参考基因组是否有NW_018085246.1

# 查看参考基因组的染色体编号
grep ">" pig.fa | cut -d " " -f 1 | cut -d ">" -f 2 | sort | uniq -c > pig.fa.chr.txt

# 查看包含NW_018085246对应的行数
grep -n NW_018085246 pig.fa
# 31088069:>NW_018085246.1 Sus scrofa isolate TJ Tabasco breed Duroc unplaced genomic scaffold, Sscrofa11.1 Contig556, whole genome shotgun sequence

# 输出31088069-31088070行
sed -n '31088069,31088070p' pig.fa

转化参考基因组，并建立索引

在比对之前，需要加入lambda序列数据，lambda序列是在重硫酸盐实验前人为加入的一段序列用来评估重硫酸盐的转化效率，这个一般由为你测序的公司提供，所以我们这一步就将其合并到参考基因组中

看到了有这种说法，但是我暂时并不知道这个lambda序列

# 安装bismark和bowtie2
conda install -c bioconda bismark bowtie2 -y

# 查看bowtie2的路径
which bowtie2

# 转化参考基因组
nohup bismark_genome_preparation \
    --bowtie2 \
    --parallel 10 \
    /data/wangji/wgbs/genome > bismark_genome_preparation.log 2>&1 &

重要参数说明：
--bowtie2：调用bowtie2/hisat2建立基因组索引
--path_to_aligner：比对软件所在文件夹的全路径
--verbose：输出详情
--parallel：设置线程，索引建立是并行运行，因此实际线程要×2
--large-index：大基因组索引建立
--yes：如果有安全类问题则自动选择yes，比如覆盖某个已存在的文件
<path_to_genome_folder>：基因组所在文件夹路径，即/data2/zhangzipeng/wang/wgbs/genome

比对

比对的中间文件详解见：bismark 识别甲基化位点-比对篇

cat sample.txt | while read id
do
echo "bismark --bowtie2 -p 8 -o ./2_bamfile --genome /work/home/wklgroup01/wgbs/genome -1 ./1_trim_galore/${id}_1_val_1.fq.gz -2 ./1_trim_galore/${id}_2_val_2.fq.gz &"
done > bismark.sh

nohup sh bismark.sh &

bismark -N 0 -L 20 \
    --un --ambiguous --sam \
    --bowtie2 \
    --path_to_bowtie2 ~/miniconda3/envs/wgbs/bin/ \
    --parallel 10 \
    -o ./bamfile \
    --fastq \
    --prefix SRR6457892 \
    --genome /data/wangji/wgbs/smk/pig_genome \
    -1 ./SRR6457892_1_val_1.fq.gz \
    -2 /data2/zhangzipeng/wang/wgbs/cleandata/SRR6457892_2_val_2.fq.gz

nohup bismark --bowtie2
--multicore 20
--gzip
-p 10
-o bamfile_ensembl
--genome ./genome/ncbi
-1 ./trim_reads/SRR5582006_1_val_1.fq.gz
-2 ./trim_reads/SRR5582006_2_val_2.fq.gz > bismark_ensembl.log 2>&1 &
重要参数说明：

--multicore ：这个决定了以多少个线程同时运行，运行的时候会产生相同数量的temp文件，temp文件会有正常的和C_to_T转换的
--gzip：输出文件压缩
-p:这个参数是bowtie2的参数，指定线程数

并行的一些方法：

cat srr | while read id
do
echo "bismark --bowtie2 --path_to_bowtie2 ~/miniconda3/envs/wgbs/bin/ -o ./mapping --fastq --prefix ${id} --genome /data2/zhangzipeng/wang/wgbs/genome -1 ${id}_1.fastq.gz -2 ${id}_2.fastq.gz"
done > bismark.sh

cat srr | while read id; do
  (bismark --bowtie2 --path_to_bowtie2 ~/miniconda3/envs/wgbs/bin/ -o ./mapping --fastq --prefix ${id} --genome /data2/zhangzipeng/wang/wgbs/genome -1 ${id}_1.fastq.gz -2 ${id}_2.fastq.gz) &
done
wait

查看bam文件

# 查看bam文件
samtools view -h SRR5582006_1_val_1_bismark_bt2_pe.bam | head -n 500
# 能够看到有@SQ	SN:NW_018085246.1	LN:362670
# 其中SN表示染色体编号，LN表示染色体长度

输出文件：

(a) test.file.R1_bismark_bt2_pe.sam 所有对齐和甲基化的信息，也即储存了成功比对的reads信息

甲基化调用字符串包含一个点“.”代表 BS-read 中不涉及胞嘧啶的每个位置，或者包含以下三个不同胞嘧啶甲基化上下文的字母之一（大写 = 甲基化，小写 = 未甲基化）：

字母	含义
z	unmethylated C in CpG context
Z	methylated C in CpG context
x	unmethylated C in CHG context
X	methylated C in CHG context
h	unmethylated C in CHH context
H	methylated C in CHH context
u	unmethylated C in Unknown context (CN or CHN)
U	methylated C in Unknown context (CN or CHN)
对于CpG, 采用字母Z的大小写来表征甲基化状态；对于CHG, 采用字母X的大小写来表征甲基化状态；对于CHH， 采用字母H 的大小写来表征甲基化状态。

(b) test.file.R1_bismark_bt2_PE_report.txt 对齐和甲基化的主要信息概括

生成报告

运行bismark2report生成Testpaired_PE_report.html报告

bismark2report

比对后的数据处理

去除重复序列

RRBS一定不要做这一步

cat sample.txt | while read id
do
echo "deduplicate_bismark --bam -p ./2_bamfile/${id}_1_val_1_bismark_bt2_pe.bam --output_dir ./3_duplicate/ &"
done > bismark.sh

nohup deduplicate_bismark -p --bam SRR5582006_1_val_1_bismark_bt2_pe.bam --output_dir ./deduplicate/ > deduplicate.log 2>&1 &

重要参数说明：
--sam/--bam：删除 Bismark 对齐产生的SAM/BAM 文件中的重复数据
-p/--paired：前一步双端数据产生的结果文件
-s/--single：前一步单端数据产生的结果文件
--samtools_path：samtools所在文件夹的全路径
--output_dir：输出文件夹路径
--multiple：指定输入文件都作为一个样本处理，连接在一起进行重复数据删除

得到输出文件：
(a) test.file.R1_bismark_bt2_pe.deduplicated.bam
(b) test.file.R1_bismark_bt2_pe.deduplication_report.txt：去重复后的结果描述文件

甲基化信息提取

运行命令

# 后续的plot文件需要Perl模块GD::Graph安装
# conda install -c bioconda perl-gdgraph -y
bismark_methylation_extractor -p --gzip --no_overlap --comprehensive --bedGraph --count --CX_context --cytosine_report --buffer_size 20G --genome_folder /work/home/wklgroup01/wgbs/genome/ -o testfile_report ~/bismark_example/03duplicate/test.file.R1_bismark_bt2_pe.deduplicated.bam &

bismark_methylation_extractor --pair-end --comprehensive --no-overlap --bedGraph --counts --report --cytosine_report --genome_folder /work/home/wklgroup01/wgbs/genome ../3_duplicate/*.bam -o ../4_methylation_result

ls *.sam | awk '{print "\"" $0 "\""}' | paste -sd ","
# Define an array of input files
files=(
"SRR25022243_1_val_1_bismark_bt2_pe.deduplicated.sam"
"SRR15132626_1_val_1_bismark_bt2_pe.deduplicated.sam"
)

# Get the input file corresponding to the current job array index
input_file=${files[$SLURM_ARRAY_TASK_ID-1]}

# Run the bismark_methylation_extractor command for the input file
bismark_methylation_extractor \
    -p --gzip --no_overlap \
    --comprehensive \
    --bedGraph \
    --count \
    --CX_context \
    --cytosine_report \
    --buffer_size 20G \
    --genome_folder /work/home/wklgroup01/wgbs/genome/ \
    -o testfile_report \
    "/work/home/wklgroup01/wgbs/deduplicate/$input_file"

nohup bismark_methylation_extractor -p --ignore_r2 6 --multicore 10 --bedgraph -o ./methylation --cytosine_report --genome_folder ./genome/ncbi ./bamfile/*.deduplicated.bam > bismark_methylation_extractor.log 2>&1 &
重要参数说明：
--comprehensive ：合并所有四个可能的特定链，将甲基化信息转换为context-dependent的输出文件
--no_overlap：避免因双端读取read_1和read_2理论上的重叠，导致甲基化重复计算。可能会删去相当大部分的数据，对于双端数据的处理，默认情况下此选项处于启用状态，可以使用--include_overlap禁用。
-p/--paired-end：前一步双端数据产生的结果文件
--bedGraph：指将产生一个BedGraph文件存储CpG的甲基化信息（不与--CX联用，仅产生CpG信息）
--CX/--CX_context：与--bedGraph联用，产生一个包含三种类型甲基化信息的BedGraph文件;与--cytosine_report联用，产生一个包含三种类型甲基化信息的全基因组甲基化报告
--cytosine_report：指将产生一个存储CpG的甲基化信息的报告（不与--CX联用时，仅产生CpG信息），默认情况坐标系基于 1
--buffer_size：甲基化信息进行排序时指定内存，默认为2G
--zero_based：使用基于 0 的基因组坐标而不是基于 1 的坐标，默认值：OFF
--split_by_chromosome：分染色体输出
--genome_folder：包含未修改的参考基因组和bismark_genome_preparation创建的子文件夹(CT_conversion/和GA_conversion/）的文件夹的路径，即~/bismark_example/01index/
:BAM 格式的 Bismark 结果文件

输出文件解读：

1、CHG/CHH/CpG_context_test.file.R1_bismark_bt2_pe.deduplicated.txt

OT – original top strand原始top链
CTOT – complementary to original top strand原始top链的互补链
OB – original bottom strand原始bottom链
CTOB – complementary to original bottom strand原始bottom链互补链

来自OT和CTOT的甲基化calls将提供原始top链上胞嘧啶甲基化位置的信息，来自OB和CTOB的甲基化calls将提供原始bottom链上胞嘧啶甲基化位置的信息。请注意，在Bismark比对步骤中指定--directional(默认模式)选项不会对CTOT或CTOB链生成任何报告，所以往往只得到OT和OB的甲基化报告信息。

col1 : 比对上的序列ID
col2 : 基因组正负链：+ -
col3 : 染色体编号
col4 : 染色体位置
col5 : 甲基化C的状态（XxHhZzUu）

X 代表CHG中甲基化的C
x  代笔CHG中非甲基化的C
H 代表CHH中甲基化的C
h  代表CHH中非甲基化的C
Z  代表CpG中甲基化的C
z  代表CpG中非甲基化的C
U 代表其他情况的甲基化C(CN或者CHN)
u  代表其他情况的非甲基化C (CN或者CHN)

CpG：甲基化C下游是个G碱基
CHH：甲基化C下游的2个碱基都是H（A、C、T）
CHG：甲基化的C下游的2个碱基是H和G

# 输出NW_018085246的甲基化信息
grep -n NW_018085246 CHG_OB_SRR5582006_1_val_1_bismark_bt2_pe.deduplicated.txt > NW_018085246.txt
# 18395:SRR5582006.5464_FCC6MP9ACXX:7:1101:20571:36789#/1	-	NW_018085246.1	42109	x
# 18404:SRR5582006.5464_FCC6MP9ACXX:7:1101:20571:36789#/1	+	NW_018085246.1	42081	X

# 查看18395-18404行
sed -n '18395,18404p' CHG_OB_SRR5582006_1_val_1_bismark_bt2_pe.deduplicated.txt

2、test.file.R1_bismark_bt2_pe.deduplicated.bedGraph.gz

# 使用zcat进行查看
zcat test.file.R1_bismark_bt2_pe.deduplicated.bedGraph.gz | head -n 10

col1 : 染色体编号
col2 : 胞嘧啶（c）起始位置信息
col3 : 胞嘧啶（c）终止位置信息
col4 : 甲基化率

3、test.file.R1_bismark_bt2_pe.deduplicated.bismark.cov.gz

由于甲基化百分数本身并不能提供在某一位置检测到的甲基化或非甲基化read的实际read覆盖率，bismark2bedGraph还会输出一个coverage文件(使用基于1的基因组坐标)，该文件具有两个附加列:

• col5 : 甲基化数目
• col6 : 非甲基化数目

# 使用zcat进行查看
zcat test.file.R1_bismark_bt2_pe.deduplicated.bismark.cov.gz | head -n 10

# 使用解压后SRR5582006_1_val_1_bismark_bt2_pe.deduplicated.bismark.cov文件输出NW_018085246的甲基化信息
grep -n NW_018085246 SRR5582006_1_val_1_bismark_bt2_pe.deduplicated.bismark.cov > NW_018085246.cov.txt
# head NW_018085246.cov.txt
# 48814275:NW_018085246.1	749	749	0	0	1
# 48814276:NW_018085246.1	802	802	100	1	0
# 48814277:NW_018085246.1	821	821	100	1	0

# 简单处理文件
## 依据冒号分割，取第二列
cut -d ":" -f 2 NW_018085246.cov.txt > NW_018085246.cov.txt.1 && mv NW_018085246.cov.txt.1 NW_018085246.cov.txt
## 增加行名：Chromosome	Start	End	Feature	methy_level

4、test.file.R1_bismark_bt2_pe.deduplicated.CX_report.txt

默认情况下只考虑CpG背景中的胞嘧啶，但是可以通过使用--CX 扩展到任何序列背景中的胞嘧啶

col1 : 染色体编号
col2 : 位置
col3 : 正负链信息
col4 : 甲基化碱基数目
col5 : 非甲基化碱基数目
col6 : 类型
col7 : 三核苷酸背景

5、test.file.R1_bismark_bt2_pe.deduplicated.M-bias.txt/...R1.png/...R2.png

#reads中每个可能位置的甲基化比例
CpG context (R1)
col1 : read position
col2 : 甲基化count（count methylated）
col3 : 非甲基化count（count unmethylated）
col4 : 甲基化率（beta）
col5 : 总coverage

6、test.file.R1_bismark_bt2_pe.deduplicated_splitting_report.txt：甲基化总体报告

使用IGV查看比对结果

参考资料：

干货分享 | IGV软件的安装和常用操作介绍
常见设置
保姆级 IGV 基因组浏览器使用指南（图文详解）
使用bismark进行甲基化分析并使用IGV对结果进行可视化

有非常多种针对甲基化结果进行处理，以实现IGV可视化的方法，这主要是因为IGV支持多种文件格式，包括：

• 拷贝数文件：seg
• 序列比对文件：bam、cram
• 基因组注释文件；bed、gtf、gff3、psl、bigbed
• 定量文件：wig、bedgraph、bigwig、tdf
  要想自己生成对应的文件格式，可以参考此文档：File Formats

首先对IGV软件导入参考基因组及基因注释文件

构建索引:IGV导入参考基因组和基因注释文件都需要先构建索引

igv里面对参考基因组进行构建索引

IGV 工具栏，tools-Run igvtools；选择index
igv里面给注释文件排序，构建索引
IGV 工具栏，tools-Run igvtools；选择sort；输入注释文件；生成sort；接着，IGV 工具栏，tools-Run igvtools；选择index；输入刚刚的sort文件；生成index

导入文件

参考基因组

IGV 工具栏，Genomes → load genome from
或者Genomes → Create genome File

gff/gtf注释文件
File → Load from File→找到注释文件即可

IGV中gtf文件的各种样式含义：

对bam进行可视化

对bam文件进行可视化，需要先对bam文件进行排序，然后进行索引才可以导入IGV中，这个过程可以使用samtools进行，也可以使用IGV自带的工具进行

几乎所有的reads中都有很多红色和绿色的部分，这个是因为默认的颜色模式下，未发生甲基化的C在正义链和反义链中会被转化为T和A碱基，被误判为突变碱基。可以在BAM文件右键选择"Color alignments by"→ "bisulfite mode" → "CG"，转换为甲基化模式。这里考虑到大多数情况下关注CG的甲基化，因此选择了CG模式，如果关注其他甲基化，可以将"CG"替换为对应的甲基化模式。

将bedGraph文件转化为tdf格式并导入

甲基化提取的输出可以使用选项--bedGraph转换为一个bedGraph和cov文件。这个步骤也可以使用独立脚本bismark2bedGraph从甲基化提取输出中完成，会获得两个文件：bedGraph和cov
BedGraph文件记录各个位点的甲基化率。BedGraph文件会比较大，不太适合直接导入IGV中，可以：

• 按照染色体分别提取bedGraph

# 例如直接提取SRR15132626_1_val_1_bismark_bt2_pe.deduplicated.bedGraph文件中NW_018085246.1
grep "NW_018085246.1" SRR15132626_1_val_1_bismark_bt2_pe.deduplicated.bedGraph > NW_018085246.1.bedGraph
# 如果是gz文件，可以使用
zcat SRR15132626_1_val_1_bismark_bt2_pe.deduplicated.bedGraph.gz | grep "NW_018085246.1" > NW_018085246.1.bedGraph

# 让目录下的所有bedGraph.gz文件只保留NW_018085246.1染色体的信息
ls *.bedGraph.gz | while read id
do
echo "zcat ${id} | grep \"NW_018085246.1\" > NW_018085246.1.bedGraph.${id}"
done > NW_018085246.1.bedGraph.sh

# 让目录下的所有NW_018085246.1文件只保留第二列大于139027小于144527的行
ls *.NW_018085246.1.bedGraph | while read id
do
echo "awk '{if(\$2>139027 && \$2<144527) print \$0}' ${id} > 1${id} & rm ${id}"
done > NW_018085246.1.bedGraph.sh

- 使用igvtools转化为tdf格式
>菜单栏中选择"Tools" → "Run igvtools"，在弹出窗口中，"Command"选择"To TDF"，"Input File"选择需要转换的BedGraph文件，其他参数保持默认，点击"Run"可以进行转化

### 将bedGraph文件转化为bw格式
[滑动窗口展示甲基化程度](https://www.jianshu.com/p/4533d37a605e)
```bash
# 排序
awk '{ print $1"\t"$2"\t"$3"\tid-"NR"\t"$4; }' GSM1032058_1_val_1_bismark_bt2_pe.deduplicated.bedGraph | sort-bed - > input.bed

# 下载各个染色体的长度文件并排序
fetchChromSizes hg38 > hg38.bounds.unsorted.txt
awk '{ print $1"\t0\t"$2; }' hg38.bounds.unsorted.txt | sort-bed - > hg38.bounds.bed

# 获取最终的bed文件
bedops --chop 1000  hg38.bounds.bed | bedmap --faster --echo --mean --delim "\t" --skip-unmapped - input.bed > final_no_stagger.bed

# 转换为bw文件
bedGraphToBigWig final_no_stagger.bed hg38.bounds.unsorted.txt final.bw

最后将final_no_stagger.bed和 final.bw导入IGV查看

针对deduplicated.bismark.cov文件进行处理，获取格式

deduplicated.bismark.cov文件的第一列到第三列是甲基化位点的坐标，第四列是位点甲基化水平，第五列是methylated counts, 第六列是unmethylated counts。有了它，我们就可以知道全基因组CHG或CpG或CHH甲基化水平的分布

## 下载包 
install.packages("BiocManager")
library(BiocManager)
BiocManager::install("tidyverse")
library(tidyverse)
## 文件导入
SRR5582006 <- read.delim("SRR25022243_1_val_1_bismark_bt2_pe.deduplicated.bismark.cov", header=FALSE)
## 文件处理
SRR5582006<-SRR5582006[,1:4]
SRR5582006$Feature<-"SRR5582006"
SRR5582006<-select(SRR5582006, V1, V2, V3, Feature, V4)
SRR5582006$V4<-SRR5582006$V4/100
colnames(SRR5582006)<-c("Chromosome", "Start", "End", "Feature", "methy_level")
## 文件导出
write.table(SRR5582006, file="SRR5582006.igv", sep = "\t", row.names = F, quote = F)

加载得到的.igv文件：
“File”--"Load from File"--选择.igv文件。如果.igv文件超过240 MB会出现以下提示，但是我们忽略它，选择"Continue"

针对deduplicated.CX_report.txt文件进行处理，获取格式

bismark_methylation_extractor会生成一个总的bedgraph文件：*_P_1_bismark_bt2_pe.deduplicated.bedGraph.gz。但该文件是整个基因组所有甲基化C位点，没有分CG/CHG/CHH信息，也不利于后续可视化，因此根据下面代码可以获取的CX_report.txt（按照染色体拆分）文件拆分CG，CHG，CHH三类甲基化后转换成bedgraph格式文件。

nohup bismark_methylation_extractor -p --parallel 40 --comprehensive --no_overlap--bedGraph --cytosine_report --CX_context --split_by_chromosome --counts--buffer_size 20G --report --samtools_path yourpath/smrtlink_v9/smrtcmds/bin --genome_folder refdir bam    -o outdir &

CX_report.txt各列分别为：
col1：染色体；
col2：染色体上具体位点（以1开始）；
col3：正负链；
col4：比对到该位点发生甲基化的reads数；
col5：比对到该位点未发生甲基化的read数；
col6：甲基化类型，有三种：CG, CHG, CHH；
col7：该位点的三核苷酸序列；

CX_report.txt需要使用bismarkCXmethykit.pl脚本进行转换，转换后出现以CG_methykit.txt，CHG_methykit.txt，CHH_methykit.txt为后缀的拆分甲基化类型的文件，文件格式如下：

再使用R脚本转换成bedgraph文件，只需保留chr列，base列和freqC列，R脚本如下：

Args <- commandArgs(T)
Args[1]
myfile=read.table(Args[1],header=T)
outfile=sub("_P_1_bismark_bt2_pe.deduplicated.CX_report.txt.chrChr19.CX_report.txt", "",Args[1])
file_out=sub("_methykit.txt", ".bedgraph", outfile )   ###输出文件
chr<-as.character(myfile$chr)
pos<-as.integer(myfile$base)
pos2<-as.integer(myfile$base)
freq<-as.numeric(myfile$freqC)
outf<-data.frame(chr,pos,pos2,freq, check.names = F)   ###合并以上四列，check.names = F可以防止将字符串转化为其他内容。
outf$chr<-factor(outf$chr)
dim(myfile)
dim(outf)  ###检查一下新文件和原文件行数是否相等
write.table(outf, file=file_out, quote=F,col.names = F, row.names = F)  ##输出

转换后生成以CG.bedgraph/CHG.bedgraph/CHH.bedgraph为后缀的文件，文件内容如下：

符合IGV需求的格式，可以用于可视化