extfieldsofstudy:
Biology
Computer Science
Medicine

text:
Practical Guidelines for the Comprehensive Analysis of ChIP-seq Data

Mapping the chromosomal locations of transcription factors, nucleosomes, histone modifications, chromatin remodeling enzymes, chaperones, and polymerases is one of the key tasks of modern biology, as evidenced by the Encyclopedia of DNA Elements (ENCODE) Project. To this end, chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-seq) is the standard methodology. Mapping such protein-DNA interactions in vivo using ChIP-seq presents multiple challenges not only in sample preparation and sequencing but also for computational analysis. Here, we present step-by-step guidelines for the computational analysis of ChIP-seq data. We address all the major steps in the analysis of ChIP-seq data: sequencing depth selection, quality checking, mapping, data normalization, assessment of reproducibility, peak calling, differential binding analysis, controlling the false discovery rate, peak annotation, visualization, and motif analysis. At each step in our guidelines we discuss some of the software tools most frequently used. We also highlight the challenges and problems associated with each step in ChIP-seq data analysis. We present a concise workflow for the analysis of ChIP-seq data in Figure 1 that complements and expands on the recommendations of the ENCODE and modENCODE projects. Each step in the workflow is described in detail in the following sections.

Introduction to ChIP-seq Technology
Chromatin immunoprecipitation followed by sequencing (ChIP-seq), first described in 2007 [1][2][3][4], allows in vivo determination of where a protein binds the genome, which can be transcription factors, DNA-binding enzymes, histones, chaperones, or nucleosomes. ChIP-seq first cross-links bound proteins to chromatin, fragments the chromatin, captures the DNA fragments bound to one protein using an antibody specific to it, and sequences the ends of the captured fragments using next-generation sequencing (NGS). Computational mapping of the sequenced DNA identifies the genomic locations of bound DNA-binding enzymes, modified histones, chaperones, nucleosomes, and transcription factors (TFs), thereby illuminating the role of these protein-DNA interactions in gene expression and other cellular processes. The use of NGS provides relatively high resolution, low noise, and high genomic coverage compared with ChIP-chip assays (ChIP followed by microarray hybridization). ChIP-seq is now the most widely used procedure for genome-wide assays of protein-DNA interaction [5], and its use in mapping histone modifications has been seminal in epigenetics research [6].

The Analysis of ChIP-seq Data Sequencing Depth
Effective analysis of ChIP-seq data requires sufficient coverage by sequence reads (sequencing depth). The required depth depends mainly on the size of the genome and the number and size of the binding sites of the protein. For mammalian transcription factors (TFs) and chromatin modifications such as enhancer-associated histone marks, which are typically localized at specific, narrow sites and have on the order of thousands of binding sites, 20 million reads may be adequate (4 million reads for worm and fly TFs) [7]. Proteins with more binding sites (e.g., RNA Pol II) or broader factors, including most histone marks, will require more reads, up to 60 million for mammalian ChIP-seq [8]. Importantly, control samples should be sequenced significantly deeper than the ChIP ones in a TF experiment and in experiments involving diffused broad-domain chromatin data. This is to ensure sufficient coverage of a substantial portion of the genome and non-repetitive autosomal DNA regions. To ensure that the chosen sequencing depth was adequate, a saturation analysis is recommended-the peaks called should be consistent when the next two steps (read mapping and peak calling) are performed on increasing numbers of reads chosen at random from the actual reads. Saturation analysis is built into some peak callers (e.g., SPP [9]). If this shows that the number of reads is not adequate, reads from technical replicate experiments can be combined. To avoid over-sequencing and estimate an optimal sequencing depth, it is important to take into account library complexity. Several tools are available for this purpose. For example, the preseq package allows users to predict the number of redundant reads from a given sequencing depth and how many will be expected from additional sequencing [10]. Similarly, the ENCODE software tools offer a quality metric called the PCR bottleneck coefficient (PBC), defined as the fraction of genomic locations with exactly one unique read versus those covered by at least one unique read.

Read Mapping and Quality Metrics
Before mapping the reads to the reference genome, they should be filtered by applying a quality cutoff (Box 1). The remaining reads should then be mapped using one of the available mappers such as Bowtie [11], BWA [12], SOAP [13], or MAQ [14]. Recent versions support gapped alignment (e.g., Bowtie2), but detection of indels is not necessary for most ChIP-seq experiments. It is important to consider the percentage of uniquely mapped reads reported by the mapper. The percentage varies between organisms, and for human, mouse, or Arabidopsis ChIP-seq data, above 70% uniquely mapped reads is normal, whereas less than 50% may be cause for concern. A low percentage of uniquely mapped reads often is due either to excessive amplification in the PCR step, inadequate read length, or problems with the sequencing platform, but with some ChIPed proteins it may be unavoidable (e.g., if the protein binds frequently in repetitive DNA). The read mappers are designed to allow a (user-settable) number of mismatches in the reads, and it is important to choose this parameter to be appropriate with the NGS platform being used (consult the manufacturer). A final potential cause of high numbers of ''multimapping'' reads is that the protein binds frequently in regions of repeated DNA. In this last case, using paired-end sequencing to reduce the mapping ambiguity may help. It should be kept in mind that multi-mapping reads will be ignored (filtered out) by most peak-calling algorithms (see section ''Peak Calling''), although they can drive the discovery of novel binding sites [15].
After mapping, the signal-to-noise ratio (SNR) of the ChIP-seq experiment should be assessed, for example via quality metrics such as strand cross-correlation [7] or IP enrichment estimation using the software package CHANCE [16] (Box 2). These measures will detect several possible failure modes of ChIP-seq: insufficient enrichment by immunoprecipitation step, poor fragment-size selection, or insufficient sequencing depth. Strand cross-correlation analysis is built into some peak callers (e.g., SPP or MACS [17] [version 2]).

Peak Calling
A pivotal analysis for ChIP-seq is to predict the regions of the genome where the ChIPed protein is bound by finding regions with significant numbers of mapped reads (peaks). A fine balance between sensitivity and specificity depends on choosing an appropriate peak-calling algorithm and normalization method (Boxes 3-6, Table S1, and [18,19]) based on the type of protein ChIPed: point-source factors such as most TFs (Box 3), broadly enriched factors such as histone marks (Box 4), and those with both characteristics such as RNA Pol II (Box 5) [20]. It is strongly recommended that mapped reads from a control sample be used (e.g., from input DNA), although some peak callers can use GC content or mappability as information necessary to assess the level of non-specific or background binding. Duplicate reads (same 59 end) can be removed before peak calling to improve specificity (Box 7). Although some peak callers support both single and paired-end reads (e.g., MACS), others

Box 1. Quality metrics of sequence reads
Preprocessing of ChIP-seq data will, in general, be similar to that of any other sequencing data and will assess the quality of the raw reads to identify possible sequencing errors or biases (FastQC can be used for an overview of the data quality). Phred quality scores are used to describe the confidence of each base call in each sequence tag, are logarithmically linked to error probabilities, and can be used to filter low-quality reads. After this filtering step, it may also be necessary to trim the end of reads that are of low quality (see sickle, https://github.com/ najoshi/sickle). Additionally, library complexity is a common quality measure for ChIP-seq libraries (preseq package [10] or PCR bottleneck coefficient [PBC] from ENCODE tools, https://code.google.com/p/phantompeakqualtools/), and library complexity is linked to many factors such as antibody quality, over-cross-linking, amount of material, sonication, or over-amplification by PCR. The latter can be corrected by systematic identification and removal of redundant reads, which is implemented in many peak callers as it may improve their specificity. Readers may be interested in the Galaxy toolbox, which offers access to many of the tools described here [50].
are specifically designed to improve sensitivity and specificity in paired-end sequencing (e.g., SIPeS [21]). Existing peak callers have many user-settable parameters that can greatly affect the number and quality of the peaks called. For instance, the enrichment metric for most peak callers, such as p-value or FDR, could be hugely affected by the statistical model used, the sequencing depth, or the actual number of binding sites in the genome. Thus, using the same p-value or FDR threshold does not ensure that the numbers of peaks called are comparable across libraries and different peak callers [22]. A better approach is to threshold the irreproducible discovery rate (IDR) [23], which, along with motif analysis, can also aid in choosing the best peak-calling algorithm and parameter settings (see sections ''Assessment of Reproducibility'' and ''Motif Analysis'').

Assessment of Reproducibility
To ensure that experimental results are reproducible, it is recommended to perform at least two biological replicates of each ChIP-seq experiment and examine the reproducibility of both the reads and identified peaks [7,24]. The reproducibility of the reads can be measured by computing the Pearson correlation coefficient (PCC) of the (mapped) read counts at each genomic position [25]. The range of PCC is typically from 0.3-0.4 (for unrelated samples) to .0.9 (for replicate samples in high-quality experiments). Low values typically suggest one or both replicates may be of low quality. However, this quantity can be dominated by a small number of very highly enriched regions, so it may not reflect the reproducibility for regions that are less enriched [26]. Thus it is important to remove the artefact regions with high ChIP signals, such as regions near centromeres, telomeres, satellite repeats, and ENCODE and 1000 Genomes blacklisted regions, before computing the PCC. To measure the reproducibility at the level of peak calling, IDR analysis (Box 8) [23] can be applied to the two sets of peaks identified from a pair of replicates. This analysis assesses the rank consistency of identified peaks between replicates, and outputs the number of peaks that pass a user-specified reproducibility threshold (e.g., IDR = 0.05). It has been reported that using a reproducibilitybased metric (e.g., IDR) rather than an enrichment-based metric (e.g., FDR or pvalue) makes the numbers of peaks declared more comparable across experiments [7]. In addition, IDR analysis can also be used for comparing and selecting peak callers [8,23] and identifying experiments with low quality [7].

Differential Binding Analysis
Comparative ChIP-seq analysis of an increasing number of protein-bound regions across conditions or tissues is expected with the steady raise of NGS projects. For example, temporal or developmental designs of ChIP-seq experiments can provide different snapshots of a binding signal for the same TF, uncovering stage-specific patterns of gene regulation [27,28]. With this in mind, one should note that the simple binary overlap of two sets of peaks (e.g., BEDTools [29]) does not represent the optimal approach when comparing peaks [25].
Two alternatives have been proposed. The first one-qualitative-implements hypothesis testing on multiple overlapping sets of peaks [30], therefore extending the two-set overlap approach mentioned above. The second one-quantitativeproposes the analysis of differential binding between conditions based on the total counts of reads in peak regions or on the read densities, i.e., counts of reads overlapping at individual genomic positions (Table S3 and [31,32]). The direct calculation of differentially bound regions

Box 2. Quality metrics of read counts
Strand cross-correlation analysis [7] assesses data quality by measuring the degree of immunoprecipitated (IP) fragment clustering in ChIP-seq experiments. It is developed based on the observations that (1) a high-quality ChIP-seq experiment often shows a significant clustering of enriched DNA sequence tags at the locations bound by the protein of interest, and that (2) the enriched sequence tags on the forward and reversed strands are positioned at a distance from the binding site center that depends on the fragment size distribution [9]. This method quantifies the degree of clustering by computing the cross-correlation between the two strands, i.e., the Pearson correlation between the strand-specific read density profiles as a function of the shift (k) applied to one of the two strands ( Figure S1). The cross-correlation typically peaks at the shift corresponding to the fragment length and the shift corresponding to the read length. The ratio between the cross-correlation at the fragment length and the background crosscorrelation, referred to as normalized strand cross-correlation coefficient (NSC), and the ratio between cross-correlation at the fragment length and the crosscorrelation at the read length, referred to as relative strand cross-correlation coefficient (RSC), jointly reflect signal-to-noise ratio in the ChIP-seq data. Very successful ChIP experiments generally have NSC.1.05 and RSC.0.8 [7], although there can still be significant biological information present in ChIP-seq data not meeting these criteria. Readers may refer to [7] for prototypical profiles of crosscorrelation illustrated on ENCODE data.
The software CHANCE [16] assesses IP strength by estimating and comparing the IP reads pulled down by the antibody and the background, using a method called signal extraction scaling [77]. For each sample, it first bins the genome into nonoverlap bins both for the IP and the Input, then partitions the bins into a signal region and a background region by comparing the cumulative distributions of tag counts in the bins of the IP and the Input. It next computes a p-value for significance of enrichment according to the percentage allocation of reads in each type of regions. Based on the empirical p-value distribution computed from a set of ENCODE IP-Input and Input-input experiments on human data, it estimates a q-value by treating the two types of experiments as true positives and false positives, respectively. The q-value thus is interpreted as the fraction of comparisons with ENCODE data that show differential enrichment at the level of the user's data but turn out to be technical replicates of the Inputs. The software determines the success of the experiment based on the q-values, and also reports some descriptive quality statistics, such as the percentage increase in mean tag density in IP compared to Input and the percentage of the genome classified as signal region. Because the q-values are computed based on human data, users should be aware that the q-values may not be relevant if their data are generated from other organisms.
CHANCE also provides a graphical visualization of IP strength with genome coverage, by plotting the empirical cumulative percentage of tags covered by the bins that are sorted in an increasing order of read density for both the IP and the Input. By examining and comparing the IP and Input curves, one may identify quality issues, such as insufficient sequencing depth, amplification bias, and weak IP enrichment.
between treatment samples without controls (i.e., using one of them as a control) is not recommended because highly enriched regions could be identified due to artefacts or different chromatin structure and not due to true binding events.
Typically, both methodologies assume that significant (see section ''Peak Calling'') and reproducible (see section ''Assessment of Reproducibility'') peaks have been found in advance independently for each condition. In order to increase sensitivity for detecting differentially bound regions (at the expense of increasing the number of false positives), more relaxed thresholds can be used to find peaks at each condition. Then, depending on the biological question, the sets of peaks called in any of the conditions can be considered separately, or collapsed into one or more meaningful lists of consensus peak regions. One can use the qualitative approach to get an initial overview of differential binding. However, peaks identified in all conditions will never be declared as differentially bound sites by this approach based just on the positions of the peaks [33]. The quantitative approach works with read counts (e.g., DBChIP [33]) or read densities (e.g., MAnorm [34]) computed over peak regions, and has higher computational cost, but is recommended as it provides precise statistical assessment of differential binding across conditions (e.g., p-values or q-values linked to read-enrichment fold changes). It is strongly advised to verify that the data fulfill the requirements of the software chosen for the analysis. For instance, DIME [35] assumes that a significant proportion of peaks are common to the conditions under comparison, MAnorm assumes that peaks that are common in both conditions do not change significantly, while other methodologies may expect a constant number of peaks across conditions [25]. Importantly, with some tools only two conditions can be submitted simultaneously for comparison (e.g., MAnorm), and some may perform better depending on the protein ChIPed (e.g., ChIPDiff [36] for histone marks and POLYPHEMUS [37] for RNA Pol II).

Peak Annotation
The aim of the annotation is to associate the ChIP-seq peaks with functionally relevant genomic regions, such as gene promoters, transcription start sites, intergenic regions, etc. In the first step, one uploads the peaks and reads (in an appropriate format, e.g., BED or GFF for peaks, WIG or bedGraph for normalized read coverage; see Text S1 and [38][39][40][41]) to a genome browser, where regions can be manually examined in search for associations with annotated genomic features. If comparable data (e.g., ChIP-qPCR) is available, it can be compared with the ChIP-seq peaks and reads manually in the browser as well. A systematic analysis can also be performed using tools in packages such as BEDTools to compute the distance from each peak to the nearest landmark (e.g., TSS), or to identify the genes within a given distance of a peak. The output of such ''location analyses,'' obtained for instance using CEAS [42] or the Bioconductor package ChIPpeakAnno [43], can be further correlated with expression data (e.g., to determine if proximity of a gene to a peak is correlated with its expression) or subjected to a gene ontology analysis (e.g., to Box 3. Peak calling: Punctate-source transcription factors Nowadays, since ChIP-seq data of point-source factors are the most abundant type, most peak callers are designed and fine-tuned for these factors. Existing peak callers differ from each other in terms of signal smoothing and background modeling. Because DNA around interaction sites is more easily sheared, the ends of ChIPed DNA fragments would form footprints on DNA whose size is more related to protein-DNA interaction than to size selection during library preparation. Those peak callers able to capture this experiment-specific information can greatly improve accuracy of prediction. For example, peak callers SPP [9] and MACS [17] (version 2) use cross-correlation to find the lag between reads mapped to the minus and the plus strand as the size of actual protein-DNA interacting regions. After smoothing, background models are then used to remove noise either directly from the control sample or from features of the genome sequence such as GC content or mappability (BEADS [84]). Peaks are finally called above a user-defined SNR level. Models used for the statistical assessment of enriched regions (peaks) range from Poisson (CSAR [85]), local Poisson (MACS), negative binomial (CisGenome [56]) to zero-inflated negative binomial (ZINBA [86]), or even extend to more sophisticated machine learning modeling techniques such as Hidden Markov Model (HPeak [87] and BayesPeak [88]).
Most peak-calling algorithms apply a window-based method to detect peaks, so nearby binding events may be erroneously merged. To improve the spatial resolution of binding event predictions, several peak callers use peak shape as a clue. PeakSplitter [89] can look for local maxima in a broader region containing several sub-peaks. GPS [67] builds a probabilistic model of the distribution of ChIP-seq reads at given peak candidate regions to deconvolve nearby homotypic events. The R packages polyaPeak and NarrowPeaks can analyze the shape of the peaks to re-rank and narrow down the final peak list, respectively. These approaches are highly recommended as a post-processing step after general peak calling for point-source factors.

Box 4. Peak calling: Broad enriched regions from histone marks
Due to the increasing interest in epigenetic regulation, epigenetic marks such as histone modifications, DNA methylation, and chromatin remodeling factors are being explored through ChIP-seq. Some of these marks are enriched strongly in narrow genomic regions (e.g., H3K4me3 at gene promoters), and the peak callers appropriate for point-source factors (discussed in Box 3) can be used. However, most histone marks tend to have more broadly spreading and weaker patterns (e.g., H3K27me3). Several peak callers are specifically designed for predicting broad regions from ChIP-seq data, including SICER [90], CCAT [91], ZINBA, and RSEG [92]. Other peak callers including SPP, MACS (version 2), and PeakRanger [93] can also be used with this type of ChIP-seq data by using their options to increase ''bandwidth'' or to relax the ''peak cutoff.'' For broad marks, the pattern of enrichment should be described as ''domains'' instead of ''peaks'' because there are no clearly defined peak summits. An alternative representation of the pattern of mapped reads is hierarchical: combining multiple levels of enrichment. For example, MACS (version 2) and Scripture [94] (originally designed for RNA-seq) can make narrow calls for strong enrichment inside broader calls for weak enrichment associated with domain boundaries.
determine if the ChIPed protein is involved in particular biological processes). Gene ontology analysis can be done using DAVID [44], GREAT [45], or GSEA [46]. Sometimes, the reads densities relative to a specific annotated feature are plotted and compared across different samples, thus revealing protein-binding pattern differences between them [47].

Motif Analysis
Motif analysis is useful for much more than just identifying the causal DNAbinding motif in TF ChIP-seq peaks. When the motif of the ChIPed protein is already known, motif analysis provides validation of the success of the experiment.
Even when the motif is not known beforehand, identifying a centrally located motif in a large fraction of the peaks by motif analysis is indicative of a successful experiment. Motif analysis can also identify the DNA-binding motifs of other proteins that bind in complex or in conjunction with the ChIPed protein, illuminating the mechanisms of transcriptional regulation. Motif analysis is also useful with histone modification ChIP-seq because it can discover unanticipated sequence signals associated with such marks. Table S4 and [48,49] list a small sample of the publicly available tools for motif analysis.
Motif analysis is applied to the genomic regions identified by peak-calling algorithms.
Hence, the first step in motif analysis is to assemble a set of genomic sequences in FASTA format corresponding to all the significant ChIP-seq peaks [50][51][52][53][54]. The second step in motif analysis is motif discovery and it is advisable to input the peak sequences to two or more of the many algorithms able to discover sequence motifs in unaligned DNA sequences [55][56][57][58], as the algorithms have complementary strengths and weaknesses. Some motif discovery algorithms form part of pipelines that perform several motif analysis steps (e.g., MEME-ChIP [57] and peak-motifs [58]), including word-based motif discovery algorithms and motif enrichment algorithms that can identify motifs present in only a small fraction of the peaks. Following motif discovery, comparing the discovered motifs with known DNA motifs using motif comparison software [59,60] is useful to confirm the presence of the ChIPed TF motif if its (or its TF-family) binding motif is known. The results will also provide hints about other TFs that bind near the ChIPed TF. Next, central motif enrichment analysis will determine if other known DNA motifs are enriched near the centers (or summits) of the ChIP-seq peaks [61]. It can also be useful to perform local motif enrichment analysis on regions centered on genomic landmarks such as transcription start sites overlapped by ChIP-seq peaks [61]. Additionally, a motif spacing analysis detects preferred distances and arrangements of pairs of motifs that can be indicative of physical interactions between TFs [62]. Finally, motif prediction maps and visualizes the genomic locations of the motifs in each of the ChIP-seq regions [63,64]. In this step, the discovered or enriched motifs are used to scan the ChIPseq peak regions, and the coordinates of the matches are uploaded to a genome browser for visualization.

Outlook
The challenges of ChIP-seq require novel experimental, statistical, and computational solutions. Ongoing advances will allow ChIP-seq to analyze samples containing far fewer cells, greatly expanding its applicability in areas such as embryology and development where large samples are prohibitively expensive or difficult to obtain. Nano-ChIP-seq can analyze a sample as small as 10,000 cells [65]. No less critical is to trim today's peaks that are much wider than the actual transcription factor binding sites. This is necessary to distinguish artefacts from bona fide joint binding events: most transcription factors competitively, cooperatively, or co-bind with other Box 5. Peak calling: Mixed signals There are also some factors (such as RNA Pol II) that bind to DNA in regions with bigger variation. It is known that some RNA Pol II complexes are stalled while others are moving along with active transcription [95]. In the first case, data ideally should be treated as for a point-source factor, whereas in the second case, the data should be treated as for factors with broad marks. An ideal algorithm should accommodate both patterns, which means peak calling should be more general. Some tools have options for both narrow and broad peak calling, such as SPP, MACS, ZINBA, and PeakRanger. However, with careful parameter tweaking any algorithm suitable for broad peak detection would work for this type of data.

Box 6. Normalization
Whether comparing one ChIP sample against input DNA (sonicated DNA), ''mock'' ChIP (non-specific antibody, e.g., IgG) in peak calling, or comparing a ChIP sample against another in differential analysis, there are linear and non-linear normalization methods available to make the two samples ''comparable'' (Table  S2). Although many methodologies focus on normalization to a control sample, none of them make the distinction on the type of control samples used. An intuitive and commonly used linear normalization technique is called sequencing depth normalization. In this method the number of reads is multiplied by a scale factor to make the total reads in different samples the same (see [9,96] for details). A slight modification of the method is used in PeakSeq [24], where a scale factor is estimated in a region (,10 Kb) using linear regression. Many other existing methods also use a normalization factor to linearly scale samples, focusing on normalization against control samples (see for example CisGenome [56], MACS [17], and USeq [97]). Another scaling normalization method known as RPKM (Reads per Kilobase of sequence range per Million mapped reads) proposed in [98] adjusts for biases due to the higher probability of reads falling into longer regions.
A non-linear normalization adjusts for biases with non-linear trend. In a method described in [28] the data is normalized with respect to mean and variance using locally weighted regression (LOESS). It is based on the assumption that the effect of biological condition change does not cause global binding alterations. This assumption can be applied, for example, when comparing samples with different stages of disease progression, or on samples before and after a certain treatment (see section ''Differential Binding Analysis''). A modified version of this non-linear normalization is implemented as MAnorm [34], assuming that peaks common in the two conditions do not undergo global changes. The R package called POLYPHEMUS [37] has also been developed, implementing two normalization methods: (1) the non-linear method described in [28] and (2) a Quantile normalization that makes the distribution in different samples the same. Normalization issues are, at present, not fully exploited although they might have a substantial impact on the results [28,37,99].
transcription factors, the transcriptional machinery, or cofactors. The effects of context-dependent regulatory mechanisms can fundamentally differ from the effects of individual binding events [66]. To address this issue, the Genome Positioning System (GPS) resolves closely spaced peaks using a segmented expectation maximization algorithm [67]. A promising experimental method for localizing narrow peaks is ChIP-exo that uses bacteriophage l exonuclease to digest the ends of DNA fragments not bound to protein [68].
The number of false positive peaks can be reduced both experimentally and computationally. Improving antibody specificity is a long-term endeavor, and despite impressive progress, still a quarter of histone modification antibodies fail the specificity test [69]. Another way to eliminate massive amounts of false positive peaks is to limit the regulatory binding sites to nucleosome-depleted regions, which are accessible for regulator binding. These regions are mapped by DNase I hypersensitivity sequencing (DNase-seq) and similar techniques: Thurman et al. found that 94% of the human transcription factor binding sites fell into DNase hypersensitivity regions with only a few exceptions like the transcription factors ZNF274, KAP1, and SETDB1, which also bind to closed chromatin [70]. False positive peaks are also due to unrealistic p-values (and hence FDRs) coming from unrealistic statistical models used in most methods [71]. The computational analysis of peak calling is still in its infancy, expanding the diverse and condition-specific performance of the methods [72,73], therefore we recommend using several methods for peak calling.
Perhaps the most important novel developments are related to the detection and analyses of distal regulatory regions, which are distant in sequence but brought close in 3-D space by DNA bending. To reveal such 3-D mechanisms of transcriptional regulation, two major techniques have emerged: chromatin interaction analysis by paired-end tags (CHIA-PET) [74] and chromosome conformation capture assays such as circular chromosome conformation capture (4C) [75] or chromosome conformation capture carbon copy (5C) [76].
Biological functions of binding sites are not necessarily indicated by the reproducibility of peaks or FDR/IDR values (Box 8, [7,23,77,78]). This issue re-emerged during the ENCODE Project that produced unprecedented regulatory information [66,79] under rigorous quality standards [7]. DNA-protein binding is dynamic, and the measured strength of a binding event depends (among other things) on the fraction of cells in the (often inhomogeneous) sample where it occurs, as well as the proportion of the time it is occupied in a given cell. Hence, ''weak'' binding sites, regardless of what significance threshold is used, may have strong biological functions [80][81][82]. ChIP-seq will also detect indirect DNA binding by the protein (via another protein or complex), so predicted sites not containing the motif may also be functional. Finally, binding does not necessarily imply function, so it will remain necessary to use additional information (such as expression or chromatin conformation data) to reliably infer the function of individual binding events [83].
The diverse experimental and computational methods discussed here are revolutionizing our understanding of the complex networks that, by regulating transcription, impact translation and almost all biological processes.

Supporting Information
Text S1 Standard graphing track data formats for genome browser visualization. (DOCX) Figure S1 Assessment of read quality using strand cross-correlation. Strand cross-correlation is computed as the Pearson correlation between the positive and the negative strand profiles at different strand shift distances, k. The cross-correlation (panel A) usually peaks at two distances of shift, one corresponding to the read length, and one to the average fragment length of the library. The absolute and relative height of the two peaks is useful for assessing IP enrichment. Adapted from Landt et al. [7]. (TIF) Figure S2 The irreproducible discovery rate (IDR) framework for Box 7. Duplicated reads Duplicate (identical) reads present a challenge because they can arise from independent DNA fragments or by PCR amplification of a single fragment. In the former case, the duplicate reads are signals, in the latter case they are noise (experimental artefact). A safe solution is to keep a fixed number of hits per genomic location (considering different strands as different locations) according to sequencing depth, and in this way better specificity (fewer false positive peaks) can be achieved [8]. However in terms of estimating the protein's affinity for a given genomic region, it is more reasonable to consider all hits. This can be done in a well-designed pipeline with certain steps before and after peak calling. For example, one can remove a certain number of duplicates to call confident peaks, and then put duplicates back to refine properties of these peaks such as peak height and boundaries.

Box 8. Irreproducible discovery rate (IDR)
Given a set of peak calls for a pair of replicate data sets, the peaks can be ranked based on a criterion of significance, such as the p-value, the q-value, or the fold enrichment. Significant peaks generally are ranked more consistently across the replicates than the peaks with low significance. This provides an indicator of the transition from real signal to noise. IDR [23] quantifies this transition by classifying peaks into a reproducible and an irreproducible group, where the peaks in the reproducible group should be ranked higher and more consistently across replicates than the irreproducible group ( Figure S2). It assigns each signal a reproducibility index, which estimates its probability to be reproducible, and also reports the expected rate of irreproducible discoveries in the selected peaks (referred to as IDR) in a fashion analogous to that of false discovery rate (FDR). An R package for computing IDR is given in [23]. Prototypical examples illustrated using ENCODE data may be found in [7]. When using IDR, a relatively relaxed peak-calling threshold is advised because the IDR algorithm requires sampling of both signal and noise distributions to assess the reproducibility of peaks.
A major advantage of the IDR method is that it is independent of the peak-calling algorithms and can be applied to a variety of significance criteria, across labs and platforms. It has been shown that it produces a stable threshold that is more consistent across laboratories, antibodies, and analysis protocols (e.g., peak callers) than FDR measures [7]. assessing reproducibility of ChIP-seq data sets. Panel A shows a scatterplot of the significance scores of peaks identified in two replicate ChIP-seq experiments. The IDR method classifies peaks into reproducible (black) and irreproducible (red) groups, and computes for each peak the probability that the peak belongs to the irreproducible group. It ranks and selects peaks according to this probability, and computes IDR, the expected rate of irreproducible discoveries in the selected peaks. Panel B shows the estimated IDR at different rank thresholds when the peaks are sorted by the original significance score. (TIF) Table S1 Examples of peak callers employed in ChIP-seq. The list includes tools that allow the processing and post-processing of diverse types of narrow read-enriched regions (peaks), broad enriched regions (domains), and mixed signals such as in RNA Pol II ChIP-seq. (DOCX)   . The table gives examples of  publicly available software tools for per-forming motif analysis on ChIP-seq peaks or nearby genes. The tools are grouped by the major task (''category''), and checkmarks indicate the specific steps that each tool performs. Web-based motif discovery input size limits-ChIPMunk: unknown; CompleteMOTIFS: 500,000 base pairs; MEME-ChIP: 50,000,000 base pairs; peak-motifs: no limit; Cistrome: 5,000 peaks. (DOCX)


extfieldsofstudy: Biology Medicine

text:
Transcriptional Response of Bovine Monocyte-Derived Macrophages after the Infection with Different Argentinean Mycobacterium bovis Isolates

Infection of bovines with Mycobacterium bovis causes important financial hardship in many countries presenting also a risk for humans. M. bovis is known to be adapted to survive and thrive within the intramacrophage environment. In spite of its relevance, at present the information about macrophage expression patterns is scarce, particularly regarding the bovine host. In this study, transcriptomic analysis was used to detect genes differentially expressed in macrophages derived from peripheral blood mononuclear cells at early stages of infection with two Argentinean strains of M. bovis, a virulent and an attenuated strains. The results showed that the number of differentially expressed genes in the cells infected with the virulent strain (5) was significantly lower than those in the cells infected with the attenuated strain (172). Several genes were more strongly expressed in infected macrophages. Among them, we detected encoding transcription factors, anthrax toxin receptor, cell division and apoptosis regulator, ankyrin proteins, cytoskeleton proteins, protein of cell differentiation, and regulators of endocytic traffic of membrane. Quantitative real-time PCR of a selected group of differentially expressed genes confirmed the microarrays results. Altogether, the present results contribute to understanding the mechanisms involved in the early interaction of M. bovis with the bovine macrophage.

Introduction
Mycobacterium bovis is an important aerobic pathogenic bacterium and the causative agent of most cases of tuberculosis in bovines (BTB). M. bovis also poses a zoonotic risk to herd keepers and communities that have close interactions with infected livestock [1]. e impacts of BTB infection are manifold: not only the losses due to the depopulation of herds and the social problems of removing cattle from small holdings, but also the vast economic consequences, particularly in the developing countries, that result from herd restrictions, disruptions to trade, and reduced agricultural productivity [2]. e pathogenesis of BTB results from a complex interaction between the pathogen and the host [3]. e gateway of the mycobacteria into the host is through macrophages, where the bacteria can survive a long period of time without the appearance of any symptoms [4]. In response to factors that are still unknown, the population of bacteria begins replication and escapes the infected macrophages causing illness in a small proportion of individuals. e macrophage plays a key role as an effector cell, activating innate immune responses and initiating acquired immune response. Innate immunity provides an important early defense against M. bovis and is critical in determining the clinical outcome. In order to survive the immune response, mycobacteria have developed a variety of mechanisms that prevent them from being killed by their host cells [4]. In addition, the level of pathogenicity of mycobacterial species in different hosts can be related to their effects on the host's immune response. It has been observed that the efficiency of phagolysosomal fusion and the rate of bacteria proliferation in macrophages vary among different M. tuberculosis clinical isolates [5]. Furthermore, the initial interaction between macrophages and mycobacteria is thought to play a key role in determining the outcome of infection [4].
Large-scale transcriptional gene expression analysis has been used to de�ne the repertoire of genes expressed in hostpathogen interactions of tuberculosis disease and many other infectious diseases [4]. ese transcriptional approaches have signi�cantly contributed to understanding the mechanisms involved in these interactions. Interactions of pathogenic Mycobacterium species with macrophages from different hosts have been assessed using gene expression pro�les [6][7][8][9][10]; however, there is small number of studies where the early responses of bovine macrophages to M. bovis have been previously evaluated [11].
For this reason, the aim of the present study was to compare the early gene expression pro�le of bovine monocytederived macrophages (MDMs) obtained from peripheral blood mononuclear cells (PBMCs) following infection with either an attenuated or a hypervirulent M. bovis strain. ese distinct M. bovis strains were previously studied according to their virulence and replication rate in murine and bovine macrophages [12,13]. As a result, it has been found that 464 genes were signi�cantly differentially expressed during the early cellular events of macrophage infection. Analysis of the data sheds light on the crosstalk between intracellular pathogens and their host cells, highlighting possible novel mechanisms of bacterial evasion from immunological elimination.

Bacterial Strains and Culture
Media. e M. bovis strains used in this work were isolated from an Argentinean wild boar (04-303) and from an Argentinean bovine (04-534), respectively. In a previous study by Aguilar and colleagues [12] based on a mouse model of progressive pulmonary tuberculosis, M. bovis strain 04-534 showed an attenuated phenotype, whereas strain 04-303 was the most virulent in mice. erefore, these strains were selected since they represent the extremes in terms of virulence among the ones studied. In addition, the replication rate of strain 04-303 in bovine macrophages has been shown to be signi�cantly higher than that of the strain 04-534, as was described by Blanco and colleagues [13].
Strains were grown in either Middlebrook 7H9 medium or Middlebrook 7H11 medium, both supplemented with albumin/dextrose/catalase, 0.4% pyruvate, and 0.05% Tween 80 to avoid the formation of clumps. Culturing of all M. bovis strains was performed in a biosafety containment level (CL3) laboratory.

PBMCs Isolation and MDMs
Infection. e animals used in this study were selected from the experimental herd of INTA and were tested negative for bovine tuberculosis infection for the single intradermal tuberculin test. Six hundred microliters of blood were taken in sterile conditions according the instructions of the Committee for Institutional Care and Use of Animal Experimentation (CICUAE) of INTA. PBMCs were separated from heparinized blood by centrifugation over histopaque 1077 (Sigma, USA) following manufacturer's protocol. To derive monocytes, PBMCs were seeded in T75 �asks (2 × 10 8 ) containing RPMI1640 complete medium (Invitrogen, USA) supplemented with 10% of autologous plasma for 16 h at 37 ∘ C, 5% CO 2 . Nonadherent cells were removed by washing with Phosphatebuffered saline (PBS). Only adherent cells were maintained in culture for 5 days at 37 ∘ C, 5% CO 2 to obtain MDMs [14]. Cell viability was con�rmed by trypan blue exclusion assay [15]. e M. bovis cultures were harvested, washed to eliminate rest of bacteria culture medium, suspended in RPMI medium, vortexed, sonicated for 1 min in an ultrasonic cleaner, and passed through a syringe needle (25 gauge) to disaggregate bacteria clumps. Viable bacteria were counted with Live/Dead BacLight Bacterial Viability Kit (Molecular probes, Invitrogen) and then used to infect the cultured bovine macrophages. e infection was performed at a multiplicity of infection (moi) of 5 bacteria/cell. Infected cells were incubated at 37 ∘ C, 5% CO 2 for 4 h and then washed three times with fresh RPMI 1640 medium to eliminate extracellular bacteria. Two independent sets of RNA samples were prepared; one was for the microarrays and the other for the reverse transcription quantitative real-time polymerase chain reaction (RT-qPCR) analysis. ree infections were performed for each set.

Isolation of RNA from Macrophages.
Four hours aer infection, the cells were scraped and lysed with 1 mL of chilled Trizol (Invitrogen, Carlsbad, USA). Trizol was removed from the �asks and the lysate homogenized. Cellular RNAs extractions were then performed according to the manufacter's instructions for Trizol reagent (Invitrogen, USA). RNA pellets were suspended in 30 L RNase-free water. RNA samples were further puri�ed using RNeasy Mini Elute Cleanup Kit (Qiagen, Melbourne, Australia) eluted from the column in 20 L RNase-free water. e RNA quality and integrity was checked before hybridizations using a Bioanalyser 2100 (Agilent).

Microarrays
Hybridization. Ten microarray slides were processed in total, representing the biological replicates per group: three for the cell infected with the M. bovis virulent strain (04-303), three for cells infected with the M. bovis attenuated strain (04-534), and four for the noninfected cells. Affymetrix GeneChip Bovine Genome Array platform (Affymetrix, Santa Clara, USA) was used in this study. e array contains 24,027 probe sets representing over 23,000 transcripts from Bos taurus and include approximately 19,000 annotated UniGene clusters. e experiment was designed to be compliant with Minimum Information About a Microarray Experiment (MIAME) standards. Each RNA sample was processed and hybridized to individual slides. Target preparation including veri�cation of RNA quality assessed using a Bioanalyser 2100 (Agilent Technologies, Santa Clara, USA) and microarray processing procedures were carried out at the Affymetrix facility at the School of Agronomy of the University of Buenos Aires, as described in the Affymetrix GeneChip Expression Analysis Manual (Affymetrix), and scanning was performed with a Microarray Scanner 3000 7G (Affymetrix).

Statistical Analysis of Microarrays
Data. e analysis of expression data was performed with different packages from the Bioconductor project (http://www.bioconductor.org/) [15], an extension for bioinformatics of the R statistical language (http://www.r-project.org/) [16]. e data quality of ten CEL �les obtained from the Affymetrix chips were assessed with packages affyPLM [17] and simple affy. On the basis of the boxplots of normalized unscaled standard errors, relative log expressions, and the clustering of between array distances (measured by median absolute deviations), one chip of each treatment was removed from the analysis. e remaining good-quality data was normalized using function rma of the package Affy [18], which essentially reproduces the standardization procedure of the Affymetrix MicroArray Suite (MAS) soware.
To �lter the arrays, the gene�lter package was used: �rst those probesets with low variance across samples were removed, and then probesets with little variation between treatments or high levels of noise were �ltered out using an Anova �lter with a relaxed cutoff of . e detection of probesets with differential expression between macrophages infected with virulent or attenuated strains and the uninfected macrophage control was performed with linear models and empirical Bayes methods to correct for multiple comparisons as implemented in the Linear Models for Microarray Data (LIMMA) package [19]. e false discovery rate was set at . All microarray data were MIAME compliant and have been submitted to the NCBI Gene Expression Omnibus (GEO) database with experiment series accession number GSE39819.

Validation of Microarrays
Results by RT-qPCR. e validation procedure was performed using PBMCs collected (three times) from different animals from those sampled for microarrays experiments. e cells were infected as mentioned above.
DNA-free RNA (1 g) was mixed with 50 ng of random primers (Invitrogen) in 20 L of �nal volume and reverse transcribed to total cDNA with SuperScript II reverse transcriptase (Invitrogen, USA) following the manufacturer's instructions. One microliter of the cDNA was used as template for each real time quantitative PCR (qPCR) reaction.
All primers were designed using Primer Express Soware to span an intron-exon boundary and to anneal only to cDNA synthesized from spliced mRNAs. Primer sequences are listed in Table 1. qPCR reactions were performed as described by Blanco and colleagues [13]. Brie�y, S�BR green QuantiTec Mastermix (Qiagen) was employed, and reactions were made on Applied Biosystem 7000 SDS using standard cycling conditions. All reactions were performed in duplicate, and qPCR data were analyzed using the 2ddCT with efficiency correction as described previously [21] to assess differences on gene expression in macrophages within groups, RNA polymerase II (RPII) gene was used as the control gene, and data from control group was used as the calibrator [22]. Relative Expression Soware Tool (REST) beta 9 soware (http://www.gene-quanti�cation.de/rest-2009.html) was used for �nal calculations and statistical analysis [21].

DNA Microarray Analysis.
Bovine monocyte-derived macrophages were infected with either a virulent (Argentinean wild boar 04-303) or an attenuated M. bovis strain (Argentinean bovine 04-534) for 4 hours, and cellular transcriptomes were analyzed and compared to the reference transcriptome of uninfected cells.
Comprehensive gene expression pro�les for the three attenuated and the three virulent M. bovis macrophage infections were generated with high-density oligonucleotide bovine arrays (Affymetrix) and 24,072 probe sets, which in total interrogated the expression levels of approximately 23,000 transcripts. Aer the 2-step �ltering, 464 probesets were retained. In order to assess the statistical sig-ni�cance of the expression pro�les for the 464 probesets obtained in the previous step, LIMMA's linear models and empirical Bayes corrections for multiple comparisons were applied. e number of probe sets with a statistically sig-ni�cant fold change was highly in�uenced by the �ltering and the statistical testing settings. Considering a threshold of 0.05 for the adjusted P value, 310 unique probesets showed a differential expression that was statistically sig-ni�cant (see Supplementary Material S1 available online at http://dx.doi.org/10.1155/2013/458278). e results showed that the number of upregulated genes (161) when MDMs were infected with both strains versus control was signi�cantly higher than the number of downregulated genes (15). Furthermore, the number of differentially expressed genes in the cells infected with the virulent strain (5) was signi�cantly less than those in the cells infected with the attenuated strain (172). In Table 2, these results are summarized.

Description of Highly Differentially Expressed Genes.
Among the �ve exclusively differential expressed genes in the cells infected with the virulent strain (04-303), three were upregulated genes (Table 3 and Supplementary Material S1). One is the interferon-gamma-inducing factor, IL18. IL-18 together with IL-12 cytokines is primarily produced by dendritic cells and macrophages in response to Toll-like receptor (TLR) signalling interaction with tubercle bacilli [23]. Another upregulated gene is the ankyrin repeat domain 17, ANKRD17. e ankyrin superfamily is composed of proteins that are ubiquitously expressed and typically found within the membrane-associated cytoskeleton. Ankyrins contain repeat domains that interact with other ankyrin repeats. ese ankyrin interactions act in many diverse functions T 1: Primers sequences used for RT-qPCR.

Genes
Forward ACCAAGTTCCACTCGCTGAA GCTCATCTCCTCCTGGCTAA  MPTN  TTGAAGCCACTGACAACCAG  AACAGACAGACAGGCAGCAG  RQCD1  GCCAAAGGACAGCAGGAAC  CAGCGAAAGCCATAAGCAA  ANTXR2  CCTCCTCCACCACCACCT  TTATCACCCCAACGAACCTC  XIST  GGGGTGTTTGACCGTTACAT  TCCCTCCTTTCACTTTGTCC  RPII GCACCACGTCCAATGACAT GTGCGGCTGCTTCCATAA T 2: Analysis of probesets with differential expression. and include regulation of transcription, cell cycle, cell fate determination, cytoskeletal integrity, cellular mechanosensation, and endocytosis [24,25]. TCDD-inducible poly (ADPribose) polymerase, TiPARP, is the last upregulated gene in this group. Poly(ADP-ribose) polymerases catalyse the covalent attachment of ADP-ribose units from NAD+ to itself and to a limited number of other DNA binding proteins, which decreases their affinity for DNA. Poly(ADPribose) polymerase is a regulatory component induced by DNA damage. 2,3,7,8-Tetrachlorodibenzo-p-dioxin (TCDD) causes pleiotropic effects in mammalian species through modulating gene expression. e gene TiPARP is a target of TCDD that may contribute to multiple responses to TCDD by modulating protein function through poly ADP-ribosylation [26].
Genes exclusively upregulated in the cells infected with the attenuated strain (04-534) are involved in different cellular processes (Table 3). DNA transcription and RNA metabolism are some one of these processes and some representative genes that include DNA binding genes, and RNA metabolism genes or transcription factors are present in this group: TAR DNA binding protein (TARDBP); RNA binding motif protein 18 (RBM18); Rho-associated, coiledcoil containing protein kinase 2 (ROCK1); and the gene CREB/ATF bZIP transcription factor (CREBZF) which is a member of the CREB/ATF family of transcription factors that play an important role in cellular growth, metabolism, and survival [27].
e gene clathrin interactor 1 (CLINT1) was found also upregulated in this group. It was reported that CLINT1 might participate in the formation of clathrin-coated vesicles at the level of the trans-Golgi network (TGN). Clathrincoated vesicles participate in receptor-mediated endocytosis, in the transport of lysosomal enzymes from the TGN to the endosomes/lysosomes, and in the sorting of receptors within the endosomal system [28].
Another important overexpressed gene is the interferon alpha receptor 1 (IFNAR1). Although type I interferon is key to the immune response to viral pathogens, its role in bacterial infections is less understood. Recently mice lacking this receptor (IFNAR −/− ) have enhanced their resistance to infection with Listeria monocytogenes, an intracellular pathogen like Mycobacterium sp. [29]. Importantly, upregulation of type-1 interferon ( / ) mediated signalling has been recently described during active TB disease in patients and also in the lungs of M. tuberculosis infected mice and in infected human macrophages [30].
e group of genes that was signi�cantly upregulated (fold change >2) both in the macrophage infected with virulent and attenuated M. bovis strains, encode proteins involved in transcription process, in intracellular trafficking, and in the regulation of cell processes. Established functions have been previously assigned for several of these differentially expressed genes (Table 3).
Anthrax toxin receptor 2 (ANTXR2) is one of the genes up-regulated in both macrophage infections. It has been demonstrated that ANTXR2 is expressed in primary mononuclear phagocytes and facilitates the entry of Bacillus anthracis exotoxins into the cytosol of susceptible cells [31].
Other genes upregulated in both macrophage infections were myotrophin (MPTN) the EH domain containing 4 (EHD4). MPTN is another member of the ankyrin superfamily previously mentioned. Members of the EHD protein family (EHD1-4) have emerged as critical regulators of endocytic membrane traffic as well as cell surface receptors [32]. Remarkably, an EH domain has been found upregulated in human alveolar macrophages infected with M. tuberculosis H37Rv [8].
Genes encoding regulator proteins or intracellular traf-�cking proteins were also found differentially expressed. Such are the examples of the nuclear receptor corepressor 1 (NCOR) which has a role in the regulation of in�ammatory response genes in macrophages [33] and the kinesin family member 3 (KIF3A) which regulates the vesicular transport in leucocytes and monocytes [34]. It has been reported that kinesins KIF5B and KIF5C are directly associated with RanBP2, a zinc �nger, and RAN binding domain containing protein [35]. Interestingly, in the present work it was found that together with KIF3A, a zinc �nger RAN binding domain containing 2 (ZRANB2) was also upregulated in macrophages infected with both M. bovis strains. Further studies are necessary to demonstrate a functional connection between both proteins in the infected macrophages.

BioMed Research International
In addition, it has been reported that CARP/CCAR1 induces apoptosis in a manner dependent on phosphorylation of a speci�c tyrosine residue of CARP1 and downstream activation of p38 MAPK and caspase-3 and -9 in HBC cells [36]. Although this function of CARP1 has not been yet reported in monocytes, it is consistent with the apoptotic activity induced in the macrophage upon infection with some M. tuberculosis isolates.
Another overexpressed gene is ZBTB1 that encodes zinc �nger and BTB domains. Established functions for zinc �nger-BTB domains containing proteins have been reported, such as the promyelocytic leukemia zinc �nger protein (PLZF) that induces monocytic differentiation in hematopoietic cells [37].
Other proteins with a probable role in transcription, splicing and different aspects of cellular RNA metabolism whose genes were overexpressed aer macrophage infections are LOC515193, similar to SAFB like transcription modulator; a YTH protein family member; RQCD1, a protein required for cell differentiation, and the X-inactive speci�c transcript (XIST gene) which is the only gene known to be speci�cally transcribed from the inactive X chromosome in female somatic cells and is necessary and sufficient for the initiation and spread of X inactivation. e role of this last gene in the early cellular events of macrophage infection is unclear.
Although all these differential expressed proteins may play roles in cellular processes induced aer mycobacterial infections, understanding the relationship of these protein functions with M. bovis infection needs further investigation.
Gene Set Enrichment Analysis. e probesets differentially expressed in macrophage cells infected with the virulent and the attenuated strains versus the control cells were submitted to the DAVID web application for their functional clustering. is procedure clusters genes according to their similarities in functional annotations (Supplementary Material S2).
For the virulent versus control contrast, only one cluster with 5 genes was found, and their annotations were related to nucleotide binding and more speci�cally to ATP binding. A detailed description of the genes included in this cluster is provided below.
For the attenuated versus control contrast, seven clusters were found with DAVID's functional gene clustering. e �rst cluster contained 5 probe sets whose common annotations were RNA recognition motif and nucleotide binding.
e second cluster contained 7 probe sets with proteolysis-related annotations. e third cluster contained 5 probe sets enriched in regulators of small GTPases and related annotations, including GTPase-activating proteins (GAP) and guanidine exchange factor (GEF). GAPs are a family of regulatory proteins whose members can bind and activate small GTPases and others G proteins. Regulation of small GTPases is important because they participate in crucial cellular processes like cellular trafficking and cell cycling. A representative of this cluster is RASA2 that encodes a member of the GAP1 family. is GAP stimulates the GTPase activity of normal RAS p21. In this cluster, a gene containing a TBC domain and IQGAP1 is also present. IQGAP1 is a protein scaffold that binds to actin �laments and also to Rac and Cdc42, members of the Rho family of small GTPases [38], while TBC domains have been expected to function as a certain Rab-GAP [39]. On the other hand, the function of GEF proteins is to activate Rho GTPases, which are important regulators of multiple cellular activities and, most notably, reorganization of the actin cytoskeleton [40]. FGDA, a gene also present in cluster three, encodes a GEF protein, which activates Cdc42 that regulates a variety of cellular functions including cell shape, migration, endocytosis, and cell cycle.
e fourth cluster contained 17 probesets, and 3 of them were also part of the only cluster found for the virulent strain versus control. In addition to the nucleotide binding and ATP binding, there were other enriched annotations, like protein kinase and phosphorylation. Genes encoding members of mitogen-activated protein kinases (MAPKs) and cyclin-dependent kinases (CDKs) are present in this cluster. Kinases are known to regulate the majority of cellular pathways, especially those involved in signal transduction. In particular, MAP kinases are mainly involved in mechanisms of immune responses, from innate immunity to activation of adaptive immunity and to cell death, while CDKs (PRPF4B and CDK2) form a family of protein kinases mostly involved in regulating the cell cycle. Other groups of kinases included in the cluster are ToR (TOR1B) that participates in signaling pathways of autophagy [41] and tribbles (TRIB3) that regulate MAPKs [42].
e �h cluster contained 4 probe sets, all with annotations "purine nucleoside triphosphate biosynthetic process" and similar. ree of these clustered genes encode ATPases, one V-ATPase and two P-type cation transport ATPases that are involved in acidi�cation of intracellular compartments and calcium homeostasis. e V-ATPase is a multisubunit enzyme that mediates acidi�cation of eukaryotic intracellular organelles, and the role of this enzyme in the antimicrobial mechanism against pathogenic mycobacteria has been established. e P-type ATPases remove bivalent calcium ions from eukaryotic cells against very large concentration gradients and play a critical role in intracellular calcium homeostasis [43]. Calcium homeostasis participates in apoptosis in different systems and regulates important cellular events triggered upon infection of macrophages with pathogenic mycobacteria [44] in turn; mycobacteria have developed mechanisms to evade these antimicrobial responses [45,46]. e sixth cluster had 4 probesets include transmembrane and transport involved domains.
Finally, the seventh cluster included 4 probesets with "zinc binding annotations" and zinc �ngers domains.
Validation of Microarray Analysis by RT-qPCR. e microarray data for the genes whose expression varied at least two fold between conditions were validated for a subset of �ve genes by RT-qPCR. For this purpose, PBMCs were collected (three times) from different animals to those sampled for the microarrays experiments. e RT-qPCR results corroborated the microarray data for all genes tested. As shown in Figure  1, the fold changes determined by RT-qPCR were oen greater than the fold changes for the same genes determined by microarray analysis. us, the concordance between the microarray results and the RT-qPCR results was complete, supporting the statistical approach that was adopted in this study.

Discussion
In the present study, the early (four hours) gene expression pro�le induced in bovine blood MDMs aer infection with two M. bovis strains with distinct level of virulence was compared to that of uninfected macrophages [12]. PBMCs represent an accessible tissue for the development of improved diagnostics, and previous studies have shown that immune responses generated in the peripheral blood of bovines infected with tuberculosis re�ect those at the site of disease [47]. Comparative gene expression pro�le of virulent straininfected macrophages and attenuated strain-infected macrophages did not identify signi�cant differences suggesting that during the early events of intracellular replication both strains have equivalent �tness in spite of their differences in virulence [12] and in replication rate both in mice and bovine macrophages [13]. In agreement with this, a microarray analysis study on bovine alveolar macrophages infected with a virulent and an attenuated isogenic M. bovis strains has detected small differences in gene expression between both strains and has only identi�ed IL-8 gene as signi�cantly overexpressed in macrophages infected with the virulent strain [8]. As well, other differential expressed genes identi�ed in this microarray analysis study were not con�rmed as signi�cantly different in RT-qPCR validations [8]. Unexpectedly, in the work presented here, the expression of IL-8 did not show a signi�cant difference among the condition analyzed, suggesting differences in both infection models studied: alveolar and monocyte. e present work shows that the expression of few genes was considerably changed aer M. bovis infection, irrespectively of the strain used to infect the macrophages. However, at least �ve genes were highly upregulated in macrophages infected with the virulent M. bovis strains, which gives strength to the results obtained here.
As this study was limited to the analysis of M. bovismacrophage interactions, the possibility that the differentially expressed genes here identi�ed were part of a shared macrophage activation program induced by intracellular bacteria cannot be discarded. In this regards, Nau and colleagues have performed a comparative examination of the transcriptional responses of human macrophages to a variety of bacteria, including M. tuberculosis [48]. e authors have found that macrophages responded to bacteria with a robust, shared pattern of gene expression, encoding receptors, signal transduction molecules, and transcription factors. Remarkably, the set of differentially expressed genes identi�ed by these authors did not include any of the highly upregulated genes identi�ed in the present study, suggesting that there are differences in response according to different species, with bovine macrophages responding differently than human macrophages.
A previous global transcriptional pro�le analysis between macrophages and dendritic cells in humans has revealed that the expression of more than two thousand genes was altered during 48-hour infection with M. tuberculosis. is altered gene expression pro�le included genes involved in sensing pathogens, intracellular signalling and trafficking, cell motility, and cytoskeleton remodelling, presenting profound differences between cell type in genes involved in oxidative stress, intracellular vesicle acidi�cation and trafficking [7]. In concordance with these previous �ndings, the present study found that differential expressed genes were clustered in functional annotations involved in different aspects of phagocytosis, such as cellular trafficking, endocytocis, cell migration, morphology, and cell cycling. Furthermore, it was found that the most signi�cantly upregulated genes (fold change >4) in macrophages infected with M. bovis include those involved in DNA interaction, nuclear receptor, vesicular transport, and apoptosis regulation.
Silver and colleagues have reported that cytokines IL-23, IL-6, and TNF-were signi�cantly upregulated in human alveolar macrophages infected during 24 hours with H37Rv, while IL-12 was not [9]. In the present study, overexpression of these cytokines was not observed, and the reason for this discrepancy could be the different origin of the macrophages used in both studies. In fact, in the study by Silver and colleagues, it is shown that IL-23 is not produced in blood monocytes [9]. Remarkably, it was found in the present work that among the upregulated genes in macrophages infected with the virulent M. bovis strains was that encoding IL-18. IL-18 together with IL-12 cytokines, is primarily produced by dendritic cells and macrophages in response to Toll-like receptor (TLR) signalling interaction with tubercle bacilli [21]. In addition, both M. bovis strains induced the upregulation of IL-13 receptor. It has been demonstrated that IL-13 can subvert 1-mediated immunity and promote inappropriate activation of macrophages, abrogating IFN-induced autophagy-mediated killing of intracellular mycobacteria in murine and human macrophages [49]. One possible reason for the lack of other cytokine expression in our studies is the short period of infection. It has been reported that the peak of cytokine productions in bovine macrophages is 24 hours aer M. bovis infections [8].
Altogether, the present results contribute to de�ne the transcriptional responses induced during the early events that occur in M. bovis macrophage interactions. However, it is important to take into consideration that some of the gene expression changes observed in this study may not be speci�c for M. bovis infection and may represent a shared transcriptional program induced by infectious diseases. us, further studies are needed to investigate this possibility.


text:
Evaluating GWAS-Identified SNPs for Age at Natural Menopause among Chinese Women

Background Age at natural menopause (ANM) is a complex trait with high heritability and is associated with several major hormonal-related diseases. Recently, several genome-wide association studies (GWAS), conducted exclusively among women of European ancestry, have discovered dozens of genetic loci influencing ANM. No study has been conducted to evaluate whether these findings can be generalized to Chinese women. Methodology/Principal Findings We evaluated the index single nucleotide polymorphisms (SNPs) in 19 GWAS-identified genetic susceptibility loci for ANM among 3,533 Chinese women who had natural menopause. We also investigated 3 additional SNPs which were in LD with the index SNP in European-ancestry but not in Asian-ancestry populations. Two genetic risk scores (GRS) were calculated to summarize SNPs across multiple loci one for all SNPs tested (GRSall), and one for SNPs which showed association in our study (GRSsel). All 22 SNPs showed the same association direction as previously reported. Eight SNPs were nominally statistically significant with P≤0.05: rs4246511 (RHBDL2), rs12461110 (NLRP11), rs2307449 (POLG), rs12611091 (BRSK1), rs1172822 (BRSK1), rs365132 (UIMC1), rs2720044 (ASH2L), and rs7246479 (TMEM150B). Especially, SNPs rs4246511, rs365132, rs1172822, and rs7246479 remained significant even after Bonferroni correction. Significant associations were observed for GRS. Women in the highest quartile began menopause 0.7 years (P = 3.24×10−9) and 0.9 years (P = 4.61×10−11) later than those in the lowest quartile for GRSsel and GRSall, respectively. Conclusions Among the 22 investigated SNPs, eight showed associations with ANM (P<0.05) in our Chinese population. Results from this study extend some recent GWAS findings to the Asian-ancestry population and may guide future efforts to identify genetic determination of menopause.

Introduction
The occurrence and timing of menarche and menopause play major roles in a woman's life and reproduction-related events. Menopause is the final reproductive event and unambiguously marks the ovarian aging and exhaustion of the primordial follicle pool in women [1]. Menopause at a younger age often predicts an earlier induction of subfertility, sterility, and transition to cycle irregularity and vice versa [2].
Generally, age at natural menopause (ANM) ranges between 40-60 years with a median of 49 to 52 years, depending on the population [3]. Population-based studies report that the beginning of transition from regular to irregular menstrual cycles is 6-7 years ahead of ANM. This probably reflects the end of female fertility, which is nearly 10 years ahead of ANM [4]. Menopause age variation influences the levels of serum estrogen, follicle stimulat-ing hormone, and progesterone, all which affect the well-being of women. Series studies report that the variation of menopause is associated with several major age-related diseases [5], such as cardiovascular disease [6,7], breast cancer [8], osteoporosis [9], and depression [10]. So, identifying factors which determine variations of ANM may provide insights into the pathogenesis of these diseases.
ANM is a complex trait, influenced by multiple environmental and genetic factors [11]. ANM has high heritability (70%). Genome-wide association studies (GWAS) have identified multiple genetic loci which influence ANM [12][13][14][15][16]. All of these GWAS were exclusively conducted in European-ancestry populations. Given the considerable differences in genetic architecture including allele frequencies, linkage disequilibrium (LD) structure and genetic diversity across ethnic groups, it is important to investigate whether GWAS-identified variants are associated with ANM in non-European populations. In this current study, we report a systematic evaluation of the GWAS-identified loci in association with ANM among Chinese.
This analysis included women from Shanghai, China, who were participants in multiple ongoing GWAS: the Shanghai Breast Cancer Genetics Study (SBCGS) [17], the Shanghai Endometrial Cancer Genetics Study (SECGS) [18], the Shanghai Diabetes Genetics Study (SDGS) [19], the Shanghai Colorectal Cancer Genetics Study [20], the Upper GI Cancer Genetics Study, and the Panscan Cancer Genetics Study [21]. Among these studies, 4,106 reported detailed information of menopause. We excluded 573 women including 390 individuals whose menopause were disease-related, 158 individuals who took hormone replacement treatment, and 25 individuals with other non-nature cause menopause. Thus, 3,533 post-menopausal women were finally included in this study.
All parent studies are population-based and similar study protocols were used across studies to collect blood or buccal cell samples and relevant exposure information. Interviews were conducted in-person, and structured epidemiological questionnaires were administered by trained interviewers. Age at natural menopause (ANM) is defined as the age when menstrual periods ceased permanently and naturally (in years). Anthropometric characteristics including birth date, education, occupation, and age at menarche were collected in addition to information regarding individual habits, history of hormone replacement therapy, prior use of oral contraceptives, and post-disease status. Written, informed consent was obtained from all participants prior to interview, and study protocols were approved by the institutional review boards of all institutions involved in the study.

SNP Selection and Genotyping
To date, 28 SNPs have been reported to be associated with ANM in four GWASs with P,5610 28 among European-ancestry women [12][13][14][15]. These SNPs represent 19 independent loci with low pairwise linkage disequilibrium (LD) (r 2 ,0.2) based on the HapMap Asian data (Table S1). We also searched flanking 500kb of each of the 28 index SNPs using the HapMap data. There were three common SNPs (MAF.0.01) which are in strong LD (r 2 $0.8) with the index SNP in European-ancestry, but not in Asian-ancestry populations (r 2 ,0.2). Three such SNPs, rs4667673, rs2517388, and rs12294104, were also included in this study (Table S1).
For SNPs which were not directly genotyped through the GWAS chip, imputed data were extracted. Using the Han Chinese (CHB) and Japanese (JPT) data from HapMap Phase 2 (release 24) and Phase 3 (release 2) as the references, we used MACH 1.0 software (http://www.sph.umich.edu/csg/abecasis/MaCH/) to conduct imputation within each GWAS. The imputation quality RSQR was more than 0.91 for all SNPs included in the present study, with the exceptions of rs7246479 (0.79) and rs2277339 (0.77). Dosage data (with imputation uncertainty taken into account) were investigated in relation to ANM.

Statistical Analysis
ANM was presented as the means 6 SD between different groups. Student's t-test and analysis of variance (ANOVA) were applied to evaluate the differences of ANM by non-genetic factors. Linear regression was used to derive beta coefficients and 95% CIs In LD with GWAS-identified SNP rs10183486 in CEU (r 2 = 0.863) but not in CHB (r 2 = 0.005). b In LD with GWAS-identified SNP rs2517388 in CEU (r 2 = 0.843) but not in CHB (r 2 = 0.169). for associations between genotypes and ANM under the additive model. All regression models include adjustment for birth cohort, education, and disease status. P-values reported for single SNPs are one-tailed since the hypothesis tested is whether the direction of association was consistent between Asian-ancestry and European-ancestry women. Program Quanto (http://hydra.usc.edu/ gxe/) was used to estimate statistical power.
The binomial sign test was conducted to assess the consistency of the associations previously reported and seen in our study. Genetic risk scores (GRS) were used to summarize multiple SNPs [23,24]. Two genetic risk scores were constructed in current study. The GRSsel included 8 SNPs which were associated with ANM (P,0.05) in our study; the GRSall included all 22 SNPs evaluated. Relationships between GRS quartiles and ANM were evaluated by linear regression and ANOVA. All statistical analyses were conducted using SAS, version 9.2 (SAS Institute, Inc., Cary, N.C.). Table 1 presents the mean ANM by education, occupation, number of live births, use of oral contraceptives, age at menarche, and disease status. Significant differences of ANM were observed across each of these categories with the exception of age at menarche. Table 2 shows the results of SNPs included in this study. 22 SNPs showed the same direction of association as previously reported (binomial sign test P,1610 24 ). Eight SNPs were nominally significant with P#0.05: rs4246511 (RHBDL2), rs365132 (UIMC1), rs2307449 (POLG), rs1246111 (NLRP11), rs12611091 (BRSK1), rs1172822 (BRSK1), rs2720044 (ASH2L), and rs7246479 (TMEM150B). Notably, SNPs rs4246511 (RHBDL2), rs365132 (UIMC1), rs1172822 (BRSK1), and rs7246479 (TMEM150B) showed associations with P#0.002 and remained significant even after Bonferroni correction. The perallele effect on ANM was 0.287, 0.284, 0.600, 0.487 years for rs4246511, rs365132, rs1172822, and rs7246479. These estimates are comparable with previous results in European-ancestry women. No significant associations were found for the other 14 SNPs. Among them, three have a MAF,0.05, including rs4886238, rs7123626, and rs7333181. Considering the disease of the subjects could also show overlap in the genetic architecture with menopausal age, the association analysis in the healthy controls was further performed as a sensitivity analysis and the results indicated that six SNPs of the above eight SNPs associated with ANM still showed association with P,0.05 (Table S2).

Results
GRS analyses also showed significant associations (Table 3). A dose-response association was observed between the number of effect alleles and ANM. The average ANM across increasing quartiles of the GRSsel were 49.21, 49.34, 49.60 and 49.90 (P = 3.24610 29 ). Similar association was observed for GRSall with women in the highest quartile starting menopause approximately 0.9 years later than those in the lowest quartile (P = 4.61610 211 ).

Discussion
In this study, we evaluated 22 SNPs in GWAS-identified genetic loci for menopause age originally discovered in Europeans in our Chinese population. Significant associations were observed for 8 SNPs and association directions were all same as originally reported. Furthermore, six of eight SNPs were replicated to be correlated with ANM in the healthy controls. Of the 8 replicated SNPs, four were located in 19q13.42, including rs12461110 at exon 5 of NLRP11gene, rs7246479 at exon 8 of TMEM150B gene, rs1172822 and rs12611091 at introns 19 and 5, respectively, of the BRSK1 gene. SNP rs12461110 results in a missense change from Proline to Leucine. The NLRP11 gene was reported to be involved in multiprotein complexes called inflammasomes and implicated in the activation of proinflammatory caspases [25]. SNP rs7246479 results in amino acid change from Leucine to Phenylalanine alanine. The TMEM150B gene encodes transmembrane protein 150B, which belongs to the DRAM/TMEM150 family and is involved in apoptosis pathway. SNP rs1172822 may affect the binding of the transcription factor CdxA and C/EBP binding sites (FASTSNP, http://fastsnp.ibms.sinica.edu.tw), whereas rs12611091 has no known or predicted function. BRSK1 is highly expressed in the human brain and moderately expressed in mammalian ovaries [26]. The BRSK1 was previously published by multiple studies in European ancestry populations and was recently replicated in Hispanic women [27]. Further functional exploration may shed insight to the role of the BRSK1 gene or other genes impacting menopause. SNP rs4246511 is located at intron 3 of the RHBDL2 gene, which may function as an intramembrane serine protease [28]. SNP rs2307449 is located at intron 18 of the POLG gene, which encodes the catalytic subunit of mitochondrial DNA polymerase. The enzyme is responsible for replication and repair of mitochondrial DNA, and modulates the balance between DNA synthesis and excision [29,30]. For the ASH2L gene, the previously reported SNP, rs2517388, did not replicate in our study. However, the other SNP rs2720044, which is in strong LD with rs2517388 in European-ancestry population with r 2 = 0.843 but not in Asianancestry with r 2 = 0.169, was significantly associated with ANM in our population. These results indicate that the ASH2L gene might harbor different susceptible variants to ANM across ethnic groups. The ASH2L is required during early embryogenesis [31] and recruited to the inactive X chromosome at the onset of stable X inactivation [32].
We did not replicate the association of other 14 SNPs with ANM. Similarly, in a recent study in Hispanic women, only two of 14 evaluated menopause loci were replicated [27]. For the 14 non- replicated SNPs, we have limited statistical power (,40%) to detect evidence of association. Allele frequency is different across ethnic groups, which may affect the confirmation. Among the 10 non-replicated SNPs, three showed very low MAF in our study. In addition, LD structure often differs across populations, and GWAS hits are typically markers in LD with functional alleles. Therefore, it is expected that many genetic variants identified in one population cannot be directly replicated in other populations. In addition, the effect size for the same SNP may differ across ethnic groups due to different environment exposure. This may also result in the non-replication of the SNPs in other ethnic groups. Furthermore, replication study with larger sample size would be warranted.
Strengths of the study include the focus on replicating significant GWAS-identified SNPs among Asian women and a relatively large sample size. Furthermore, we investigated SNPs in LD with the SNPs reported in European-ancestry but not in Asian-ancestry populations, which would enhance the genomic coverage of susceptible regions for ANM. Results observed in the present study should be interpreted considering several limitations. First, GWAS-identified SNPs were selected only if they met a GWAS significance threshold (P,5610 28 ), and some SNPs with less significant P-values may be ignored. Second, different studies used different genotyping platforms, and there are only a small portion of SNPs overlapping between these platforms. Some SNPs were directly genotyped in some parent studies but imputed in other parent studies. Third, no surrogate SNP was identified for rs16991615 and rs236114, and we did not investigate these two loci. Fourth, ANM was self-reported with measured error, although this error is unlikely to be associated with genotype.
In summary, we successfully replicated eight SNPs associated with age at natural menopause in Chinese women. This is the first report of the common shared genetic architecture affecting menopause between European and Chinese women. Further research on these genetic regions will enhance our understanding of the genetic determination of menopause.