测序数据质控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
