# 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进行甲基化的检测 参考文档 ### 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%，所以总体服从**(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多进程脚本 **进行多重比较检验，使用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']] ##输出到新的文件 ``` ## 参考 FDR矫正