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 )

本次流程主要软件已安装,直接测试,后续增加功能--GO与KEGG富集分析。

测试通过后,整理报告模板,docker打包并部署。

(0)

相关推荐