Sequence reads were demultiplexed, quality filtered with FastQC (Andrews, 2010) and mapped to the human genome (hg19) using Bowtie 1.2 (Langmead et al., 2009). We used option “-m 1” to retain those reads that map only once to the genome. Each individual sample contributed with the same number of reads in the ChIP-seq final merged sample. For the computation of ChIP-seq binding sites (peaks), we used MACS2 (Zhang et al., 2008) with option “-q 0.01.” We used the R package DESeq2 (Love et al., 2014) to identify differentially expressed genes from RNA-seq data following authors guidelines. Only genes with associated adjusted p values ≤ 0.05 and absolute fold change ≥ 2 were considered as differentially expressed.For the identification of Pol II differentially bound genes, we followed the same approach than for differentially expressed genes, this time by restricting the gene length to the region stretching 2kb from 500 bp downstream the TSS.

To estimate the level of Pol II recruitment, we used the so-called “pausing index” (PI) and “pause release ratio” (PRR). PI is defined as the ratio of Pol II enrichment within the promoter to that in the gene body, while PRR is an inverse of PI that restricts the gene body to the first 2kb downstream the TSS. For the estimation of both parameters, we based our strategies on Chen et al. (2015) and Day et al. (2016). To calculate PI, the level of Pol II within the promoter was computed as the sum of ChIP-seq reads in 400 bp surrounding the TSS. The level of Pol II within the gene was computed as the average number of reads in 400 bp windows throughout the gene body, from 200 bp downstream the TSS. Finally, the PI was estimated as the ratio of both measures. To calculate PRR, Pol II level within the promoter was estimated as the sum of ChIP-seq reads from 100 bp upstream to 300 bp downstream of the TSS. Gene body Pol II level was computed as the sum of reads within the region stretching from 300 bp downstream of the TSS to 2 kb downstream of the TSS. After normalizing each value by the corresponding window sizes, the PRR was estimated as the level of Pol II within the gene body divided by the level of Pol II within the promoter.

For Gene-filtering, we started from the whole set of protein-coding transcripts associated to Ensembl-annotated genes (GRCh37, release 75), from which we kept only those having a peak of Pol II (see peak calling section) overlapping with the region from 0 to 500 bp downstream of the TSS. If several transcripts of the same gene were found to match this condition, the one whose TSS was closest to the Pol II peak was selected. In order to account for potential false negatives when calling Pol II peaks, some genes with no associated Pol II peak were also selected if the window from 0 to 500 bp downstream of the TSS was found to: 1) have a value of Pol II reads per million (RPM) larger than 0.5, or 2) have a value longer than 0.5 in the difference of H3K4me3 RPM and the corresponding input. Finally, genes that were closer than 1kb of another gene or were smaller than 2kb were excluded. We ended up keeping 9,588 human genes (∼42% of total).

Regulatory regions, namely enhancers, promoters, and insulators, were defined as follows. For promoters, the whole set of transcripts associated to Ensembl-annotated genes were considered. Then, promoters were defined as ± 1 kb from the TSSs. Enhancers were defined as H3K27ac peaks not overlapping with a promoter, and insulators as CTCF peaks not overlapping with promoters and enhancers.

ChIP-seq averaged reads around TSSs were computed using the R package bamsignals (Mammana and Helmuth, 2019) and smoothed using a Gaussian smoothing kernel with the R function ksmooth, respectively. To generate the profile of TOP2A ChIP-seq signal around randomized genes, the 148 upregulated genes upon 2h of merbarone treatment were first considered. Then, 50 sets of 148 genes were randomly selected from the human genome and TOP2A ChIP-seq reads were counted around the TSSs. Finally, the median of such read counts was smoothed and plotted.

Publicly available sequencing data used in this study include ChIP-seq of several proteins and post-translational modifications: CTCF (ENCSR000DVI), H3K4me3 (ENCSR000DVK) and H3K27ac (Sanchez et al., 2018). CTCF and H3K4me3 BAM files (hg19) were batch-downloaded from ENCODE. H3K27ac and H3K27me3 raw sequencing files were processed as described above (Langmead et al., 2009).

The R package TopGO (Alexa and Rahnenfuhrer, 2019) was used to calculate the significance of Gene Ontology (GO) terms associated to differentially expressed genes after merbarone treatment. We computed enrichments using the Fisher’s exact test and the default algorithm (weight01), which is a hybrid between the ‘elim’ and the ‘weight’ algorithms described (Alexa et al., 2006). To perform hypergeometric-based tests, there is a need for defining a ‘gene universe’ (which can be conceptualized as the number of balls in an urn) and a set of ‘interesting genes’ from that universe. To define the gene universe, we started from expressed genes, which were defined as genes for which the sum of RNA-seq reads (combining the replicates) overlapping exons was larger than 10. Then, the gene universe was defined as expressed genes mapping to at least one GO term and the set of interesting genes as differentially expressed genes upon merbarone treatment.

Note: The content above has been extracted from a research article, so it may not display correctly.



Q&A
Please log in to submit your questions online.
Your question will be posted on the Bio-101 website. We will send your questions to the authors of this protocol and Bio-protocol community members who are experienced with this method. you will be informed using the email address associated with your Bio-protocol account.



We use cookies on this site to enhance your user experience. By using our website, you are agreeing to allow the storage of cookies on your computer.