测序数据质控fastp

fastp -i ../cleandata/$seq\_R1.fastq.gz -o $seq.R1.clean.fq.gz \  ##read1
-I ../cleandata/$seq\_R2.fastq.gz -O $seq.R2.clean.fq.gz     \  ##read2
-5 --cut_front_window_size 4   --cut_front_mean_quality 20 \     ## from 5'-, window_slide 4bp,if quality <20,cut, NOTE:affect deduplication
-3 --cut_tail_window_size  4   --cut_tail_mean_quality 20  \     ## from 3'-, window_slide 4bp,if quality <20,cut, NOTE:affect deduplication
--cut_right --cut_right_window_size 4 --cut_right_mean_quality 20 \   ## cut right if <threshold , NOTE:affect deduplication
--detect_adapter_for_pe     \   ##remove adaptor
-q 15 \  ## base quality >15
-u 40 \  ## keep reads which had less than 40% low quality bases
-e 20 \  ## read average quality >20
-n 5 \   ## Ns number <=5
-l 30 \  ## read length >=30
-p \     ## enable overrepresented sequence analysis. small p
-P 20 \  ## One in 1/20 reads will be computed for overrepresentation analysis, big P
-w 4 \   ## worker thread number
#-c   \  ## overlap analysis for PE data ,bismark will do this

genome preparation

bismark_genome_preparation --bowtie2 --verbose --genomic_composition ./hg19.fa

alignment

biamrk --bowtie2 \
           -q \                 ## fastq
           -N  1 \                ## mismatch <=1
           -p 8 \               ## 8 threads
          --multicore 4 \       ## use 4*4 cores 
          --non_directional \   ## alignment to 4 strands
          -o ./           \     ## output directary
        --gzip          \       ## gzip output
        --nucleotide_coverage \  ## coverage
        -1 $i.R1.fq.gz   \
        -2 $i.R2.fq.gz  

bismark_deduplication

deduplicate_bismark -p $i.bam  

picard_deduplication

samtools sort $i.bam -o $i.sorted.bam   ## sort by position

java -jar ~/anaconda2/share/picard-2.21.6-0/picard.jar \  ##picard.jar absolute path
MarkDuplicates REMOVE_DUPLICATES   \  ## function 
-I $i.sorted.bam -O $i.dedup.bam -M picard1.txt  \   ## input and output

samtools sort -n  $i.dedup.bam -o $i.dedup.sorted.bam ## resort by reads name

extraction_DNA_methylatiuon

bismark_methylaiton_extractor -p \  ##paired-end
                              --no_overlap  \  ##avoids scoring overlapping methylation calls twice
                              --comprehensive \  ##merge all four possible strand-specific methylation info into context-dependent output files
                              --report \  ##methylation summary
                              -o ./  \  ## output file
                              --gzip  \  ## gzip output
                              --multicore 4 \  ## use multicores 4*4
                              --buffer_size 20G \  ## main memory sort buffer when sorting the methylation information
                              --CX  \  ##information on every single cytosine in the genome irrespective of its context
                              --cytosine_report \  ## genome-wide methylation report for all cytosines in the genome
                              --genome_folder ../genomefolder/  \  ## reference genome path
                              --zero_based \  ## 0-based or 1-based
                              $i.dedup.sorted.bam

extract CG positions

 perl 0_cx_report2CG.pl $i.CX_report.txt.gz $i.CG.txt

merged CG postition (merge + and - strands data)

perl  1_CG2merged_CG.pl $i.CG.txt $i.merged_CG.txt 

extractor candiate DMPs data

perl candiate.pl $i.CG-merged.txt DMPs_candiate.txt >$i.dmps_data.txt