BS-seeker2 CGmaptools 处理DNA甲基化流程

--------------------1. BSseeker2----------------

安装BSseeker2

git clone git://github.com/BSSeeker/BSseeker2.git

1.1 构建针对BS-seq的index

1.BS-seeker2与Bismark的基本策略相似，都是将基因组和测序数据看作"三碱基"再进行比对；在序列的比对过程中调用bowtie或者bowtie2。
2.和bowtie/bowtie2相似，BS-Seeker2也需要首先建立一个index。但要注意的是，BS-Seeker2不能直接使用bowtie的index，而是需要创建特定的index。
3.BS-Seeker2针对WGBS和RRBS设计了不同的创建index的策略。

bs_seeker2-build.py \
    -f genome.fa \
    --aligner=bowtie2 \
    -r \                 # 加-r为构建rrbs索引，不加-r为构建WGBS索引
    -d index_dir         # -d参数为指定存放index的文件夹路径，不加-d index默认存放在/BSseeker2-master/bs_utils/reference_genomes/文件夹下

1.2 序列比对

1.测序数据一般较大，建议切割成小块进行比对，以提高比对效率。

先将质控后的双端测序的fastq.gz格式数据分成各个小块,

zcat test_clean.1.fastq.gz | split -l 16000000 - test.1_
zcat test_clean.2.fastq.gz | split -l 16000000 - test.2_
gzip test.*

2.比对时，BS-seeker2支持多种格式文件，输入文件可以为fastq、fasta、qseq或者纯序列格式及其gzip压缩文件。

bs_seeker2-align.py \
    -1 test.1_aa.gz \          # 输入切割后的左端测序文件
    -2 test.2_aa.gz \          # 输入切割后的右端测序文件
    -g genome.fa \             # 输入参考基因组序列
    --aligner=bowtie2 \        
    --bt2-p 12 \               # 设置程序运行线程数
    --bt2--end-to-end \        # 全局比对
    --bt2--very-sensitive \    
    --bt2--dovetail \
    --temp_dir=`pwd` \
    -d  index_dir \            # 输入存放index的文件夹路径。如果在建立索引时未设置-d，则此处不需加-d
    -m 1 \                   # 允许错配的碱基数，1代表只允许1个mismatch。
    --XSteve \   
    -o test_aa.bam

1.3 合并bam文件并根据染色体分离bam文件

samtools merge test_merge.bam `ls *.bam`
samtools sort test_merge.bam -o test_sort.bam

1.去除bam中的重复

java -jar /usr/local/bin/picard.jar MarkDuplicates \
        I=test_sort.bam \
        O=test_filter.bam \
        M=duplicate_stats.txt \
        REMOVE_DUPLICATES=true \
        AS=true

2.bam文件中标签XS:i:1代表BS-Seeker2利用CH位点的甲基化信息判断出亚硫酸盐转化失败的reads，将带有该信息的行舍弃

samtools view -h test_filter.bam | grep -v 'XS:i:1' | samtools view -bS - > test_filter_rmXS.bam  
samtools index test_filter_rmXS.bam

3.根据染色体分离出每条染色体的bam文件

for i in chr{1..5};do samtools view test_filter_rmXS.bam ${i} -h | samtools view -bS - > test_${i}.bam; done

1.4 Call methylation

1.可同时并行每条染色体call methylation

bs_seeker2-call_methylation.py \ 
    -i test_chr1A.bam \
    -o test_chr1A \
    -d index_dir \          # 输入存放index的文件夹路径。
    --rm-overlap \          # 确保双端测序reads重合的区域只被计算一次
    --sorted                # 代表读入的文件是已经经过排序的，可节省计算时间

--------------------2. CGmaptools----------------

安装CGmaptools

git clone git://github.com/guoweilong/cgmaptools.git

2.1 CGmaptools用法

1.将bam格式转换为CGmap格式

cgmaptools convert bam2cgmap -b AA_chr1.bam -g ara.fa --rmOverlap -o AA &

2.将ATCGmap转换为ATCGbz

cgmaptools convert atcgmap2atcgbz -c AA_chr1.ATCGmap.gz -b AA.ATCGbz

3.将CGmap转换为CGbz

cgmaptools convert cgmap2cgbz -c AA_chr1.CGmap.gz -b AA.CGbz

4.选取两个不同CGmap文件中的交集

cgmaptools intersect -1 AA_chr1.CGmap.gz -2 BB_chr1.CGmap.gz -o AA_BB.intersect.gz

-C 指定选取类型(CG\CHG\CHH)

5.合并两个ATCGmap或CGmap文件

cgmaptools merge2 atcgmap -1 AA_chr1.ATCGmap.gz -2 BB_chr1.ATCGmap.gz | gzip > AA_BB.merge.ATCGmap.gz
cgmaptools merge2 cgmap -1 AA_chr1.CGmap.gz -2 BB_chr1.CGmap.gz | gzip > AA_BB.merge.CGmap.gz

6.将多个CGmap的甲基化水平信息整合到一个文件

zcat AA_chr1.CGmap.gz BB_chr1.CGmap.gz | gawk '$8>=1' | cut -f1,3 | sort -u | cgmaptools sort -c 1 -p 2 > index
cgmaptools mergelist tomatrix -i index -f AA_chr1.CGmap.gz,BB_chr1.CGmap.gz -t AA,BB -c 1 -C 100 -o AA_BB.matrix.gz

-f 输入需要整合的CGmap文件，以","隔开
-t 输入对应文件的名称
-c 最小覆盖度
-C 最大覆盖度

7.根据染色体分离CGmap文件

cgmaptools split -i AA_chr1.CGmap.gz -p AA -s CGmap.gz

-p 输出文件名
-s 输出文件名后缀

8.输出目标区域的甲基化信息

awk 'BEGIN{OFS="\t";print 1,113000,124000,"+"}' > region.bed
zcat AA_chr1.CGmap.gz | cgmaptools select region -r region.bed > AA_region.CGmap

-R 反向选取，输出除输入文件中以外区域的甲基化信息
此步骤CGmap.gz需解压缩

9.输出目标位点的甲基化信息

gawk 'NR%100==50' index > site
cgmaptools select site -f AA_chr1.CGmap.gz -i site -o AA_site.CGmap.gz

-r 反向选取，输出除输入文件中以外位点的甲基化信息

10.Call SNP

cgmaptools snv -i AA_chr1.ATCGmap.gz -m bayes -v AA.vcf -o AA.snv --bayes-dynamicP
cgmaptools snv -i AA_chr1.ATCGmap.gz -m AA -o AA.snv

11.计算两个CGmap的差异甲基化位点(DMS)

cgmaptools dms -i AA_BB.intersect.gz -m 4 -M 100 -o AA_BB.DMS.gz -t fisher

-m 最小覆盖度
-M 最大覆盖度
-t 检验方法(fisher/chisq)

12.计算两个CGmap的差异甲基化区域(DMR)

cgmaptools dmr -i AA_BB.intersect.gz -c 1 -C 100 -o AA_BB.DMR.gz

-c 最小覆盖度
-C 最大覆盖度

13.Allele-Specific Methylation

gawk '{if(/^#/){print}else{print $0;}}' AA.vcf > AA1.vcf
cgmaptools asm -r ara.fa -b AA_chr1.bam -l AA1.vcf > AA1.asm

14.查看某一区域平均甲基化水平

zcat AA_chr1.CGmap.gz | cgmaptools mbed -b region.bed -c 1 -C 100

-c 最小覆盖度
-C 最大覆盖度

15.滑窗统计DNA甲基化水平(mbin)

cgmaptools mbin -i AA_chr1.CGmap.gz -B 50000 -c 1 -f png -t AA -p AA > mbin.AA.data

-B 滑窗大小
-c 最小覆盖度
-C 指定统计的DNA甲基化类型(CG/CHG/CHH)
-f 生成文件类型(png/pdf)
-t 图片的标题
-p 图片的名称

16.滑窗统计多个样本DNA甲基化水平(mmbin)

cgmaptools mmbin -l AA_chr1.CGmap.gz,BB_chr1.CGmap.gz -c 1 -B 2000 | \
gawk '{printf("%s:%s-%s", $1, $2, $3); for(i=4;i<=NF;i++){printf("\t%s", $i);} printf("\n");}' > mmbin

cgmaptools heatmap -i mmbin -c -o cluster.pdf -f pdf

-l 输入需要统计的CGmap文件，以","隔开
-c 最小覆盖度
-B 滑窗大小

17.将多个区域中的每个区域均分成多个片段，统计每个片段内的DNA甲基化水平(可用于统计基因上的DNA甲基化分布情况)

for CHR in 1 2 3 4 5; do (for P in 1 2 3 4 5; do echo | \
gawk -vC=$CHR -vP=$P -vOFS="\t" '{print C, P*10000, P*10000+1000, "+";}' ; done) ; done | \
cgmaptools bed2fragreg -n 5 -F 500,500,500,500,500 -T 500,500,500,500,500 > fragreg.bed

-n 将目标区域均分成n个片段
-F 在目标区域上游选取500,500,500,500,500片段
-T 在目标区域下游选取500,500,500,500,500片段

gunzip -c AA_chr1.CGmap.gz | cgmaptools mfg -r fragreg.bed -c 1 > AA.mfg

-c 最小覆盖度

-x 指定统计的DNA甲基化类型(CG/CHG/CHH)

18.统计CGmap文件中全部的DNA甲基化水平(mstat)

cgmaptools mstat -i AA_chr1.CGmap.gz -c 1 -f png -p AA -t AA > AA.mstat.data

-c 最小覆盖度
-f 生成文件类型(png/pdf)
-t 图片的标题
-p 图片的名称

19.计算目标区域DNA平均甲基化水平

cgmaptools mtr -i AA_chr1.CGmap.gz -r region.bed -o AA.mtr.gz

20.画某一区域的甲基化分布情况图(lollipop)

cgmaptools lollipop -i test.matrix.CG.gz -a test.refFlat -f pdf -c 1 -o test.pdf
最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
平台声明：文章内容（如有图片或视频亦包括在内）由作者上传并发布，文章内容仅代表作者本人观点，简书系信息发布平台，仅提供信息存储服务。

推荐阅读更多精彩内容