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/profile2.主要流程
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 )