DNA甲基化分析流程

1.对原始数据进行过滤

java -jar trimmomatic-0.32.jar PE -threads 18 -phred33 R1.fastq R2.fastq -baseout 输出文件目录/文件前缀 ILLUMINACLIP:TruSeq3-PE.fa:2:30:10 LEADING:10 TRAILING:10 SLIDINGWINDOW:4:20 MINLEN:50 >log.${ID}.trim 2>err.${ID}.trim

2.使用Bismark进行甲基化的检测

参考文档 https://github.com/FelixKrueger/Bismark/tree/master/Docs

2.1 准备基因组索引文件

• 程序会对基因组序列进行C->T and G->A 的转化

• 之后使用hisat2-build进行索引的构建

#指定基因组文件所在文件夹
./bismark_genome_preparation --hisat2 travis_files/
#之后会生成两个转化后的基因组文件，及对应的索引
├── Bisulfite_Genome
│   ├── CT_conversion
│   │   ├── BS_CT.1.ht2
│   │   ├── BS_CT.2.ht2
│   │   ├── BS_CT.3.ht2
│   │   ├── BS_CT.4.ht2
│   │   ├── BS_CT.5.ht2
│   │   ├── BS_CT.6.ht2
│   │   ├── BS_CT.7.ht2
│   │   ├── BS_CT.8.ht2
│   │   └── genome_mfa.CT_conversion.fa
│   └── GA_conversion
│       ├── BS_GA.1.ht2
│       ├── BS_GA.2.ht2
│       ├── BS_GA.3.ht2
│       ├── BS_GA.4.ht2
│       ├── BS_GA.5.ht2
│       ├── BS_GA.6.ht2
│       ├── BS_GA.7.ht2
│       ├── BS_GA.8.ht2
│       └── genome_mfa.GA_conversion.fa

2.2 进行比对

• geneme_folder 指定基因组文件夹

• --hisat2 使用hisat2进行比对

• non_direction 跟建库方式有关 默认是direction

• -L 设置种子序列长度

• -N 设置允许错配数 只在bowtie2模式下有用后面要考虑 --score-min L,0,-0.2

• --parallel 指定核心数

• -p 指定线程数

• --known-splicesite-infile 提供剪切位点文件，只在hisat2比对模式下可用

• -o 输出文件前缀

bismark  --genome_folder travis_files  --hisat2    --non_directional -1 travis_files/test_R1.fastq.gz  -2 travis_files/test_R2.fastq.gz  -L  30 -N 1 -o  test

产生可变剪切位点文件

## 使用cufflinks内的脚本将gff3转化为gtf
gff2gtf_cufflinks -T  Graimondii_221_v2.1.gene.gff3 -o  Graimondii_221_v2.1.gene.gtf
## 使用hisat2中脚本提取可变剪切
extract_splice_sites.py Graimondii_221_v2.1.gene.gtf   >Graimondii_221_v2.1.gene_splice.txt

进行hisat2模式的比对，lsf脚本

softwarwPath=~/software/Bismark-0.22.1/ 
genomePath=基因组文件夹
splicePath=可变剪切文件
R1Path=./../D4.R1.trimed.fastq.gz 
R2Path=./../D4.R2.trimed.fastq.gz

${softwarwPath}/bismark_genome_preparation  --hisat2 ${genomePath} >log.bismark_genome_preparation 2>err.bismark_genome_preparation

${softwarwPath}/bismark --genome_folder ${genomePath} --hisat2    --non_directional --known-splicesite-infile ${splicePath} -1 ${R1Path} -2 ${R2Path} -L 30 -N 1 -o ./ -p 5

2.3 过滤比对得到的bam文件

当使用全基因组进行亚硫酸氢盐处理得到的文库时，建议进行这一步操作,去除 Bismark 比对的重复reads, 因为PCR跑的循环数多了的话，会导致另外一条链上比对到的read数增加。

• excessive PCR amplification

• Sequences which align to the same genomic position but on different strands are scored individually.

# 当对双端测序数据进行去重时，bam文件需要安装read名字排好序
# bismark的输出结果是排好序的
deduplicate_bismark  -p -o 输出文件前缀 --output_dir 输出目录 --bam  比对后的BAM输出文件

### 批量脚本
softwarwPath=~/software/Bismark-0.22.1/ 
fileName=D1.R1

${softwarwPath}/deduplicate_bismark  -p -o  ${fileName} --output_dir ./ --bam ${fileName}.trimed_bismark_hisat2_pe.bam

• -p 指定为双端测序

• -o 输出文件前缀

• --output_dir 输出文件目录

• --bam 最终去重后为BAM格式

2.4 提取对应的甲基化位点

将比对的bam文件经过过滤之后，使用Bismark methylation extractor 提取程序，在默认参数下使用

bismark_methylation_extractor -p  --parallel 2  --bedGraph --gzip  -o ./   ./${fileName}.deduplicated.bam

• -p 指定双端测序数据

• --parallel 指定核数

• --bedGraph 输出bed文件

• --gzip 输出文件使用gzip压缩

因为甲基化的数据量非常大，程序根据甲基化的类型结合比对时，链的关系可以得到12个文件

OT      original top strand
CTOT    complementary to original top strand

OB      original bottom strand
CTOB    complementary to original bottom strand
# 3种甲基化 CpG CHG CHH *4 =12个文件

当然对链特异性甲基化数据不感兴趣时，也可以使用--comprehensive参数，只生成3中甲基化的报告文件，文件的数据是由四种链进行合成的

2.5基于统计上的显著性提取对应的位点

在进行建库的时候，使用了λ噬菌体的基因组来获得Bisufite处理的转化效率，大概是99.4%，接下来使用二项分布对每个位点甲基化的显著性进行检验

• 每个甲基化的C被转化的概率为99.4%；不被转化的概率为0.6%，所以总体服从(n,0.006)的二项分布

• 进行假设检验，使用p-value=0.05作为阈值，当某个甲基化位点达到阈值，才被认为是真实的

• 同时至少有3条read支持这个甲基化位点

自己写了一个python脚本进行统计，每个甲基化位点的甲基化概率

import sys
import gzip

'''
python $0 bismark输出文件 最终统计文件
最终文件名字以gz结尾
'''
with gzip.open(sys.argv[1], 'r') as inputfile:
    Alllines = inputfile.readlines()
Allindex = {}
for i in range(1, len(Alllines)):
    Alllines[i] = Alllines[i].decode()  # 对byte字符进行解码
    index = Alllines[i].strip("\n").split(
        "\t")[2]+"~"+Alllines[i].strip("\n").split("\t")[3]

    Allindex[index] = [0, 0]


for i in range(1, len(Alllines)):
    tag = Alllines[i].strip("\n").split("\t")[1]
    index = Alllines[i].strip("\n").split(
        "\t")[2]+"~"+Alllines[i].strip("\n").split("\t")[3]
    if(tag == '+'):
        Allindex[index][0] += 1
    else:
        Allindex[index][1] += 1

with gzip.open(sys.argv[2], 'w') as outputfile:
    for key in Allindex:
        outstr = key.split("~")[0]+"\t"+key.split("~")[
            1]+"\t"+str(Allindex[key][0])+"\t"+str(Allindex[key][1])+"\n"
        outputfile.write(str.encode(outstr))

输出结果

• 前两列，甲基化的位置

• 支持该位点的read数目

• 不支持该位点的read数目

Chr05   68302093        3       0
Chr05   68302156        2       0
Chr05   68302192        3       0
Chr05   68302203        3       0
Chr05   68302237        1       1
Chr05   68302282        3       0

进行二项分布检验

运行脚本Rscript bunomTest.R test.txt

• 对每个位点，假定的甲基化位点reads超过3条的进行二项分布检验

  没有超过3个read的位点会被抛弃掉，在之后的分析中会带来偏差

• 在每个没有被甲基化位点仍旧保持C的概率为0.006，H0假设，该位点是由于试剂的原因导致原本没有甲基化的C仍旧保持为C，因此当p-value小于阀值时，拒绝原假设，认为这个C位点并不是由于试剂的原因，而是真实存在的甲基化位点

# bunomTest.R脚本如下
args <- commandArgs(TRUE)
indata <- args[1]
mydata <- read.table(indata,header=F)
myindex <- NULL ;
p.values <- NULL ;
for (i in 1:dim(mydata)[1]){
   if(mydata[i,3]>=3){
    count <- mydata[i,3]+mydata[i,4] 
      p.value <- dbinom(mydata[i,3],count,0.006)
     if (p.value <= 1e-5){
      myindex <- c(myindex, i) ;
      p.values <- c(p.values, p.value) ;
   }
}
}
outdata <- cbind(mydata[myindex,],as.numeric(p.values))
outfileName <- paste(strsplit(indata,'\\.')[[1]][1],"binomTest_result.txt",sep="_")
write.table(outdata,file=outfileName,quote=FALSE,sep="\t",col.name=FALSE,row.name=FALSE)

得到的结果文件

Chr12    13159675    4    0    1.296e-09
Chr12    13159661    5    0    7.77599999999999e-12
Chr12    13159642    9    0    1.0077696e-20
Chr12    13159618    12    0    2.17678233599999e-27

前面的R脚本太慢了，数据量大的话要跑很久很久，使用我另外写的一个python多进程脚本 https://zpliu.gitbook.io/booknote/script/binomtest

进行多重比较检验，使用FDR方法

FDR：错误拒接率，意思就是错误的拒绝了H0假设，是假阳性，通过FDR校验保证数据总体假阳性不超过阈值,并且进行多重校验的数据集是需要满足二项分布的数据集的p-value和read数目超过3的。

#原始数据格式，V5为p-value值
           V1       V2 V3 V4          V5           
1       Chr13 13682239  4  0 1.29600e-09 
2       Chr13 13682246  3  0 2.16000e-07 
3       Chr13 13682254  3  1 8.58816e-07 l
4 tig00019827    35589  3  1 8.58816e-07 
5 tig00019827    35602  3  1 8.58816e-07 
6 tig00019827    35609  3  0 2.16000e-07 
testdata=read.table("./test.txt")
testdata$V6= p.adjust(testdata$V5,"fdr",dim(testdata)[1])
write.table(testdata[testdata$V6<=0.05,],"111",sep="\t",col.names=F,row.names=F,quote=F)

fdr_correction.R

在提交任务的时候，文件路径不能有.不然有错误

args <- commandArgs(TRUE)
indata <- args[1]
mydata <- read.table(indata,header=F)
mydata$V6= p.adjust(mydata$V5,"fdr",dim(mydata)[1])
outfileName <- paste(strsplit(indata,'\\.')[[1]][1],"fdr_result.txt",sep="_")
write.table(mydata[mydata$V6<=0.05,],file=outfileName,quote=FALSE,sep="\t",col.name=FALSE,row.name=FALSE)

CHH文件太大了分开跑

split -l 行数 原始文件  分割文件前缀

对每个文件分别提取支持的read与不支持的read

合并多个文件的内容

直接使用awk将多个文件的内容进行合并

merge.sh

#!/bin/bash 

cat ./split*|awk -F "\t" '{a[$1"~"$2][1]+=$3;a[$1"~"$2][2]+=$4}END{for(i in a){split(i,b,"~");print b[1]"\t"b[2]"\t"a[i][1]"\t"a[i][2]}}' >merge_end

bsub -q high -M 200G -e merge.err -R span[hosts=1] -n 1 -J merge "bash merge.sh"

统计没有甲基化的位点与发生了甲基化的位点

• 由于测序深度的原因，在计算DNA甲基化的时候，有点胞嘧啶可能没有被测到可能会导致结果偏低，因此仅仅考虑被测到的胞嘧啶。

• 区域内所有甲基化C位点总的reads数/区域内总的覆盖度

  FDR值达到显著性的支持read数/FDR达到显著性支持read数+FDR达到显著性不支持数+FDR没有达到显著性值的read数

  检测到的坐标 支持的read数，不支持的read数，FRD值
  Ghir_D02    4969712    4969712    0    12    0.994
  Ghir_A13    66114242    66114242    12    0    1.11882196473219e-26
  Ghir_D07    17177217    17177217    16    0    3.31847691703274e-35

在python中进行多重测验

import panda as pd 
import numpy as np
import statsmodels.stats.multitest as multi
##读取原始数据
data=pd.read_csv("./CpG_read_count_bionormal.txt",sep="\t",header=None,low_memory=False)
##删除一些SRR开头的行
dropRow=data[data[0].str.startswith('SRR')]
data.drop(dropRow.index,inplace=True)
##对p-value进行多重校验，获得数组
a=multipletests(np.array(data[4]))
a[1]返回FDR的值
##在原有数据框上添加FDR列
data.loc[:,'fdr']=multi.multipletests(np.array(data[4]))[1]
##查看指定位置的值
data[data[1]==67304209]
##增加同一个坐标作为bed文件
data.loc[:,5]=data[1]
##调整输出顺序
data=data[[0,1,5,2,3,4,'fdr']]
##输出到新的文件

参考

https://github.com/FelixKrueger/Bismark/tree/master/Docs

FDR矫正 https://blog.csdn.net/zhu_si_tao/article/details/71077703