# 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>

### 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>

**进行多重比较检验，使用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>

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


---

# Agent Instructions: Querying This Documentation

If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.

Perform an HTTP GET request on the current page URL with the `ask` query parameter:

```
GET https://zpliu.gitbook.io/booknote/ke-bian-jian-qie/ruan-jian-shi-yong/04-jia-ji-hua-fen-xi-liu-cheng.md?ask=<question>
```

The question should be specific, self-contained, and written in natural language.
The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.