溯源精微 BinBash WGBS分析流程

BinBash 计算平台为以个专业生物计算平台，提供超快速的WGS , Meta Genome, WGBS 等高计算量的样品进行快速计算服务。已经完成了大量的人，动植物的千例以上的数据加速工作。以下为一个WGBS的核心分析流程。

1. 质控

trim_galore -j $(nproc) --paired --basename $sam  -o /data/ $read1 $read2
#fastQC
fastqc -t $(nproc) -o /data/${sam}_fastqc \
                    $read1 \
                    $read2  \
                    /data/${sam}_R1_val_1.fq.gz   \                 
                    /data/${sam}_R2_val_2.fq.gz

   2. 数据切分

mkdir /data/${sam}_1 /data/${sam}_2
cd /data/${sam}_1
zcat /data/${sam}_R1_val_1.fq.gz  | split -d -a 10 -l 200000000 --filter='pigz  > ${FILE}.gz' 
cd /data/${sam}_2
zcat /data/${sam}_R2_val_2.fq.gz  | split -d -a 10 -l 200000000 --filter='pigz  > ${FILE}.gz' 
cd /data 
ls /data/${sam}_1 | while read id;do mkdir -p $id;done 
ls /data/${sam}_1 | while read id;do mv /data/${sam}_1/$id $id/${sam}_R1_val_1.fq.gz ;done 
ls /data/${sam}_2 | while read id;do mv /data/${sam}_2/$id $id/${sam}_R2_val_2.fq.gz ;done 
rm -rf /data/${sam}_2 /data/${sam}_1

3. 分块比对与合并

#mapping bismark 
block='X0000000000'
bismark --quiet --bam --rg_sample $sam --rg_id $sam  \
         -p $(nproc) --bowtie2 /database/hg38.fa/ -o ${block}/ \
          --fastq -1 $block/${sam}_R1_val_1.fq.gz -2 $block/${sam}_R2_val_2.fq.gz \
          --basename ${sam}.hg38.fa --samtools_path /usr/bin --temp_dir $block/
          
#dedup S005413.deduplication_report.txt; 8C64G 
samtools cat -@ $(nproc) -o $sam.bam --no-PG  X0000*/${sam}.hg38.fa_pe.bam

4. 比对结果统计与去重复

#bamqc
samtools sort -@ $(nproc)  -o /data/$sam.sort.bam -O BAM /data/$sam.bam
qualimap  bamqc -bam /data/$sam.sort.bam -outdir /data/bamqc_result -outformat PDF:HTML
deduplicate_bismark --bam --output_dir /data/ --samtools_path /usr/bin -p $sam.bam 

5. 甲基化提取

#人WGBS 需要46GB左右内存
/usr/bin/time bismark_methylation_extractor --samtools_path /usr/bin --parallel ${thread} \
       --comprehensive -p --no_overlap --report --gzip --cytosine_report --CX_context \
       --buffer_size ${mem}M --genome_folder /database/hg38.fa /data/$sam.deduplicated.bam\
       -o /data/

6. 差异甲基化

  差异甲基化部分，使用methylKit 软件包进行

#差异甲基化  methylKit R 
#处理/data/${sam}.deduplicated.CX_report.txt.gz  并转换成methylKit 能识别的格式 
分开CpG, CHH, CHG
echo -e "chrBase\tchr\tbase\tstrand\tcoverage\tfreqC\tfreqT" > $sam.methylKit.CG.txt
echo -e "chrBase\tchr\tbase\tstrand\tcoverage\tfreqC\tfreqT" > $sam.methylKit.CHH.txt
echo -e "chrBase\tchr\tbase\tstrand\tcoverage\tfreqC\tfreqT" > $sam.methylKit.CHG.txt

zcat ${sam}.deduplicated.CX_report.txt.gz | \
    awk -v OFS="\t" -v sam=$sam  '{ if($3=="+"){s="F"}else{s="R"};total=$4+$5;
                        print $1"."$2,$1,$2,s,total,100*$4/total,100*$5/total >> sam".methylKit."$6".txt"}'
                
 
                        
#数据很大，需要大内存                                                                                      
library(methylKit)
library(graphics)

file.list <- list(system.file("extdata", "test1.myCpG.txt",package= "methylKit"), 
                  system.file("extdata",
"test2.myCpG.txt",package= "methylKit"), 
                  system.file("extdata","control1.myCpG.txt",package= "methylKit"), 
                  system.file("extdata","control2.myCpG.txt",package= "methylKit"))   
myobj <- read(file.list, sample.id= list("test1","test2", "ctrl1", "ctrl2"),
               assembly = "hg38", treatment=
c(1,1, 0, 0), context = "CpG")  
getMethylationStats(myobj[[2]], plot= T, both.strands = F)                                   
getCoverageStats(myobj[[2]],plot = T, both.strands = F)

#过滤代码过滤甲基化数据并丢弃覆盖率低于10X的碱基，并且丢弃每个样本中覆盖率超过99.9％的碱基。
filtered.myobj <- filterByCoverage(myobj,lo.count = 10,lo.perc = NULL,hi.count = NULL, hi.perc = 99.9)  
#整合样品
meth <- unite(myobj, destrand= FALSE)
head(meth)
#样品相关性
getCorrelation(meth, plot= True)
#PCA 
PCASamples(meth,screeplot=True)
PCASamples(meth)

#聚类
clusterSamples(meth,dist='correlation',method='',plot=True)

#差异甲基化位点  
myDiff <- calculateDiffMeth(meth)

#在q值计算之后，我们可以根据q值和甲基化差异截止百分比选择差异甲基化区域/碱基。
#以下位选择q值<0.01且甲基化差异百分比大于25％的碱基。
#如果指定type =“hyper”或type =“hypo”选项，则会出现超甲基化或低甲基化区域/碱基。

#get hyper methylated bases
myDiff25p.hyper <- get.methylDiff(myDiff,difference = 25,qvalue = 0.01,type = "hyper")
# get hypo methylated bases
myDiff25p.hypo <- get.methylDiff(myDiff,difference = 25,qvalue = 0.01,type = "hypo")

#get all differentially methylated bases
myDiff25p <- get.methylDiff(myDiff,difference = 25,qvalue = 0.01)

#我们可以基于基因注释来注释我们的差异甲基化区域/碱基。
#在这个例子中，从bed文件中读取基因注释，并用这些信息注释差异甲基化区域。
#在启动子/内含子/外显子/基因间区域中差异甲基化区域的百分比。
#使用Bioconductor.org提供的GenomicFeatures软件包可以获取类似的基因注释    
gene.obj<- read.transcript.features(system.file("extdata","refseq.hg38.bed.txt", 
                                    package= "methylKit"))

diffAnn<-annotate.WithGenicParts(myDiff25p,gene.obj) 
plotTargetAnnotation(diffAnn,precedence=TRUE,main="differential methylation annotation")   
 


#CpG岛注释    
cpg.obj <- read.feature.flank(system.file("extdata","cpgi.hg38.bed.txt", 
                              package= "methylKit"),
feature.flank.name = c("CpGi","shores"))

diffCpGann <- annotate.WithFeature.Flank(myDiff25p,
                     cpg.obj$CpGi,cpg.obj$shores, 
                     feature.name = "CpGi",
                     flank.name= "shores")
#甲基化区域MDR 按照以上操作即可
tiles <- titleMethylCounts(myobj,win.size = 1000)
head(titles[[1]])
MDRdif <-calculateDiffMeth(tiles)

#差异甲基化区域对应基因的Go , Pathway富集分析

                                                                                                                                  

                                                             

                                 公众号                                                                     微信