Polygenic, enrichment and heritability analysis of three histone modification marks across cell types.

<p>We performed three different analyses to test for cell type specific effects on the genetic risk for MI/CAD. Analyses were conducted on SNPs residing in the three histone marks (H3K27ac, H3K4me3, H3K9ac) that were present in the different cell types. (A) Polygenic risk score analysis. We performed polygenic risk score association analysis on SNPs with MIGen discovery association <i>P</i><0.05. Negative logarithm of <i>P</i> values from association testing of the polygenic risk score performed in the WTCCC CAD was shown. Cell types were sorted based on the strength of polygenic association. Orange vertical line represents a significant level with 5% alpha error. (B) Enrichment of association. Enrichment analyses were performed by comparing the proportion of significant variants passing a specific association <i>P</i> threshold of a variant set with that of a baseline set. Different association <i>P</i> thresholds 5×10⁻⁷, 5×10⁻⁶, 5×10⁻⁵ from the CARDIoGRAM study were tested [<a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1005622#pgen.1005622.ref024" target="_blank">24</a>]. The variant sets in this analysis were SNPs in the specified histone marks that were present in the indicated cell type. For the baseline set, we test SNPs in regions that are outside of these histone marks within 10 kilobases (kb) of the protein coding regions of the genome. To reduce the effects of linkage disequilibrium, these baseline SNPs were selected to be 5 kb away from the histone marks. In the plot, each triangular point represents the strongest enrichment result for each mark in each cell type across the three possible association <i>P</i> thresholds. (C) Heritability analysis. Heritability analysis was performed within histone marks in the MIGen study. Each point in the plot represents the variance in liability generated from a joint model involving two variance components using the Genome-wide Complex Trait Analysis software [<a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1005622#pgen.1005622.ref022" target="_blank">22</a>,<a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1005622#pgen.1005622.ref023" target="_blank">23</a>]. The two variance components include 1) SNPs in the specified histone mark that was present in the indicated cell type and 2) all other SNPs outside of these regions. The variance in liability is an estimate from the ratio of genetic variance to phenotypic variance for the specified variance component (i.e. the specified variance component is all SNPs within the specified histone mark) whereas the <i>P</i> value is from the likelihood ratio test of a reduce model with the specified genetic variance component dropped from the full model, from the restricted maximum likelihood method in the Genome-wide Complex Trait Analysis software [<a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1005622#pgen.1005622.ref022" target="_blank">22</a>,<a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1005622#pgen.1005622.ref023" target="_blank">23</a>]. MI, myocardial infarction; CAD, coronary artery disease; SNP, single nucleotide polymorphism.</p

Hierarchical clustering of 45 MI/CAD GWAS SNPs and specific cell types for a histone modification mark (H3K27ac).

<p>We mapped 45 MI/CAD GWAS SNPs, as well as SNPs in high linkage disequilibrium (<i>r</i>²≥0.8), to H3K27ac in different cell types. Hierarchical clustering was based on the presence or absence of a SNP residing in H3K27ac in different cell types and was performed using the heatmap function in R (R Project for Statistical Computing). We observed unique patterns between the different GWAS loci and cell types. For example, 12 of the 45 GWAS loci were expressed in more than 80% of the cell types, whereas 13 of the 45 GWAS loci were expressed in less than 20%. Red color indicates a lead SNP or tag SNPs (linkage disequilibrium value of <i>r</i>²≥0.8) residing in H3K27ac in different cell types (See <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1005622#pgen.1005622.s009" target="_blank">S9</a> and <a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1005622#pgen.1005622.s010" target="_blank">S10</a> Figs for H3K9ac and H3K4me3, respectively). MI, myocardial infarction; CAD, coronary artery disease; GWAS, genome-wide association study; SNP, single nucleotide polymorphism.</p

Contributions of three genomic compartments to the polygenicity of MI/CAD.

<p>Polygenic risk score analysis was performed across three different genomic compartments. The top bar plot represents the strength of association for the polygenic risk score analysis whereas the bottom bar plot represents the number of SNPs within each of the compartments. The strongest polygenic association signals were within noncoding regions adjacent to protein-coding genes (“genic noncoding”). MI, myocardial infarction; CAD, coronary artery disease; SNP, single nucleotide polymorphism. Genic coding, variants that code amino acid sequence within ±10 kilobases of the 3′ or 5′ untranslated regions of a gene. Genic noncoding, variants that do not code amino acid sequence within ±10 kilobases of the 3′ or 5′ untranslated regions of a gene. Intergenic, variants that are beyond ±10 kilobases of the 3′ or 5′ untranslated regions of a gene.</p

Heritability of MI/CAD explained by three genomic compartment sets.

<p>Heritability estimates were inferred independently first in MIGen and WTCCC CAD from a single model involving three variance components (“genic coding”, “genic noncoding” and “intergenic”) using the GCTA software [<a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1005622#pgen.1005622.ref022" target="_blank">22</a>,<a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1005622#pgen.1005622.ref023" target="_blank">23</a>]. Heritability estimates shown here are from a meta-analysis of the Variance and standard error (V-SE) from these models using as weights the inverse variance from these models.</p><p>¹Variance and V-SE are estimates from the ratio of genetic variance to phenotypic variance for the specified variance component whereas the <i>P</i> value (V-P) is from the likelihood ratio test of a reduce model with the specified genetic variance component dropped from the full model, from the restricted maximum likelihood method in the GCTA software [<a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1005622#pgen.1005622.ref022" target="_blank">22</a>,<a href="http://www.plosgenetics.org/article/info:doi/10.1371/journal.pgen.1005622#pgen.1005622.ref023" target="_blank">23</a>].</p><p>²Enrichment of variance was calculated as the % variance of total divided by % SNPs of total. MI, myocardial infarction; CAD, coronary artery disease; SNP, single nucleotide polymorphism.</p><p>³<i>P</i> value from difference in the observed variance minus the expected variance (variance of whole genome as sum multiplied by % SNPs of total). Genic coding, variants that code amino acid sequence within ±10 kilobases of the 3′ or 5′ untranslated regions of a gene. Genic noncoding, variants that do not code amino acid sequence within ±10 kilobases of the 3′ or 5′ untranslated regions of a gene. Intergenic, variants that are beyond ±10 kilobases of the 3′ or 5′ untranslated regions of a gene.</p><p>We calculated the SNP-heritability in three genomic compartment sets for MI/CAD in a meta-analysis of the MIGen and WTCCC CAD studies using the Genome-wide Complex Trait Analysis (GCTA) software. We observed increased enrichment in variance in both “genic coding” and “genic noncoding” regions.</p

Genomic context of cis-regulated allele-specific methylation events.

<p>Illustrations showing genomic context and individual CpG methylation levels at three separate MPRs (rs10491434 (left), rs6569648 (middle) and rs943049 (right)). The chromosomal location of each amplicon is demonstrated with an ideogram. and the RefSeq genes (orange) surrounding the amplicon (red line) shown below. A section (grey box) contracts to the amplicon region itself to show the relative positions of the SNPs (black lines) and CpGs (blue lines) within the amplicons themselves (red bars); methylation levels of the alternate (grey circles) and reference (yellow circles) alleles within samples heterozygous for the SNP are graphed below each CpG. Asterixes (*) mark the CpGs illustrated within <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0098464#pone-0098464-g003" target="_blank">Figure 3</a>.</p

Candidate cis-regulated ASM variants in phenotypically implicated SNPs.

<p>Mean number of candidate cis-regulated ASM variants in ten random LD pruned (ld<0.3) sets of candidate cis-regulated ASM SNPs that are also in high LD (LD>0.8) with a Genome Wide Association Study (GWAS) derived variant from the National Human Genome Research Institute (NHGRI) database or an eQTL in monocytes or peripheral blood monocytes (PBMCs), with accompanying standard deviations of the mean (sd).</p

Types of allele-specific methylation candidates.

<p>Plots showing number of different categories of ASM candidates within the microarray study sample population. Of the 242533 MPRs for which there were at least 6 heterozygotes within the population (left pie chart) we detected some level of ASM in at least 116045 (left pie chart, green), of these we detected ASM in at least 6 samples in 12032 MPRs (left pie chart, dark green). Of these 12032 ASM MPRs, we detected cis-regulated ASM in 9750 MPRs (right pie chart, solid blue or solid yellow), and random or stochastic ASM in 2282 MPRs (right pie chart, mixed blue and yellow). Representation of patterns of allelic-choice in ASM within the microarray study sample population. ASM allelic choice is shown at 28 ASM and 2 non-ASM MPRs for the 42 samples in our initial microarray sample population. Non-heterozygous samples (white), samples with biallelic methylation (grey), and samples with ASM (blue and yellow) with methylation at either Allele-A (yellow) or Allele-A (blue) are shown. MPRs are organized in columns to show those determined to have no ASM (first two columns), cis-regulated ASM (both for Allele-A (3rd to 11th columns) and Allele-B (15th to final columns) or random ASM (12th to 14th columns).</p

Microarray based detection of allele-specific methylation.

<p><b>A)</b> A simplified representation of the Methyl-Sensitive Restriction Enzyme (MSRE) based Allele-Specific Methylation (ASM) assay. DNA is MSRE treated (left panels) and MSRE sites with methylated CpGs protected from digestion (upper panels, Allele-A) while its homologous chromosomal region with unmethylated CpGs are not (lower panels, Allele-B). The DNA is digested with StyI and NspI to form 200–1200 bp fragments, linkers ligated and DNA amplified to create amplicons that are hybridized to the array. Only regions with protected MSRE sites (methylated CpG) are amplified and can hybridize to show signal on the array (final panel). <b>B)</b> Bioinformatic detection of Allele-Specific Methylation (ASM) from Affymetrix SNP 6.0 arrays signals after MSRE digestion. In the scatter plot on the left, 4 different expected states after MSRE digest at a heterozygous region are compared to the typical distribution of probe intensities observed within the HapMap samples for the same MPR (here portrayed by light grey squares): biallelic methylation (dark grey circles), monoallelic A methylation (blue circles), monoallelic B methylation (yellow circles) and finally biallelic lack of methylation (red circles). The primary calling method relies on feature extraction by way of conversion of 2-dimensional A and B probe intensity data (scatter plot) from heterozygotes to log2(A/B) values and is compared against the typical log2(A/B) distribution observed for this MPR within the HapMap samples (histogram, light grey). Simply put, MPRs diverging from this distribution after MSRE treatment are called ASM. Using this method, biallelic unmethylated states have the potential to result in false positive ASM calls as any log2(A/B) value would be based on background noise, so are filtered out by removing MPRS with low total intensities (highlighted here with a red quarter-circle, for further information on how this filter was devised, see <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0098464#pone.0098464.s008" target="_blank">Figure S8</a>).</p

Confirmation of ASM in a subset of candidate cis-regulated ASM variants.

<p>Results from microarray and next-generation bisulfite sequencing ASM assays of ten variant-containing regions. The table shows the CpG position, whether it is found within an MSRE site and the Bonferroni adjusted chi-square based p-values of the association of methylation at this CpG with either the reference of alternate alleles. The allele with the highest number percentage of methylated reads was designated the most frequently methylated allele (REF = reference, ALT = alternate) and compared to that from the microarray data; for all CpGs that were within MSRE sites and showed significant association of methylation with an allele in the sequencing assay the methylated allele matched that of the microarray assay. For more details, see <a href="http://www.plosone.org/article/info:doi/10.1371/journal.pone.0098464#pone.0098464.s020" target="_blank">Table S3</a>.</p

Illustration of the GSE48091 gene-expression data-set used in Example-A (see main text).

<p>Each row corresponds to a patient, and each column to a ‘gene’ (i.e., gene-expression measurement): the color of each pixel codes for the intensity of a particular measurement of a particular patient (see colorbar to the bottom).<i>M</i>_<i>D</i> = 340 of these patients are cases, the other <i>M</i>_<i>X</i> = 166 are controls; we group the former into the case-matrix ‘<i>D</i>’, and the latter into the control-matrix ‘<i>X</i>’.</p