BS-Seq(3) -- 数据比对并提取甲基化位点
2023-03-19 本文已影响0人
Z_bioinfo
这里我们使用的软件是bismark ，但是需要调用bowtie2，所以先将这两个软件下载好，解压运行。
构建甲基化比对基因组文件
#软件安装
conda install -c bioconda bismark bowtie2
#基因组
wget https://hgdownload.cse.ucsc.edu/goldenpath/hg19/bigZips/hg19.fa.gz
gunzip hg19.fa.gz
bismark_genome_preparation /home/jimmy/my_data/ref_data/hg19
运行成功后，会在fasta 所在目录生成如下的目录结构
├── Bisulfite_Genome
│ ├── CT_conversion
│ │ ├── BS_CT.1.bt2
│ │ ├── BS_CT.2.bt2
│ │ ├── BS_CT.3.bt2
│ │ ├── BS_CT.4.bt2
│ │ ├── BS_CT.rev.1.bt2
│ │ ├── BS_CT.rev.2.bt2
│ │ └── genome_mfa.CT_conversion.fa
│ └── GA_conversion
│ ├── BS_GA.1.bt2
│ ├── BS_GA.2.bt2
│ ├── BS_GA.3.bt2
│ ├── BS_GA.4.bt2
│ ├── BS_GA.rev.1.bt2
│ ├── BS_GA.rev.2.bt2
│ └── genome_mfa.GA_conversion.fa
可以看到，分别对CT和GA转换的基因组建立了bowtie2的索引。
比对
# paired-end data
bismark --genome /home/jimmy/my_data/ref_data/hg19 -1 read1.fastq.gz -2 read2.fastq.gz -p 4 -o ./ 2>test.log
# single-end data
bismark --genome /home/jimmy/my_data/ref_data/hg19 test_data.fastq -p 4 -o ./ 2>test.log
# 输出
├── test_data_bismark_bt2.bam
└── test_data_bismark_bt2_SE_report.txt
查看BAM文件 samtools view test_data_bismark_bt2.bam | head -n5 ：
SRR020138.15024317_SALK_2029:7:100:1672:902_length=86 16 chr1 57798677 42 50M * 0 0 TTCTTTCCCATCCCATAAATCCTAAAAATAATAAAAAATCATCCCCAAAT @@:AC@<=+@?+8)@BCCCA=6BCCCCCCCCCCCCCCCCACB=<88BCCA NM:i:11 MD:Z:14G2G6G0G0G0G4G1G1G0G10G1 XM:Z:..............z..h......hhhh....h.h.hh..........h. XR:Z:CT XG:Z:GA
SRR020138.15024318_SALK_2029:7:100:1672:137_length=86 0 chr12 129774096 8 50M * 0 0 AAAAAAAAAAAAAAGAAAAAAAAGAAAAAGAAAAGGAAAAGTAAAAAAAA =@CAA=@B@CB=98%:AB?>@56/=3<=<)>B@:*=:=61%,<A@@1+12 NM:i:2 MD:Z:41C5G2 XM:Z:.........................................h........ XR:Z:CT XG:Z:CT
SRR020138.15024319_SALK_2029:7:100:1672:31_length=86 0 chr2 10166575 42 50M * 0 0 ATTTTGTTATAGAGTGGGGTATTTTCGGGAAGAAGGAGGAGGAGTGTATT BCCCCBCCCCA?:=ACCBCABCCCCCBCCA??5=9@4BB@;??B@BABBA NM:i:8 MD:Z:1C1C5C9C2C0C22C1C1 XM:Z:.h.x.....x.........h..hh.Z....................h.x. XR:Z:CT XG:Z:CT
SRR020138.15024320_SALK_2029:7:100:1672:1164_length=86 16 chr5 28344472 8 50M * 0 0 CACAAAATATCAACACCCCTAAACCCCACATTATTCAAAAATCAATTATA @@@BBBA@A9=A@<?::2:<CB@?=:BBAC??CB@@BBBBC>:ACABCAB NM:i:11 MD:Z:4G1G1G3G9G9G5G0G0G1T4G2 XM:Z:....x.h.h...x.........h.........h.....hhh......h.. XR:Z:CT XG:Z:GA
SRR020138.15024321_SALK_2029:7:100:1672:433_length=86 0 chr14 38711099 42 50M * 0 0 TTTTGAGTAGAGAAGTTAGTATTTTAGGGAATTTTTGATTTTTTTAAGTT BCCBB?B@@A>@-4BBB:7@BBBCBBC@@=A@BCACA;BCBBCBB@@@BB NM:i:14 MD:Z:0C0C13C0C6C0C6C0C1C3C0C1C0C1C5 XM:Z:hh.............hx......hx......hh.x...hh.hh.h..... XR:Z:CT XG:Z:CT
甲基化call字符串对于BS-read中不涉及胞嘧啶的每个位置都用一个点 . 代替，或包含以下不同的胞嘧啶甲基化的字母 (大写=甲基化，小写=未甲基化) ：
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)
. # 该位置不是胞嘧啶
去除重复
deduplicate_bismark -s --bam test_data_bismark_bt2.bam
# 输出
├── test_data_bismark_bt2.deduplicated.bam
└── test_data_bismark_bt2.deduplication_report.txt
#-p paired-end
# -s single-end
甲基化位点提取
默认情况下，软件会自动根据 甲基化的C的类型 (CpG, CHG, CHH) 和 比对到四条链上 (OT, OB, CTOT, CTOB) 两个因素生成结果文件。
OT – original top strand
CTOT – complementary to original top strand
OB – original bottom strand
CTOB – complementary to original bottom strand
bismark_methylation_extractor -s --gzip --bedGraph --buffer_size 10G --cytosine_report --comprehensive --genome_folder /home/jimmy/my_data/ref_data/hg19 test_data_bismark_bt2.bam 2>extracor.log
#参数：
#-p paired-end
# -s single-end
#--comprehensive 将可能的甲基化信息都输出，包括CHG,CHH,CpG
#--no_overlap：对于双端读取，read_1和read_2理论上是可能重叠的。这个选项避免了重复的甲基化计算了2次。虽然这消除了对序列片段中心更多甲基化的计算偏差。
#--bedGraph:输出bedGraph文件
#--cytosine_report：输出全基因组所有的CpG。
#--counts 计算bedGraph中有每个C上甲基化reads和非甲基化reads的数目
#--buffer_size 指定内存
#--report ：产生甲基化的summary 文件
#--genome_folder:参考基因组的路径
#-o:输出文件路径
生成处理报告
bismark2report --dir ./ --alignment_report test_data_bismark_bt2_SE_report.txt
--dir指定结果的输出目录
--alignment_report 指定比对产生的report文件的路径
所有结果文件
CHG_context_test_data_bismark_bt2.deduplicated.txt.gz
CHH_context_test_data_bismark_bt2.deduplicated.txt.gz
CpG_context_test_data_bismark_bt2.deduplicated.txt.gz
test_data_bismark_bt2.deduplicated.bedGraph.gz
test_data_bismark_bt2.deduplicated.bismark.cov.gz
test_data_bismark_bt2.deduplicated.CpG_report.txt.CpG_report.txt.gz
test_data_bismark_bt2.deduplicated_splitting_report.txt
test_data_bismark_bt2.deduplicated.cytosine_context_summary.txt
test_data_bismark_bt2.deduplicated.M-bias.txt
SW480-WGBS_bismark_bt2.deduplicated_splitting_report.html
结果解读
--cytosine_report参数会根据当前目录下的信息文件生成一个HTML格式的报告文件，即test_data_bismark_bt2_SE_report.html文件，它包括了比对信息，甲基化信息，M-bias等，可以对数据有一个大概的认知（下图只展示了一部分）：
image.png
image.png
同时因为使用了--comprehensive，所以结果合并正反链的数据后会输出CpG/CHG/CHH三种类型的甲基化文件，包含了胞嘧啶所有的组合形式，但实际上我们自然最关注的是CpG位点的甲基化。其中
CpG_context_test_data_bismark_bt2.deduplicated.txt即CpG甲基化位点的文件。
head CpG_context_test_data_bismark_bt2.deduplicated.txt
Bismark methylation extractor version v0.19.0
SRR020138.15024317_SALK_2029:7:100:1672:902_length=86 - 1 57333019 z
SRR020138.15024319_SALK_2029:7:100:1672:31_length=86 + 2 10026473 Z
SRR020138.15024331_SALK_2029:7:100:1673:1282_length=86 + 11 78025243 Z
SRR020138.15024343_SALK_2029:7:100:1673:202_length=86 + 10 121617231 Z
SRR020138.15024357_SALK_2029:7:100:1673:879_length=86 - 4 75173715 z
SRR020138.15024361_SALK_2029:7:100:1673:235_length=86 - 2 130768889 z
SRR020138.15024368_SALK_2029:7:100:1673:123_length=86 + 10 106402850 Z
SRR020138.15024376_SALK_2029:7:100:1673:1908_length=86 - 12 124097382 z
SRR020138.15024380_SALK_2029:7:100:1673:397_length=86 + 8 95038627 Z
# 第一列为测序信息
# 第二列为甲基化状态 + 代表甲基化 -代表未甲基化
# 第三列代表chromosome
# 第四列代表location
# 第五列代表methylation call，简单来说大写的就是甲基化的（因为还有CHG，CHH的数据，分别对应x, X , h, H）
test_data_bismark_bt2.deduplicated.bismark.cov文件则给了每个位点的甲基化比例，为下一步确定CpG岛提供了基础，其数据形式如下：
$head test_data_bismark_bt2.deduplicated.bismark.cov
1 975476 975476 100 1 0
1 975488 975488 100 1 0
1 975490 975490 100 1 0
1 2224487 2224487 100 1 0
1 2224489 2224489 100 1 0
1 2224514 2224514 100 1 0
1 2224520 2224520 100 1 0
1 2313220 2313220 0 0 1
1 9313902 9313902 100 1 0
1 9313914 9313914 100 1 0
# 第一列代表chromosome
# 第二，三列代表location
# 第四列代表甲基化百分比
# 第五列代表甲基化数目
# 第六列代表未甲基化数目
test_data_bismark_bt2.deduplicated.CpG_report.txt.CpG_report.txt文件则是背景信息：
$head test_data_bismark_bt2.deduplicated.CpG_report.txt.CpG_report.txt
1 10469 + 0 0 CG CGC
1 10470 - 0 0 CG CGA
1 10471 + 0 0 CG CGG
1 10472 - 0 0 CG CGC
1 10484 + 0 0 CG CGG
1 10485 - 0 0 CG CGG
1 10489 + 0 0 CG CGC
1 10490 - 0 0 CG CGG
1 10493 + 0 0 CG CGC
1 10494 - 0 0 CG CGG
# 第一列为chromosome
# 第二列为location
# 第三列为strand
# 第四，五列为甲基化和非甲基化的碱基数目
# 第六列为CG
# 第七列为背景（第一个C延伸两个碱基）
其它参数
bam2nuc
bismark2summary
coverage2cytosine
NOMe_filtering
filter_non_conversion
# 有需要的可以自信探索 --help or manual
此处根据测序数据得到了甲基化位点的信息，但是后续DML以及DMR的确定还需要R包的使用，以及后续的可视化还以探索以下包：
image.png
Methy-Pipe: An Integrated Bioinformatics Pipeline for Whole Genome Bisulfite Sequencing Data Analysis
2014年出的一个pipeline，分析作图一条龙，有兴趣的同学安排一下。
