OxBS-seq data analysis

Prior to alignment paired-end reads were adaptor-trimmed and filtered using FASTQ Toolkit in Basespace (Illumina, San Diego, CA). End-trimming removed 3 bp and 4 bp from the 5′ and 3′ end of each paired-end read, respectively. Only reads with a Q score ≥ 25 and matching length criteria were used for mapping. Alignment of trimmed bisulfite-converted sequences was carried out using MethylSeq in Basespace (Illumina) (Krueger and Andrews 2011) against the mouse reference genome (GRCm38/mm10). Methylation call percentages for each CpG and non-CpG (CH) site within the genome were calculated by dividing the methylated counts over the total counts for that site in the oxidative bisulfite-converted libraries (OXBS). Genome-wide CpG and CH methylation levels were calculated separately. Hydroxymethylation levels in CpG (hmCG) and CH (hmCH) contexts were calculated by subtracting call levels from the oxidative bisulfite-converted (OXBS) libraries from the bisulfite-converted (BS) libraries. BAM files generated by MethylSeq (Basespace, Illumina) were run through MethylKit in R (Akalin et al. 2012) to generate context-specific (CpG/CH) coverage text files. Methylation percentages for the X and Y chromosomes were calculated separately as previously described [(Hadad et al. 2016)]. Briefly, coverage files were filtered to X and Y chromosomes to calculate sex-chromosome-specific methylation and hydroxymethylation levels in the same manner as genome-wide levels. Bisulfite conversion efficiency for C, mC, and hmC was estimated using CEGX spike-in control sequences (Supplemental Fig. 1). Untrimmed fastq files were run through CEGX QC v0.2 which output a fastqc_data.txt file containing the conversion mean for C, mC, and hmC.