1.对原始数据进行过滤
Copy 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 的转化
Copy #指定基因组文件所在文件夹
./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 进行比对
non_direction 跟建库方式有关 默认是direction
-N 设置允许错配数 只在bowtie2模式下有用后面要考虑 --score-min L,0,-0.2
--known-splicesite-infile 提供剪切位点文件，只在hisat2比对模式下可用
Copy 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 产生可变剪切位点文件
Copy ## 使用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脚本
Copy 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.
Copy # 当对双端测序数据进行去重时，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 2.4 提取对应的甲基化位点
将比对的bam文件经过过滤之后，使用Bismark methylation extractor 提取程序，在默认参数下使用
Copy bismark_methylation_extractor -p --parallel 2 --bedGraph --gzip -o ./ ./${fileName}.deduplicated.bam 因为甲基化的数据量非常大，程序根据甲基化的类型结合比对时，链的关系可以得到12个文件
Copy 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作为阈值，当某个甲基化位点达到阈值，才被认为是真实的
自己写了一个python脚本进行统计，每个甲基化位点的甲基化概率
Copy 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)) 输出结果
Copy 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位点并不是由于试剂的原因，而是真实存在的甲基化位点
Copy # 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) 得到的结果文件
Copy 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 进行多重比较检验，使用FDR方法
FDR：错误拒接率，意思就是错误的拒绝了H0假设，是假阳性，通过FDR校验保证数据总体假阳性不超过阈值,并且进行多重校验的数据集是需要满足二项分布的数据集的p-value和read数目超过3的。
Copy #原始数据格式，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
在提交任务的时候，文件路径不能有 . 不然有错误
Copy 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文件太大了分开跑
Copy split -l 行数 原始文件 分割文件前缀 对每个文件分别提取支持的read与不支持的read
合并多个文件的内容
直接使用awk将多个文件的内容进行合并
merge.sh
Copy #!/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
Copy 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数
Copy 检测到的坐标 支持的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中进行多重测验
Copy 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']]
##输出到新的文件 参考