Docker封装生物信息学甲基化流程
CpG CHG CHH含义
p代表磷酸二酯键，CpG指的是甲基化的C的下游是1个G碱基。H代表除了G碱基之外的其他碱基，即A, C, T中的任意一种，CHG代表甲基化的C下游的2个碱基是H和G, CHH表示甲基化的C下游的两个碱基都是H。
Bismark运行原理：
Bisulfite将序列正负链的C全部转换为T，所以也要将基因组序列进行转换。
基因组正负链转换很特别。
基因组的正链C->T，才能匹配原正链的reads
基因组的负链C->T相当于正链G->A，才能匹配原负链的reads
然后一条序列可以比对一个基因组的位置（即真实的基因组位置）
输出真实的基因组序列即可
bismark_methylation_extractor结果文件
默认情况下，软件会自动根据两个因素生成结果文件
1.甲基化的C的类型
就是前面提到的CpG, CHG, CHH 3种类型
2.比对情况
包括比对到四条链上OT, OB, CTOT, CTOB 4种情况
所以会生成 3 X 4 = 12 个文件，对于链特异性文库来说，会生成3 X 2 = 6 个文件，这6个文件内容是类似的，都是记录了甲基化的C的染色体位置。
在bismark中，将基因组的正链定义为top strand ， 简称OT, 负链定义为bottom strand, 简称OB; 亚硫酸氢盐处理后，正负链之间并不是完全的反向互补的，将OT链的反向互补链定义为CTOT, 将OB链的反向互补链定义为CTOB。对于链特异性文库而言，由于插入序列为单链，只需要比对OT和OB两条链即可，大大减少了运算量，所以目前illumina的标准BS-seq protocol构建的文库都是链特异性文库，新版的bismark默认的运行方式也是针对链特异性文库的。
1.流程文件 dna-meth.tar
导入镜像
docker load -i dna-meth.tar
2.未找到流程脚本--已下载
3.缺少基本命令
4.运行pipeline脚本
5.docker中时间与宿主机时间不一致（相差8小时）解决方法
docker cp /etc/localtime ef8f1c5e7533:/etc/localtime
1.环境配置
docker start 5fbb826d9572
docker exec -it 5fbb826d9572 /bin/bash
docker stop 5fbb826d9572
docker commit 5fbb826d9572 dna-meth:v4
docker build -t dna-meth:v4 .
中文格式设置
export LANG="C.UTF-8"
source /etc/profile
2.主要流程
1.参考基因组建立索引
2.数据QC
3.测序比对--bismark调用bowtie2
4.去除重复序列--bismarkDedup
5.统计比对率与去除重复之后的比对率
6.bismark_methylation_extractor 从去除重复的bam文件中提取出甲基化统计信息
1、质控
$fastqc -t $fastqc_threads --outdir Raw_FASTQC $R1.fastq.gz $R2.fastq.gz
java -jar $trimmomatic PE -phred33 -threads $trimmomatic_threads ../$R1.fastq.gz ../$R2.fastq.gz \
"$R1-paired.fastq" "$R1-unpaired.fastq" "$R2-paired.fastq" "$R2-unpaired.fastq" \
ILLUMINACLIP:$adapters:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:20 MINLEN:50 AVGQUAL:$base_quality_cutoff \
2> $sample.trimmomatic.Q$base_quality_cutoff.txt
$fastqc -t $fastqc_threads --outdir Trimmed_FASTQC Trimmed/$R1-paired.fastq Trimmed/$R2-paired.fastq
2、BisMark比对
$bismark_path/bismark -p 4 --non_directional --genome $genome_path --path_to_bowtie2 $bowtie_path --samtools_path $samtools_path -1 $1 -2 $2
# 比对到基因组
$bismark_path/deduplicate_bismark -p --samtools_path $samtools_path --bam $1
# 去重复
3.甲基化分析
$bismark_path/bismark_methylation_extractor --gzip --buffer_size $buffer_size --genome_folder $genome_path --parallel $parallel_processes $1
$bismark_path/bismark_methylation_extractor --gzip --buffer_size $buffer_size --genome_folder $genome_path --parallel $parallel_processes --comprehensive $1
$bismark_path/bismark2bedGraph --CX_context -o $2 $1
$bismark_path/bismark_methylation_extractor --gzip --buffer_size $buffer_size --cytosine_report --genome_folder $genome_path --parallel $parallel_processes $1
4.DMR区域分析
#!/usr/bin/env Rscript
args = commandArgs(trailingOnly=TRUE)
if (length(args)<3) {
stop("At least two inputs are needed input file, outFile and config JSON", call.=FALSE)
}
suppressPackageStartupMessages(library(dmrseq))
suppressPackageStartupMessages(library("BiocParallel"))
suppressPackageStartupMessages(library(rjson))
register(MulticoreParam(4))
df = read.table(args[1])
files <- c()
for(file in df$V1){
files <-c(files,file)
}
in_json = fromJSON(file = args[3], method = "C", unexpected.escape = "error", simplify = TRUE )
outFileName = args[2]
test = read.bismark(files = files,rmZeroCov = TRUE, strandCollapse = FALSE, verbose = TRUE)
#sampleNames<- c('control_R1','control_R2','treat_R1','treat_R2')
sampleNames <- c(strsplit(in_json$samplenames, " ")[[1]])
#cellType<- c('control','control','treat','treat')
cellType<- c(c(strsplit(in_json$CellType, " ")[[1]]))
pData(test)$CellType <- cellType
pData(test)$Replicate <- sampleNames
loci.idx <- which(DelayedMatrixStats::rowSums2(getCoverage(test, type="Cov")==0) == 0)
sample.idx <- which(pData(test)$CellType %in% unique(cellType))
test.filtered <- test[loci.idx, sample.idx]
testCovariate <- "CellType"
regions <- dmrseq(test.filtered,testCovariate=testCovariate)
#regions <- dmrseq(test.filtered, cutoff = 0.05, testCovariate=testCovariate)
out <- as(regions, "data.frame")
write.table(out, file = outFileName , sep="\t", col.names=TRUE, row.names=FALSE, quote=FALSE )
