ATAC- and ChIP-seq pre-processing

Raw reads were trimmed with Trimmomatic (v0.36) and aligned to strain-specific genome using STAR (v2.5.4b) with ‘–alignIntronMax 1′ to turn off splice awareness and –alignEndsType ‘EndToEnd’ to prevent soft clipping. For ATAC-seq the reads were filtered for true pairs and MAPQ > 30 with Samtools (v1.9) (42). For ChIP-seq the reads were filtered for MAPQ > 20. For both, duplicates were removed after filtering using Picard (v2.8.13, https://Broadinstitute.github.io/picard/) MarkDuplicates. ATAC-seq peaks were called with HMMRATAC (v1.2.9) (12) using ‘–means 75200400600 –upper 20 –lower 10′. H3K27ac, CTCF and TCF7L2 ChIP-seq peaks were called with MACS2 (v2.2.7.1) using ‘-g mm -B –call-summits –keep-dup auto’, and ‘-q 0.05′ for H3K27ac, ‘-q 0.01′ for TCF7L2, and ‘-q 0.001′ for CTCF (43). After peak calling, summits were widened to peaks with Bedtools (v2.27.1) using ‘slop -b 75′ to prevent the summits of the 129 strain that overlapped the B6 deletions to be lost during the coordinate conversion (44). Consensus peak sets for ATAC-seq and TCF7L2 ChIP-seq were formed from the peak centres using kernel density estimation -based approach adapted from Tuoresmäki et al. (45). The consensus peak set for CTCF was created by merging the narrowPeaks from MACS2 using Bedtools merge. For H3K27ac, we used ChIP-R (v1.2.0) to identify group-wise the regions marked by H3K27ac (46). The consensus peak set for H3K27ac was created by merging the group-specific regions using Bedtools merge. We also generated harmonized H3K27ac summits using kernel density estimation -based approach and selected the summits overlapping group-specific H3K27ac regions for valley-based active nucleosome free region (NFR) detection. NFRs between two H3K27ac summits (±1500 bp) were defined as active. Peaks were filtered for the ENCODE blacklisted regions (47). For both ATAC and ChIP-seq, peaks and filtered BAM-files for the 129-strain were shifted to B6 coordinates using the g2gtools ‘convert’ function.