Y Chromosome

TBA [19] (link) alignments of the human genome (hg18) to 43 other vertebrate species were obtained from the UCSC genome browser [20] (link), [21] (link) together with a phylogenetic tree with the generally accepted topology (Fig S1) and neutral branch lengths estimated from 4-fold degenerate sites. Both the tree and alignments were projected to the 34 mammalian species. The alignment was compressed to remove gaps in the human sequence, and GERP++ scores were computed for every position with at least 3 ungapped species present, or approximately 88.9% of the 3.08 billion positions on the 22 autosomes and X/Y chromosomes. We used the HKY85 [13] (link) model of evolution with the transition/transversion ratio set to 2.0 and nucleotide frequencies estimated from the multiple alignment.
To limit memory requirements and allow parallelization of the constrained element computation, each chromosome was broken up into regions of approximately 2 megabases, with long segments where no RS score was computed chosen as boundaries. These boundary segments contain no information usable by GERP++ and because the algorithm never annotates constrained elements spanning them, excluding such segments did not sacrifice any predictive ability. These boundary regions made up approximately 6.8% of the human genome, including a 30.2 megabase region that made up more than half of chromosome Y. Constrained element predictions were generated using default parameters and a 5% false positive cutoff measured in terms of number of predictions; the estimated nucleotide-level false positive rate was under 1%. As additional validation, we computed overlap between our predictions and a set of ancestral repeats (L2) annotated by RepeatMasker. We found the overlap to be in line with what we expected given our estimated false positive rates: about 5% of the repeats overlap a predicted CE, with around 1.6% nucleotide-level overlap.
Gene, noncoding RNA, and PhastCons conserved element annotations were obtained from the UCSC genome browser's [20] (link), [21] (link) Known Genes [22] (link), RNA Genes, and Conservation [4] (link) tracks respectively. To avoid skewed statistics due to alternative splicing, gene annotations were resolved to a consistent nonoverlapping set where any segment belonging to multiple conflicting annotations was assigned a single annotation in the following order of priority: coding exon, 5′ UTR, 3′ UTR, intron. For meaningful comparison against phastCons, separate GERP++ scores and constrained elements were generated according to the same procedure as above but using only placental mammal data (ignoring platypus and opossum in the alignment and projecting them out of the phylogenetic tree).
PolII binding regions were defined as 50 bp upstream and downstream of PolII binding ‘peaks’ as identified from ChIP-seq experiments performed by the ENCODE Consortium [3] (link). A 100 bp window allows capture of the likely PolII binding site and its flanking sequence. We obtained data from nine ChIP-seq experiments conducted in two labs (the Snyder lab at Yale and the Myers lab at Hudson Alpha) on six cell types. Data was downloaded through the DCC at UCSC (ftp://encodeftp.cse.ucsc.edu). All data have passed publication embargo periods. Overlap statistics were calculated as described above for other annotation sets and averaged across all nine experiments.

DMRcate is the R package implementation of this method, and is available from the Bioconductor repository [69 ]. All the user needs to provide to the DMRcate workflow is an Illumina probe ID-indexed matrix of methylation measurements and, for DMR finding, a model matrix (optionally with an additional contrast matrix) reflecting the experimental design. DMRcate’s experimental design idiom is lifted wholesale from limma, hence it can fit any given contrast from a limma model. The workflow assumes that the user has already normalised the data according to their preferred method, and removed bad-quality probes via detection P values and low bead count. As a further preprocessing option, DMRcate provides an optional filtering function, removing probes whose reported methylation level may be confounded by SNPs, and/or by cross-hybridisation [45 (link)]. Further investigation of a 450K data set used in this study reveals that SNPs within two nucleotides of the target CpG site have a perturbed beta distribution (Additional file 1: Figure S1). If the samples are from mixed-sex groups, the option of removal of probes hybridising to X and Y chromosomes is also provided.
Initial output consists of a data frame describing each region, ranked by its corresponding P value. Useful information such as genomic coordinates, gene associations and number of constituent CpGs per region are also reported. The user may specify any positive bandwidth they like. Longer bandwidths allow for interrogation on a broader, even chromosomal, scale, while shorter bandwidths potentially allow identification of focal regions of DM. Post-fitting, the user has the option of filtering out any region that does not have at least one constituent CpG site with a beta fold change greater than a specified threshold. As an alternative to using genomic coordinates, DMRcate has a consecutive option that assumes all assayed CpGs are equally spaced. Wrappers for GenomicRanges object and whole-genome BedGraph production are provided. For visualisation, a separate plotting function for individual DMRs is also provided. An example of the graphical output of DMRcate can be found in the online vignette, and as part of Additional file 1: Figure S5. A complete description of DMRcate’s functionality and user options is available from the manual [70 ].

IDAT files were loaded into the R (2.14) environment using the Bioconductor (2.9) minfi package (1.0.0) [25 ]. The detection P-values for all probes were then calculated for the data using functionality provided in minfi. Probes on the × and Y chromosomes were removed at this stage. Two versions of the data were used in subsequent analyses: the raw data and SWAN data. Probes with a detection P-value >0.01 in one or more samples were then excluded. The differential methylation analysis was performed for both datasets on the subset of 18,678 probes that overlapped with the RRBS data using the 'dmpFinder' minfi function. The 'dmpFinder' function uses an F-test to identify positions that are differentially methylated between two groups. The tests are performed on M-values (log₂(Methylated/Unmethylated)) as recommended in Du et al. [35 (link)]. Variance shrinkage was used due to the small sample size. In 'dmpFinder', the sample variances are squeezed by computing empirical Bayes posterior means using the limma package [36 ]. Example R code for performing a differential methylation analysis using minfi can be found in Additional file 2.
True positives were defined to be CpGs that had an absolute difference in β value >0.25 between the kidney and rectum RRBS samples. Additionally, for the ROC analysis, which was performed using the ROCR package [31 (link)], true negatives were defined as those CpGs found to have an absolute difference in β value <0.05 between the RRBS samples.

Data from the consensus coding sequence (CCDS) database (release 22) [63 (link)] were used to estimate the sizes of exonic, intronic and intergenic regions for individual chromosomes. First we retrieved the start and end exon and intron positions of each coding gene transcript. The sum of the lengths of all exons on a given chromosome was used as a size of the exonic region. We excluded overlapping gene sequences to make sure that each nucleotide was counted only once. The similar approach was used to estimate the total size of intronic regions. Noncoding genes were counted as part of the genic region. Start and end positions of noncoding genes were retrieved from the NCBI Gene database [64 (link)]. The size of the intergenic region was computed as the total size of the chromosome minus the combined size of exons, introns and noncoding genes. Genes located on the Y-chromosome and mitochondrial genes were excluded from the analysis.

The Drosophila melanogaster Oregon R (modENCODE, FBst0025211) and Phs-hidY (FBst0024638; Grether et al. 1995 (link)) fly stocks were obtained from the Bloomington Drosophila Stock center (Cook et al. 2010 (link)). The Drosophila melanogaster Canton S flies were obtained from the Bloomington Drosophila Stock center (RRID:BDSC_64349). Flies were grown on either Cornmeal, Sucrose, and Yeast Media (LabExpress fly media) or Glucose, Active Yeast, and proprionic acid media (10:5 fly media) at 25 ${}^{\circ }$ C, 60% RH. 5 females and 4 males were placed per vial for Oregon R and approximately 10 female and male flies per vial for Canton S. Oregon R flies were transferred to fresh vials daily and progeny from previous vials were collected for experiments. Canton S flies were transferred to fresh vials every 3–4 days and progeny were selected for experiments. An Oregon R stock with a Y chromosome carrying a Phs-hid insertion was generated (http://flystocks.bio.indiana.edu/Browse/misc-browse/hs-hid˙method.html), by crossing double balancer females with Y-lethal males crossed with Y-lethal males (FBst0024638, Grether et al. 1995 (link)). We introduced the Y chromosome into the Oregon R background using double balanced 2nd and 3rd chromosomes. The 4th chromosome was not tracked. For our experiments we used Oregon R (FBst0025211) males and virgin females obtained from the Oregon R Y-lethal stock after 37 ${}^{\circ }$ C heat shock using a water bath for 2 h during third instar. Canton S flies sex were sorted after mating, for female only experimentation. We verified the absence of males during collection of individual flies. When males were present in the vial, none of those flies were used for experiments. This occurred in about 1 out of every 10 heat shocked vials. Newly enclosed flies were placed in batches of 35 per vial for use in experimental assays. Vials contained either CDF (Lee and Micchelli 2013 (link)), LabExpress fly media, or 10:5 fly media . Flies were aged for 3–4 days at 25 ${}^{\circ }$ C 60% RH. After aging on solid food, flies were loaded in the WAFFL (as described below) to proceed with the defined assays.

Macrophages from male fetuses were identified on the basis of sum log Y chromosome expression above a threshold (0 in mouse and 0.4 in human). For mouse data, differentially expressed genes were defined between male fetal macrophages and all other cells, and any significant genes were used for hierarchical clustering of all cells into two groups. The group with an overrepresentation of male fetal macrophages was defined as fetal. To measure fetal correlation as a confirmation of identity of non‐Y chromosome‐expressing macrophages assigned to the fetal cluster, expression of each cell was compared with the average expression of the same set of differentially expressed genes across all confirmed male fetal macrophages. This analysis was repeated independently on all macrophages included in analysis (from Sham‐ KO‐ and WT‐infected conditions) at each time point to identify a fetal expression profile unique to that time point. Because analysis at each time point included cells processed from either 6 or 8 placentas, each fetal analysis included some cells from a male fetus to provide a ground truth at that time point. For human data, clusters were defined by hierarchical clustering on informative genes, defined as highly expressed and variable genes (fano‐factor and mean expression) and the group with an overrepresentation of male fetal macrophages was defined as fetal.