PLAC-seq data analysis

PLAC-seq reads (paired-end, 50 bases) were aligned against the mm10 genome using BWA⁶⁴ -mem. Reads mapped to the same fragment were removed and PCR duplicate reads were removed using Picard MarkDuplicates. Filtered reads were binned at 10 kb size to generate the contact matrix. Individual bins that were overlapped with H3K4me3 peaks on TSSs were used for downstream differential contact analysis. For the peak calling, we used MAPS⁷⁵ with default settings and FitHiChIP⁷⁶ with coverage bias correction with default settings in 10 kb resolution.

For differential contact analysis, the raw contact counts in 10 kb resolution bins that have the same genomic distance were used as inputs for comparison. To minimize the bias from genomic distance, we stratified the inputs into every 10-kb genomic distance from 10 kb to 150 kb, and the other input bins with longer distances were stratified to have uniform size of input bins that were equal to that of 140–150 kb distance bins. Since each input showed negative binomial distribution, we used edgeR⁷⁰ to get the initial set of differential interactions. We only used bins that have more than 20 contact counts in each sample of two replicates for downstream analysis. The significances of these differential interactions are either due to the difference in their H3K4me3 ChIP coverage or 3D contacts coverage. Therefore, the chromatin contacts overlapping with differential ChIP-seq peaks (FDR < 0.05, logFC < 0.5) were removed and only the chromatin contacts with the same level of H3K4me3 ChIP-seq peaks were processed. In this differential analysis, we used all bins for inputs that included non-significant interactions that were not identified by MAPS or FitHiChIP peak caller, because the majority of short-range interactions were not identified as significant peaks due to their high background and the changes in the short-range interactions might be also critical for gene regulation. We identified a large number of differentially changed short-range interactions even though many of them were not identified as significant peaks, and we observed a clear correlation between these differentially changed interactions and the changes of active enhancer levels on their anchor regions during neural differentiation, suggesting these interaction changes might reflect the biological changes. We used significance level with change direction (−/+ log₁₀(p-value)) instead of fold change to show the changes of interactions, because fold change tends to be small value especially in short-range interactions even though the change is actually significant for biological aspects.