# DNA甲基化分析流程

## 1.对原始数据进行过滤

```bash
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](https://github.com/FelixKrueger/Bismark/tree/master/Docs)

### 2.1 准备基因组索引文件

* 程序会对基因组序列进行C->T and G->A 的转化
* 之后使用hisat2-build进行索引的构建

```bash
#指定基因组文件所在文件夹
./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 输出文件前缀

```bash
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
```

产生可变剪切位点文件

```bash
## 使用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脚本

```bash
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.

```bash
# 当对双端测序数据进行去重时，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 提取程序，在默认参数下使用

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

* -p 指定双端测序数据
* \--parallel 指定核数
* \--bedGraph 输出bed文件
* \--gzip 输出文件使用gzip压缩

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

```bash
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%，所以总体服&#x4ECE;**(n,0.006)的二项分布**
> * 进行假设检验，使用p-value=0.05作为阈值，当某个甲基化位点达到阈值，才被认为是真实的
> * 同时至少有3条read支持这个甲基化位点

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

```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数目

```bash
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位点并不是由于试剂的原因，而是真实存在的甲基化位点

```r
# 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)
```

得到的结果文件

```bash
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](https://zpliu.gitbook.io/booknote/script/binomtest)

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

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

```r
#原始数据格式，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

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

```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文件太大了分开跑

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

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

### 合并多个文件的内容

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

`merge.sh`

```bash
#!/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
```

```bash
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数

    ```bash
    检测到的坐标 支持的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中进行多重测验

```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](https://github.com/FelixKrueger/Bismark/tree/master/Docs)

FDR矫正 [https://blog.csdn.net/zhu\_si\_tao/article/details/71077703](https://blog.csdn.net/zhu_si_tao/article/details/71077703)
