Skip to content

03.BSSEQ_extract_methylation.md

Latest commit

 

History

History
94 lines (73 loc) · 6.9 KB

03.BSSEQ_extract_methylation.md

File metadata and controls

94 lines (73 loc) · 6.9 KB
Raw

BS-Seq: Extract methylation

🏠 home

Note

Extract methylation from the trimmed down reads

Prepare diploid genomes

Bismark-0.24.2/bismark_genome_preparation --bowtie2 --parallel 12 ${ref_dir} > logs/bismark_genome_preparation_${ref_name}.log 2> logs/bismark_genome_preparation_${ref_name}.err

Run bismark on diploid genomes

--ambiguous: to get the reads multi-mapped --ambig_bam: to get a random alignment

for i in $(ls down/${ref_name}*fq.gz | grep R1); do name=$(basename $i .fq.gz | sed 's:_R1.*$::g'); Bismark-0.24.2/bismark --parallel 8 --output_dir alignment/bismark/${name}.on.${ref_name}.bismark --bowtie2 -p 4 -N 0 --ambiguous --ambig_bam ${ref_dir} -1 down/${name}_R1.fil.fq.gz -2 down/${name}_R2.fil.fq.gz > logs/${name}.on.${ref_name}.bismark.log 2> logs/${name}.on.${ref_name}.bismark.err

Output:

  • a bam file
  • a report with percentages of methylated C per context.
  • ambiguous reads + ambiguous bam

Run ambiguous reads per haplotype

for i in $(ls alignment/bismark/${ref_name}*on.${ref_name}.bismark/*ambiguous_reads_1.fq.gz); do name=$(basename $i .fq.gz | sed 's:_R1.*$::g'); Bismark-0.24.2/bismark --parallel 8 --output_dir alignment/bismark/${name}_ambig.on.${hap_name}.bismark --bowtie2 -p 4 -N 0 genomes/${hap_name}/ -1 alignment/bismark/${name}.on.${ref_name}.bismark/${name}_R1.fil.fq.gz_ambiguous_reads_1.fq.gz -2 alignment/bismark/${name}.on.${ref_name}.bismark/${name}_R2.fil.fq.gz_ambiguous_reads_2.fq.gz > logs/${name}_ambig.on.${hap_name}.bismark.log 2> logs/${name}_ambig.on.${hap_name}.bismark.err; done

Concatenate the bam

extract haplotype bams from diploid bam

mkdir alignment/bismark/${sample}.on.${hap1}.bismark_hybrid; samtools view -b -h -L genomes/${hap1}/${hap1bed} -@ 24 alignment/bismark/${sample}.on.${ref}.bismark/${sample}_R1.fil_bismark_bt2_pe.bam > alignment/bismark/${sample}.on.${hap1}.bismark_hybrid/${sample}.on.${hap1}_bismark_bt2_pe.bam 

mkdir alignment/bismark/${sample}.on.${hap2}.bismark_hybrid; samtools view -b -h -L genomes/${hap2}/${hap2bed} -@ 24 alignment/bismark/${sample}.on.${ref}.bismark/${sample}_R1.fil_bismark_bt2_pe.bam > alignment/bismark/${sample}.on.${hap2}.bismark_hybrid/${sample}.on.${hap2}_bismark_bt2_pe.bam

create fofn files of file names for concatenating

echo -e "alignment/bismark/${sample}.on.${hap1}.bismark_hybrid/${sample}.on.${hap1}_bismark_bt2_pe.bam\nalignment/bismark/${sample}_ambig.on.${hap1}.bismark/${sample}_R1.fil.fq.gz_ambiguous_reads_1_bismark_bt2_pe.bam" > alignment/bismark/${sample}.on.${hap1}.bismark_hybrid/hybrid_mapping_files.${sample}.on.${hap1}.txt

echo -e "alignment/bismark/${sample}.on.${hap2}.bismark_hybrid/${sample}.on.${hap2}_bismark_bt2_pe.bam\nalignment/bismark/${sample}_ambig.on.${hap2}.bismark/${sample}_R1.fil.fq.gz_ambiguous_reads_1_bismark_bt2_pe.bam" > alignment/bismark/${sample}.on.${hap2}.bismark_hybrid/hybrid_mapping_files.${sample}.on.${hap2}.txt

concatenate

bam1=$(cat alignment/bismark/${sample}.on.${hap1}.bismark_hybrid/hybrid_mapping_files.${sample}.on.${hap1}.txt | tr '\n' ' ')
bam2=$(cat alignment/bismark/${sample}.on.${hap2}.bismark_hybrid/hybrid_mapping_files.${sample}.on.${hap2}.txt | tr '\n' ' ')

samtools merge --threads 24 -n -f -o alignment/bismark/${sample}.on.${hap1}.bismark_hybrid/${sample}.on.${hap1}_bismark_hybrid_bt2_pe.bam ${bam1}; samtools sort --threads 24 -n -o alignment/bismark/${sample}.on.${hap1}.bismark_hybrid/${sample}.on.${hap1}_bismark_hybrid_bt2_pe.sorted.bam alignment/bismark/${sample}.on.${hap1}.bismark_hybrid/${sample}.on.${hap1}_bismark_hybrid_bt2_pe.bam

samtools merge --threads 24 -n -f -o alignment/bismark/${sample}.on.${hap2}.bismark_hybrid/${sample}.on.${hap2}_bismark_hybrid_bt2_pe.bam ${bam2}; samtools sort --threads 24 -n -o alignment/bismark/${sample}.on.${hap2}.bismark_hybrid/${sample}.on.${hap2}_bismark_hybrid_bt2_pe.sorted.bam alignment/bismark/${sample}.on.${hap2}.bismark_hybrid/${sample}.on.${hap2}_bismark_hybrid_bt2_pe.bam

Deduplicate

run the deduplicate

mkdir -p alignment/bismark/${sample}.on.${hap1}.bismark_hybrid/dedup; Bismark-0.24.2/deduplicate_bismark --paired --bam --output_dir alignment/bismark/${sample}.on.${hap1}.bismark_hybrid/dedup alignment/bismark/${sample}.on.${hap1}.bismark_hybrid/${sample}.on.${hap1}_bismark_hybrid_bt2_pe.sorted.bam > logs/${sample}.on.${hap1}_bismark_hybrid.bismark_dedup.log 2> logs/${sample}.on.${hap1}_bismark_hybrid.bismark_dedup.err

mkdir -p alignment/bismark/${sample}.on.${hap2}.bismark_hybrid/dedup; Bismark-0.24.2/deduplicate_bismark --paired --bam --output_dir alignment/bismark/${sample}.on.${hap2}.bismark_hybrid/dedup alignment/bismark/${sample}.on.${hap2}.bismark_hybrid/${sample}.on.${hap2}_bismark_hybrid_bt2_pe.sorted.bam > logs/${sample}.on.${hap2}_bismark_hybrid.bismark_dedup.log 2> logs/${sample}.on.${hap2}_bismark_hybrid.bismark_dedup.err

Output:

  • a deduplicated bam
  • a report with the number of duplicated alignments removed

Extract methylation

Bismark-0.24.2/bismark_methylation_extractor --paired-end --gzip --multicore 24 --buffer_size 12G --CX_context --cytosine_report --comprehensive --bedGraph --ignore 15 --ignore_r2 15 --ignore_3prime 2 --ignore_3prime_r2 2 --output alignment/bismark/${sample}.on.${hap1}.bismark_hybrid/dedup/bismark_methylation_extractor --genome_folder genomes/${hap1}/ alignment/bismark/${sample}.on.${hap1}.bismark_hybrid/dedup/${sample}.on.${hap1}_bismark_hybrid_bt2_pe.sorted.deduplicated.bam > logs/${sample}.on.${hap1}_bismark_hybrid.bismark_extract.log 2> logs/${sample}.on.${hap1}_bismark_hybrid.bismark_extract.err

Bismark-0.24.2/bismark_methylation_extractor --paired-end --gzip --multicore 24 --buffer_size 12G --CX_context --cytosine_report --comprehensive --bedGraph --ignore 15 --ignore_r2 15 --ignore_3prime 2 --ignore_3prime_r2 2 --output alignment/bismark/${sample}.on.${hap2}.bismark_hybrid/dedup/bismark_methylation_extractor --genome_folder genomes/${hap2}/ alignment/bismark/${sample}.on.${hap2}.bismark_hybrid/dedup/${sample}.on.${hap2}_bismark_hybrid_bt2_pe.sorted.deduplicated.bam > logs/${sample}.on.${hap2}_bismark_hybrid.bismark_extract.log 2> logs/${sample}.on.${hap2}_bismark_hybrid.bismark_extract.err

Output: https://felixkrueger.github.io/Bismark/bismark/methylation_extraction/

  • a bedgraph (bedGraph.gz)
<chromosome> <start position> <end position> <methylation percentage>
  • a coverage file (bismark.cov.gz)
<chromosome> <start position> <end position> <methylation percentage> <count methylated> <count unmethylated>
  • a genome-wide cytosine report (CX_report.txt.gz)
<chromosome> <position> <strand> <count methylated> <count unmethylated> <C-context> <trinucleotide context>
  • a default output per context (XXX_context_XXX.txt.gz)
  • M-bias R1 and R2 png + M-bias.txt
  • a summary report (splitting_report.txt)

preprocess <- Previous Step | home | Next Step -> split