当前位置：首页 > news >正文

生信软件49 - 全基因组亚硫酸氢盐测序（WGBS）比对与甲基化水平调用工具BSseeker2

news 2025/8/2 22:40:43

全基因组亚硫酸氢盐测序（WGBS）比对与甲基化水平调用工具

1. BS-Seeker2简介

BS-Seeker2是一个多功能的管道，用于准确和快速比对亚硫酸氢盐处理的reads。BS-Seeker 2的下游分析建议用CGmapTools (https://cgmaptools.github.io/)。支持的库类型: 全基因组亚硫酸氢盐测序（WGBS）与还原代表性亚硫酸氢盐测序（RRBS）。要求Python（版本2.6 +），pysam软件包（版本0.6.x）。支持的对齐工具： bowtie, bowtie2, soap。

github 地址： https://github.com/BSSeeker/BSseeker2

# 安装
# 下载最新版
wget -c https://github.com/BSSeeker/BSseeker2/archive/refs/tags/BSseeker2-v2.1.8.tar.gz
tar zxvf BSseeker2-v2.1.8.tar.gz# 常见python2 conda环境
conda create -n BSseeker2 python=2.7.18
conda activate BSseeker2

2. Reads过滤模块： FilterReads.py

可选和独立的模块。一些读段在PCR期间将被过度扩增。这个脚本进行比对之前获得唯一比对的reads，在进行比对之前决定是否过滤读取。

python FilterReads.py \-i sample.fq \ # 输入qseq/fastq/fasta/sequence文件-o sample.filter.fq

3. 基因组索引构建模块：bs_seeker2-build.py

为BS-Seeker 2建立索引的模块。

# 使用bowite2为WGBS创建基因组索引
python bs_seeker2-build.py \-f genome.fa \--aligner=bowtie2 \-p ~/path/bowtie2/ \  # 指定bowtie2工具路径# 使用Bowtie2为WGBS文库创建基因组索引，片段长度[40bp, 400bp]
python bs_seeker2-build.py \-f genome.fa \-l 40 \-u 400 \--aligner=bowtie2

4. 比对模块：bs_seeker2-align.py

用于在3-字母转化的基因组上映射读段的模块。

# Usage: bs_seeker2-align.py {-i <single> | -1 <mate1> -2 <mate2>} -g <genome.fa> [options]# WGBS 文库 ; 比对到基因组
python bs_seeker2-align.py \-i WGBS.fa \--aligner=bowtie2 \-o WGBS.bam \-f bam \-g genome.fa# 
# WGBS 文库 ; 比对到基因组, 双端PE模式
python bs_seeker2-align.py \-i WGBS.fa \-m 3 \--aligner=bowtie2 \-o WGBS.bam \-f bam \-g genome.fa \--bt2--end-to-end

4.1 BS_Seeker格式

read10	 1	+FW	chr1+0000169137	TC_CGGGGGTTATATGAGTGTGACGGCTGTAGCGTTAGGTGACGATGTCATCTCCGCGTTCCAAGCGTTATGTGCGCACTGAGGGACACATCCACGTTCCCGG_GG	CGGGGGTTATATGAGTGTGATGGTTGTAGCGTTAGGTGATGATGTTATTTTTGCGTTTTAAGCGTTATGTGCGTATTGAGGGATATATTTACGTTTTTGA	X-------------------x--y-----X---------x-----z--z-yx-X---zz---X--------X-z-y-------z-z--zz-X---zyx--	0	77	169135	169235
read102	 1	+FW	chr1+0000169137	TC_CGGGGGTTATATGAGTGTGACGGCTGTAGCGTTAGGTGACGATGTCATCTCCGCGTTCCAAGCGTTATGTGCGCACTGAGGGACACATCCACGTTCCCGG_GG	CGGGGGTTATATGAGTGTGATGGTTGTAGCGTTAGGTGATGATGTTATTTTTGCGTTTTAAGCGTTATGTGCGTATTGAGGGATATATTTACGTTTTTGA	X-------------------x--y-----X---------x-----z--z-yx-X---zz---X--------X-z-y-------z-z--zz-X---zyx--	0	77	169135	169235
read104	 0	+FW	chr1+0000325341	-C_CGGCAAACACCACGCCCCGCGATATGGCAGGATTCATGCCGACTAATGGAAAACACACCAGATGCTGGAAAGAGATAAAGGAGAGCGTTACTGCAATACT_GT	CGGTAAATATTACGTTTCGCGATATGGTAGGATTTATGTCGATTAATGGAAAATATATCAGATGTTGGAAAGAGATAAAGGAGAGCGTTATTGTAATATT	X--z---z-zz-X-zzyX-X-------y------z---yX--z----------z-z-zY-----y--------------------X----y--z----y-	0	154	325339	325509
read105	 0	+FW	chr1+0000238994	-C_CGGCCACACAGTGAAAGGCTGGGCTGTGAGAGCTTCGGTGGAAACCAGGCCTTCACCACTTCTTCTCCCTTCAAGCCACACACAGCTGTTGCAAGTTCCG_G-	CGGTTACATAGTGAAAGGTTGGGTTGTGAGAGTTTTGGTGGAAATTAGGTTTTTATTATTTTTTTTTTTTTTAAGTTATATATAGTTGTTGTAAGTTTCG	X--zz-Z-y---------y----y--------z--x--------zy---zz--z-zz-z--z--z-zzz--z---zz-z-z-y--y-----z-----yX-	0	118	238992	239093

 # 格式说明
（1） 读取ID（来自seq/fastq/qseq/fasta文件中的标题列，或原始输入的序列号）
（2） 第6列和第7列中基因组序列与BS读取列表之间的不匹配数量。不包括读取Ts到基因组Cs之间的亚硫酸氢盐转化位点。
（3） 读取的链可能来自（+FW、+RC、-RC、-FW）
（4） 映射位置的坐标，表示[染色体]、[映射的链（“+”或“-”）]和[沃森链上映射的基因组序列的0基5'端坐标]。
（5） BS读取序列从5'到3'：如果读取被唯一映射为FW读取，则显示原始读取。如果读取与RC读取一样被唯一映射，则会显示它们的反向补码。（6） 甲基化位点的序列总结：甲基化的CG/CHG/CHH位点标记为X/Y/Z（大写），而非甲基化的CG/CHG/CSH位点标记为X/Y/Z（小写）。本栏直接总结自第6栏和第7栏。
（7） XS标签，当读取被识别为亚硫酸氢盐处理未完全转化时为1，否则为0
（8） my_region_serial，仅用于RRBS的标签，映射片段的序列id
（9） my_region_start，仅用于RRBS的标记，映射片段的开始位置
（10） my_region_end，仅用于RRBS的标记，映射片段的结束位置

5. 甲基化水平调用模块：bs_seeker2-call_methylation.py

该模块从映射结果调用甲基化水平。

 # 对于WGBS测序python bs_seeker2-call_methylation.py \-i WGBS.bam \-o output \--db /path/bowtie2/ # bowtie2索引所在目录

5.1 输出文件

wig 文件

WIG文件格式。第2列的负值表示负链上的胞嘧啶。

variableStep chrom=chr1
3000419	0.000000
3000423	-0.2
3000440	0.000000
3000588	0.5
3000593	-0.000000

CG map文件

chr1	G	3000851	CHH	CC	0.1	1	10
chr1	C	3001624	CHG	CA	0.0	0	9
chr1	C	3001631	CG	CG	1.0	5	5

# 格式说明
(1)染色体
（2）沃森（+）链上的核苷酸
（3）位置
（4）背景（CG/CHG/CHH）
（5）二核苷酸背景（CA/CC/CG/CT）
（6）甲基化水平= C#/（C#+ T #）
（7）#_of_C（甲基化C，此处显示C的读段计数）
（8）= #_of_C + #_of_T（所有胞嘧啶，此处显示C或T的读数计数）

ATCG map文件

chr1	T	3009410	--	--	0	10	0	0	0	0	0	0	0	0	na
chr1	C	3009411	CHH	CC	0	10	0	0	0	0	0	0	0	0	0.0
chr1	C	3009412	CHG	CC	0	10	0	0	0	0	0	0	0	0	0.0
chr1	C	3009413	CG	CG	0	10	50	0	0	0	0	0	0	0	0.83

（1） 染色体
（2） Watson（+）链上的核苷酸
（3） 职位
（4） 上下文（CG/CHG/CHH）
（5） 二核苷酸上下文（CA/CC/CG/CT）（6） 此处映射的Watson链的读取次数，Watson链上支持A
（7） 此处映射的Watson链的读取次数，Watson链上支持T
（8） 这里映射的Watson链的读取次数，支持Watson链上的C
（9） 此处映射的Watson链的读取次数，Watson链上支持G
（10） 此处映射的Watson链的读取次数，支持N（11） 此处映射的Crick链的读取次数，Watson链上支持A，Crick链上支持T
（12） 此处映射的Crick链的读取次数，Watson链上支持T，Crick链上支持A
（13） 此处映射的Crick链的读取次数，Watson链支持C，Crick链支持G
（14） 此处映射的Crick链的读取次数，Watson链上支持G，Crick链上支持C
（15） 此处映射的Crick链的读取次数，支持N
（16） 对于Watson链，甲基化水平=#C/（#C+#T）=C8/（C7+C8），对于Crick链，=C14/（C11+C14）；
“nan”表示在这个位置没有支持C/T的读数。