Read shifting.

For single-end ChIP-Sequencing data, CLUES generates a series of shift parameters from zero to the value inputted by the user (DNA fragment size of ChIP-Seq library).

1. CLUES first shifts mapped reads towards their 3' end by the first shift parameters (the smallest one) and uses the 5' positions of shifted reads to sort the reads in the genome. The distance between two neighboring reads is transformed into a bin (S17A Fig).
2. CLUES calculates the size distribution of all bins and employs a bin-length vs bin-number plot (S17B Fig) to identify the short bins, which are the bins on the left side of the peak (or the first peak) in the plot. It records the proportion of short bins.
3. CLUES iterates the above two steps using different shift parameters until CLUES gets the highest proportion of short bins (S17B Fig). Thus, CLUES acquires an optimized-shift parameter.
4. CLUES shifts all mapped reads with the optimized-shift-parameter and then converts them into bins.

The test on 22 Ctcf ChIP-Seq datasets shows the shift parameters estimated by CLUES are highly consistent with the ones of MACS2 (S17C Fig).

For paired-end ChIP-Seq data, CLUES does not shift mapped reads. Instead, CLUES uses the middle-point of a set of paired reads to build bins directly, and it estimates the optimized-shift-parameter as half of the median size of all paired-end reads.

ERs calling.

1. CLUES generates a series of maximum-allowed-length parameters from 0.2×optimized-shift-parameter to 2×optimized-shift-parameter. The first maximum-allowed-length parameter is 0.2×optimized-shift-parameter.
2. CLUES extracts reads from Chr1 to Chr5 in a ChIP-Seq data as the training set. It builds initial windows that include 20 bins from 5' to 3' until the last bin of a chromosome with one bin per step.
3. CLUES shrinks the initial window into step windows with one bin per step. Step windows with bins larger than or equal to the maximum-allowed-length parameter are excluded.
4. CLUES employs the Poisson model to calculate read enrichment in the step windows and calculates the final window from the step windows. It involves the following steps:
   1. (i). Count reads in a step window in a case sample (R_{window_case}) and the corresponding region in the control sample (R_{window_control}).
   2. (ii). Calculate reads background of the step window (R_{window_background}) as where d_genome is the average read density of a case sample; l_window is the length of a step window; SF is the scale factor of the number of mapped reads in case sample (R_case) and the number of mapped reads in control sample (R_control); l_genome is the length of the reference genome, and MCF is the correction factor to correct the mappability of repeat regions in a genome [45].
   3. (iii). Calculate p-value of reads enrichment in the step window by Poisson distribution function in R as
   4. (iv). Select the step window with the most significant p-value as the final window.
5. CLUES merges reads-enriched final windows (corrected p-value<0.05) to get merged windows set and calculates the length distribution of all merged windows in the training set.
6. CLUES iterates steps 2–5 using different maximum-allowed-length parameters until 80% of the merged windows are longer than 2×optimized-shift-parameter or the maximum-allowed-length parameter reaches 2×optimized-shift-parameter, whichever comes first (Fig 1F). Thus, CLUES acquires the optimized maximum-allowed-length.
7. CLUES calls ERs with the optimized-maximum-allowed-length.

SERs calling.

1. CLUES calculates the distances between two neighboring ERs. It regards two ERs as fragmented ERs if the distance is less than half of the length of either ERs (Fig 1B).
2. CLUES calculates the fragment rate by distance (FR-D) of ERs in a ChIP-Seq data. If FR-D>0.01 (default value, users can set their threshold), CLUES will call SERs.
3. CLUES shuffles ERs in the genome to get simulated ERs. It calculates the distances between neighboring simulated ERs, ranks them from small to large and records the distances located at the 1st to 10th percentile of the rank list as the values of the 1st to 10th percentile of the distance between neighboring peaks (PDNP). CLUES keeps the PDNP parameters that are smaller than 10 kb.
4. CLUES extracts reads from Chr1 to Chr5 as the training set. It builds initial windows that include x bins from 5' to 3' until the last bin of Chr1 with one bin per step. x is calculated as follows: where d_genome is the average read density of case sample. l_PDNP is the value of a selected PDNP parameter, which starts from the first PDNP.
5. CLUES shrinks the initial window into step windows with one bin per step. The step windows larger than or equal to the l_PDNP are excluded.
6. CLUES employs the Poisson model to calculate the read enrichment of the step windows and obtains a final window from the step windows the same way as ERs calling.
7. CLUES merges reads-enriched final windows (corrected p-value<0.05) into merged windows.
8. CLUES trims and filters the merged windows by ERs: for a merged window overlapping a single ER, CLUES keeps the ER as the merged window; for a merged window covering multiple ERs, CLUES keeps the shortest region covering all ERs as the merged window; CLUES discards the merged windows covering no ERs. CLUES calculates FR-D of the merged windows.
9. CLUES iterates steps 3–8 using l_PDNP with different PDNP parameters (1th to 10th PDNP) until FR-D reaches 1% or the FR-D reaches a minimum, whichever comes first (Fig 1G). Thus, CLUES acquires the optimized l_PDNP.
10. CLUES repeats steps 3–8 with the optimized l_PDNP to call merged windows as SERs.

LERs calling.

CLUES clusters input enrichment signals (either ERs or SERs) to be LERs. Here, we cluster ERs as an example.

1. CLUES calculates fold-change of read enrichment of the shortest region covering two neighboring ERs in the ChIP-Seq data. The fold-change value (RE) is calculated as Fig 1C shows. CLUES regards the neighboring ERs as fragmented ERs if their RE is larger than 1.5. CLUES calculates FR-RE of ERs in the data. If FR-RE>0.01, CLUES calls LERs in the data.
2. CLUES shuffles ERs in the genome to get simulated ERs. It calculates the distances between neighboring simulated ERs, ranks them from small to large and records the distances as the values of PDNP for every 5 percentile from the 5th to 50th percentile. CLUES keeps PDNP parameters that are smaller than 100 kb.
3. CLUES repeats steps 2–5 of ERs calling to call merged windows in the training set, starting l_PDNP with the 5th PDNP.
4. CLUES trims and filters the merged windows by ERs: for a merged window overlapping a single ER, CLUES keeps the ER as the merged window; for a merged window covering multiple ERs, CLUES keeps the shortest region covering all ERs as a merged window. CLUES calculates the FR-RE of ERs in the merged windows.
5. CLUES iterates steps 3–4 above using l_PDNP with different PDNP parameters until the FR-RE decreases to 99% or l_PDNP reaches a maximum, whichever comes first (Fig 1H). Thus, CLUES acquires an optimized l_PDNP.
6. CLUES runs the above steps 3–4 with the optimized l_PDNP to call merged windows as LERs.

Replace ERs with SERs, CLUES uses the same way to call LERs from SERs.

Reporting ERs, SERs, and LERs.

CLUES reports length, reads number, p-value, q-value and fold-change value for enrichment signals of ERs, SERs, and LERs, and reports summits for ERs as follows:

1. CLUES calculates the length of an enrichment signal as l_ER.
2. CLUES counts the read number of the enrichment signal in the case-sample as R_{ER_case}.
3. CLUES calculates the p-value of the enrichment signal as follows:
   1. (i). Calculate the read background of the enrichment signal (R_{ER_background}) as here, where d_genome is the average read density of the case sample; d_ER is the read density of the corresponding region of the enrichment signal in control sample; d_1kb and d_10kb are the read density in the corresponding 1 kb and 10 kb regions around the enrichment signal in the control sample; l_genome is the length of the reference genome; MCF is the correction factor to correct the mappability of repeat regions in a genome[45]; SF is the scale factor of the number of mapped reads in the case sample (R_case) and the number of mapped reads in the control sample (R_control).
   2. (ii). Calculate the p-value of the read enrichment of the enrichment signal by a Poisson distribution function in R language as:
4. CLUES calculates the q-value from the p-value of the called enrichment signals by Benjamini & Yekutieli method [46].
5. CLUES calculates the fold-change (FC_ER) value of the called enrichment signals as: where and R_{ER_background} are calculated in the same way as CLUES did in the third step.
6. CLUES calculates the summit of an ER as follows: each read position is extended optimized-shift-parameter bases from its center for the reads in the ER; the location with the highest fragment pileup is referred to as ‘the summit’.

Ranking of ERs and broad E-signals identified by the methods.

The commands for calling ERs and broad E-signals by the methods are summarized in S3 Table.

We used a q-value of 0.05 to run CLUES and a default q-value (0.05) to run MACS2, a default q-value (0.001) to run MUSIC and a default FDR (0.0001) to run PeakRanger. Because there was no adjusted P-value option available, we used a default P-value (0.001) to run SISSRs. We ranked ERs identified by MACS2, MUSIC, PeakRanger by their adjusted P-value in ascending order and the ERs identified by SISSRs by their P-value in ascending order (for the ERs with equal adjusted-P-value or P-values, we ranked them by their read fold-change value in descending order). To rank the ERs identified by CLUES, we first ranked their q-value in ascending order and recorded the rank of each ER as Rank_{ER_qvalue}; then, we ranked the ERs again by the summit height of each ER in descending order and recorded the rank of each ER as Rank_{ER_summit}; finally, we ranked the ERs by the sum value of Rank_{ER_qvalue} and Rank_{ER_summit} of every ER in ascending order.

To identify broad H3K4me3 E-signals, we set the q-value to 0.05 to run CLUES, then we used a default q-value (0.1) to run MACS2 and a default q-value (0.001) to run MUSIC. We ranked the broad H3K4me3 E-signals identified by MACS2 and MUSIC by their q-value in ascending order (for broad H3K4me3 E-signals with equal q-values, we ranked them by their read fold-change value in descending order). We ranked the broad H3K4me3 E-signals identified by CLUES by their read number in ascending order.

To identify broad E-signals of H3K27me3 and H3K36me3 by CLUES, PeakRanger, SICER, and MUSIC, we set the q-value at 0.05 to run CLUES, the FDR at 0.05 to run SICER, a default q-value (0.001) to run MUSIC and a default FDR (0.0001) to run PeakRanger. We ranked the broad E-signals identified by SICER, MUSIC, and PeakRanger by their adjusted P-value in ascending order (for the broad E-signals with equal adjusted P-value, we ranked them by their read fold-change value in descending order). To rank the broad E-signals identified by CLUES we first ranked their fold-change value in descending order and recorded the rank of each broad E-signal as Rank_{LER_enrichment}; then we ranked the broad E-signals by their length in descending order and recorded the rank of each broad E-signals as Rank_{LER_length}; finally, we ranked broad E-signals by the sum value of Rank_{LER_enrichment} and Rank_{LER_length} of every broad E-signal in ascending order.

Comparing the PPV of CLUES with the other methods in detecting reliable-ERs.

We defined ERs with the corresponding Ctcf or Nrsf motifs within 150 bp around their summits as reliable-ERs [8]. We compared the PPV of CLUES to the other methods in detecting reliable-ERs in the ChIP-Seq data of Ctcf and Nrsf.

We conducted the comparison using the following steps:

1. We identified ERs with the various methods and ranked the ERs as described in “Ranking of ERs and broad E-signals identified by the methods”.
2. We downloaded Ctcf and Nrsf motifs from the MEME database [47] (S10 Table) and scanned the corresponding motifs 150 bp upstream and downstream of the summits of the ERs by FIMO with a default p-value (0.0001) [48].
3. We calculated reliable-ERs rated in the top 100, 200, and 300 …, ERs of each method and fitted their reliable-ERs-rate curves.
4. We calculated the areas bounded by the X axis, Y axis and the reliable-ERs-rate curves of the same number of top-ranked ERs identified by CLUES and the other methods being compared.
5. We divided the area of a rival method (MACS2, MUSIC, PeakRanger or SISSRs) by the area of CLUES. The Ratio > 1.01 indicates CLUES has lower PPV than the rival method ("L" in Fig 2A); Ratio < 0.99 indicates CLUES has higher PPV than the rival method ("H" in Fig 2A); 0.99≤ Ratio ≤ 1.01 indicates CLUES has equal PPV with the rival method ("E" in Fig 2A). The ratios of the rival methods against CLUES for 62 Ctcf and Nrsf ChIP-Seq data are summarized in S4 Table.

Comparing the number of reliable-ERs identified by CLUES and MACS2 under different signal-to-noise (SNR) values.

In each ChIP-Seq data, CLUES shifts reads according to the estimated shift parameter, then it calculates the size distribution of all bins and employs a bin-length vs. bin-number plot (S17B Fig) to identify the short bins, which are the bins on the left side of the peak (or the first peak) in the plot. It records the proportion of short bins. We used the proportion of short bins to assess the signal-to-noise ratio (SNR) of ChIP-Seq datasets at the same mapping depth. We calculated SNR of Ctcf and Nrsf ChIP-Seq data as follows:

1. We sampled 5 million reads from each data, then we calculated the proportion of small bins in the data as the value of SNR.
2. We calculated the ratio (Ratio_C/M) of the number of reliable-ERs identified by CLUES and MACS2 under default q<0.05 in each data set.
3. We associated the SNRs with Ratio_C/M between biological replicate datasets.

GO analysis.

At first, we associated top 100 and 1000 SERs and LERs in each data with genes as follows:

1. We downloaded gene location information from UCSC Genome Browser (hg19 version and mm9 version) [49].
2. A gene overlapped by or located within 1 kb of an SER was recognized as an SER-associated gene. A gene with 80% of its length covered by an LER was recognized as an LER-associated gene.

Then, we took the SERs-associated genes or LERs-associated genes as target set. We took all human (or mouse) genes as the background set. We employed the hypergeometric distribution test as enrichment function. We used GOstats to conduct GO analysis on the associated genes [50]. The database was from org.Hs.eg.db and org.Mm.Hs.eg.db, and the test method was the GOHyperGParams method. GO terms with p<10⁻⁵ are reported.

Calling Nanog LERs and Oct4 LERs in mouse ER cells.

It is known that the pluripotency and self-renewal of mouse ES cells are potentiated by core TFs Oct4, Nanog and Sox2, and their binding sites frequently cluster in close genomic proximity to activate the expression of downstream genes [25, 51–53]. Because the H3K4me3 modification is required for activating expression [54, 55], the binding clusters of Oct4 and Nanog should couple with H3K4me3 SERs. Therefore, we used the LERs calling module to learn the connecting length (l_PDNP) that allows the LERs of Nanog and Oct4 to have the most similar length compared to H3K4me3 SERs

1. We called ERs and SERs in H3K4me3 ChIP-Seq data by CLUES.
2. We called ERs and SERs in Nanog ChIP-Seq data by CLUES.
3. We called Nanog LERs from Nanog SERs as following:
   1. (i). We conducted steps 2–4 in “LERs calling” of the CLUES algorithm.
   2. (ii). We iterated steps 2–4 in “LERs calling” using l_PDNP with different PDNP parameters and recorded the l_PDNP where FR-RE > 99%.
   3. (iii). We called different Nanog LERs sets with the recorded l_PDNP.
   4. (iv). We calculated length dissimilarity between Nanog LERs and H3K4me3 SERs as follows: the dissimilarity score is the area bounded by the length distribution curves of Nanog LERs and the H3K4me3 SERs.
   5. (v). We selected Nanog LERs from different Nanog LERs sets and calculated the above dissimilarity score. We selected the Nanog LERs set with the minimum dissimilarity score as the optimized Nanog LERs set and outputted the LERs.

We similarly called Oct4 LERs.

Integrating H3K4me3 SERs, H3K27me3 LERs, Oct4 LERs and Nanog LERs in mouse ES cells.

We integrated all H3K4me3 SERs, all H3K27me3 LERs, the top 1000 Oct4 LERs (~ top 5% LERs) and the top 1000 Nanog LERs (~top 5% LERs) by their locations on the genome. A Nanog LER (or an Oct4 LER) overlapped by or located within 10 kb of an H3K4me3 SER was recognized as an H3K4me3-associated LER. If more than one Nanog (or Oct4) LERs were associated with an H3K4me3 SER, the LER with the highest rank was picked. An H3K27me3 LER overlapping by or located within 1 kb of an H3K4me3 SER was recognized as an H3K4me3-associated LER.

Cell culture.

V6.5 ESCs (the F1 hybrid of 129SvJae/C57BL/6) were cultured on gelatin-coated (Millipore) plates in Dulbecco’s modified Eagle’s medium (DMEM) supplemented with 15% FBS, 0.1 mM b-mercaptoethanol, 2 mM L-glutamine, 0.1 mM nonessential amino acid, 1000 U/ml recombinant leukemia inhibitory factor (LIF; Millipore), and 30 U/ml penicillin/streptomycin. The media were changed daily, and ESCs were split every two days. HEK293T cells were cultured in DMEM with 10% FBS and 1000 U/ml penicillin/streptomycin.

Whole-Transcriptome shotgun sequencing (RNA-Seq) and data processing.

Total RNA was isolated from 5x10⁶ mouse ES cells with the RNeasy Plus Mini Kit (QIAGEN, 74134). Poly(A)-containing mRNA molecules were purified using 1 μg total RNA as the starting material, and the library was built following NEB. Next, the Ultra RNA Library Prep Kit was used to obtain the Illumina protocols (NEB, 7530L). The library was sequenced on an Illumina HiSeq 2500 Machine. Reads were mapped to the reference genome of mm9 using TopHat with default settings [56].

Lentivirus production and transduction.

Lentivirus was produced through the co-transfection of the lentiviral vectors (Mouse GeCKO v2 Library from Addgene; the LentiCRISPRv2 plasmid was used for single-gene targeting) with the envelope plasmid (psPAX2) and VSV-G packaging plasmids (pMD2G; 4:3:3) into HEK293T cells using Lipofectamine 2000 according to the manufacturer's instructions. Media were changed 24 hours after transfection. The virus-containing supernatant was collected and filtered through a 0.45 μm low protein-binding membrane (SLHV033RB, Millipore) 48 hours after transfection. Transduction was performed on 24-well plates with 5x10⁵ mouse ES cells in each well. The spin-infection was performed in medium containing 8 μg/mL of polybrene at 1,800 rpm for 45 min at room temperature. Transduced cells were selected 24 hr after spin-infection under 1 μg/ml puromycin for 3 days.

Targeted next-generation sequencing.

Genomic DNA was extracted by the phenol/chloroform method from 10 million mouse ES cells. Sequencing libraries were acquired by two rounds of PCR amplification: round 1 served to amplify the targeted locus with 23 cycles [57], and round 2 served to add P5 and P7 adapters with the index for each sample with ten cycles [57]. Amplicons from the second round of PCR were sequenced on a HiSeq 2500 (Illumina). Primer sequences are shown in the S11 Table.

Genomic-wide CRISPR/Cas9 negative selection genetic screen.

Lentivirus of GeCKO library was transduced into mouse ES cells as described above with MOI = 0.3. After puromycin selection, 10 million cells were seeded into a 10-cm dish every generation to maintain 166x coverage of the library. A total of 10 million cells from the first generation and 10th generation were used for genomic DNA extraction for targeted sequencing of sgRNA. Additionally, 50x50 μl first-round PCR reactions were performed with 825 ng genomic DNA in each reaction to achieve 100X coverage. A 2 μl volume of mixed first-round PCR products were used as templates for the second round PCR to generate targeted sequencing library.

Prioritizing genes that contributed to ES cells proliferation from a genomic-wide negative selection genetic screen.

We prioritized genes that may contribute to mouse ES cells proliferation by the following steps:

1. We counted reads of sgRNAs in the sgRNAs pool after targeted next-generation sequencing; we only counted the reads covering 100% of the sgRNA sequence.
2. We filtered the non-targeting sgRNAs (NT-sgRNAs) in the GeCKO v2 Library and selected the sgRNAs with reads>10 at the 1st generation (a total of 874 sgRNAs) as a control pool.
3. We ranked the NT-sgRNAs in the 1st generation by their read number in descending order. We ranked them at the 10th generation in the same way.
4. We compared the reads of every sgRNA with the reads of NT-sgRNAs in the 1st generation; we calculated the corresponding ranking position of every sgRNA inside the ranked NT-sgRNA. We calculated the corresponding ranking position of every sgRNA in the 10th generation in the same way.
5. We calculated the difference in ranking position of every sgRNA between the 10th generation and 1st generation.
6. We selected the sgRNA which has the highest value of the difference from the sgRNAs set of a target gene. We used the value to rank the targeted genes in descending order.

Alkaline phosphatase (AP) Staining.

ESCs transduced with lenti-sgRNA were first selected under 1 ug/ml puromycin for three days. AP staining was performed according to the manufacturer’s instructions for the Alkaline Phosphatase Detection Kit (SCR004, Millipore).

Counting indel types of a targeted locus by CRISPR/Cas9.

FastQ reads of targeted NGS sequencing were mapped to the target locus using Bowtie2 with default parameters apart from an adjustment to relax the gap extension penalty (option:—rdg 5,1) [58]. The CIGAR string was extracted from SAM format files for frequency counting of indel types. Primer sequences are listed in the S11 Table.

Vectors and mutagenesis.

The lenticrispv2 plasmid was purchased from addgene; cDNA of Zmynd8 (MC202682) and Fam60a (MC203541) were purchased from Origene and subcloned into pyCAGIP with an additional Flag tag sequence at the C-terminus; Abt1 cDNA was directly cloned from mouse ES cells and subcloned into pyCAGIP with a Flag tag at the C-terminus.

Silent mutations against sgRNA targeting were acquired by Q5 Site-Directed Mutagenesis Kit (NEB) according to the manufacturer’s instructions on pyCAGIP[59].

The sequences of the mutant sites are as follows:

Abt1-copy2:CTACTCGGCCAAATTCCAGTTGG

Fam60a-copy2: CTACAGTAACCAGTCGGACGAAG

Zmynd8-copy2: ATTCCAGAAGCCTGTCCCCTTAG

In Vitro CRISPR/Cas9 cleavage assay.

sgRNAs were expressed from the LentiCRISPv2 plasmid according to the manufacturer’s instructions for the HiScribe T7 High Yield RNA Synthesis Kit (NEB #E2050). Cas9 Nuclease (M0386, NEB) was used to perform in vitro cleavage on the same amount of input plasmid. After overnight digestion, samples were resolved on a 0.8% agarose gel with 10 kb DNA Ladder (B600032, Sangon Biotech).

In vivo CRISPR/Cas9 cleavage assay.

The same amount of plasmids were transiently transfected into wide-type and CRISPR/Cas9-expressing mouse ES cells. At 24 hours after transfection, the cells were harvested and subjected to Western blotting analyses.

Antibodies for Western blotting analyses.

Cells were harvested and lysed in SDS-PAGE sample buffer. Equal amounts of total protein were loaded in each lane. Proteins were resolved by SDS-PAGE, transferred to 0.45 μm PVDF membranes, and probed with the indicated antibodies. The antibodies used for Western blotting were as follows: anti-α-tubulin (Sigma DM1a) and anti-Flag (a gift from Huang Lab).

Narrow down gene prioritization list of the genomic-wide negative selection genetic screen by SICER and MUSIC.

To search the minimum candidate genes, including Zmynd8, Abt1, and Fam60a, from the SICER results, we filtered genes with the following criteria: 1. marked by broad H3K4me3 E-signals; 2. marked by top 11899 Nanog or top 11899 Oct4 broad E-signals; 3. not marked by broad H3K27me3 E-signals. Then, we overlapped the candidate genes with the gene prioritization list of the genomic-wide negative selection genetic screen. We selected genes with higher or equal ranks compared to Zmynd8, Abt1, and Fam60a in the gene list. These were the genes narrowed down by SICER from the genomic-wide negative selection genetic screen.

We used the same strategy to get the genes narrowed down by MUSIC from the gene prioritization list of genomic-wide negative selection genetic screen. However, the genes were marked by top 2043 Nanog or top 2043 Oct4 broad E-signals and not by top 11899 Nanog or top 11899 Oct4 broad E-signals.