International genome-wide meta-analysis identifies new primary biliary cirrhosis risk loci and targetable pathogenic pathways.

Primary biliary cirrhosis (PBC) is a classical autoimmune liver disease for which effective immunomodulatory therapy is lacking. Here we perform meta-analyses of discovery data sets from genome-wide association studies of European subjects (n=2,764 cases and 10,475 controls) followed by validation genotyping in an independent cohort (n=3,716 cases and 4,261 controls). We discover and validate six previously unknown risk loci for PBC (Pcombined<5 × 10(-8)) and used pathway analysis to identify JAK-STAT/IL12/IL27 signalling and cytokine-cytokine pathways, for which relevant therapies exist.

Primary biliary cirrhosis (PBC) is a rare cholestatic liver disease characterized by progressive autoimmune destruction of intrahepatic bile ducts, leading to cirrhosis and liver failure in a substantial proportion of cases1. To date, four genome-wide association studies (GWAS) and two Illumina immunochip studies of PBC have confirmed associations at the human leukocyte antigen (HLA) locus and identified 27 non-HLA risk loci2345678. Consistent with GWAS data for other autoimmune diseases, results of these studies implicate immune-related genes in disease pathogenesis, but in general fail to pinpoint the disease-causal variants within the identified risk loci. To identify risk alleles that may be relevant to disease biology and treatment and illuminate additional PBC risk loci, we undertook a genome-wide meta-analysis (GWMA) combining North American, Italian and UK PBC GWAS data sets245. Functional annotation of the risk loci and pathway analyses were then performed to identify the alleles and pathways most relevant to disease cause and treatment.

Results

Discovery of new PBC risk loci

Following quality control, the combined discovery data set for GWMA consisted of 1,143,634 genotyped or imputed single-nucleotide polymorphisms (SNPs) in 2,764 cases and 10,475 controls (Supplementary Table 1). After genomic control correction and exclusion of known PBC risk loci from the final set of results, the inflation factor was λ=1.043 (Supplementary Fig. 1). Meta-analysis of this data set identified 23 loci at genome-wide level of significance (P<5 × 10−8, calculated using logistic regression of individual discovery data sets in ProbABEL followed by genomic control correction of individual discovery data sets in R and fixed-effects meta-analysis in META, see Methods). Of these, 22 had been detected in previous studies and the 23rd corresponded to a most-likely spurious signal from a single imputed SNP on chromosome 13 (Supplementary Fig. 2, Supplementary Table 2). However, we found suggestive evidence of association (P<2 × 10−5 from fixed-effects meta-analysis in META) at 41 loci not previously known to be associated with PBC. The top-scoring SNPs (or close proxies in strong linkage disequilibrium with the top-scoring SNP) from these and nine other loci (including the likely spurious chromosome 13 signal) were taken forward for genotyping in an independent panel consisting of 3,716 PBC cases and 4,261 controls. In total, 120 SNPs at 50 independent loci were taken forward for validation, of which 114 were successfully genotyped (Supplementary Data 1).

In the validation analysis, we confirmed association with SNPs at six loci not previously known to be associated with PBC (P<4.4 × 10−4, equivalent to P=0.05 with a Bonferroni correction for 114 tests, calculated using logistic regression analysis of individual validation data sets in PLINK followed by meta-analysis in META, see Methods); meta-analysis of discovery and validation cohorts at these loci reached genome-wide levels of significance (Pcombined <5 × 10−8 from fixed-effect meta-analysis in META) (Table 1, Supplementary Figs 3 and 4). Furthermore, SNPs at two additional loci achieved P values suggestive of association (P<1 × 10−3 from fixed-effect meta-analysis in META, equivalent to P=0.05 with a Bonferroni correction for testing at 50 independent loci; Table 1, Supplementary Fig. 5). Newly identified PBC risk loci overlap with those of other autoimmune disorders and harbour several immunologically relevant candidate genes, most notably chemokine ligand 20 (CCL20) and interleukin 12B (IL12B; Table 1).

Table 1

PBC risk loci identified in the current study.

a. Confirmed risk loci (validation P<4.4 × 10−4 resulting in combined P<5 × 10−8)
Locus	SNP	Position (build 38)	A1/A2	DiscoveryP	ValidationP	JointP	OR (95% CI)	Region (build 38)	Nearby genes and functional annotation*	Autoimmune overlap
2q12.1	rs12712133	102,249,813	A/G	1.62 × 10−5	7.94 × 10−5	5.19 × 10−9	1.14 (1.07–1.21)	102,118,975–102,438,307	IL1R1, IL1RL2†, FAM183DP, IL1RL1‡, IL18R1, LOC100422339, IL18RAP, MIR4772	CD, CeD
2q36.3	rs4973341	227,795,646	C/T	6.48 × 10−7	7.73 × 10−5	2.34 × 10−10	0.82 (0.74–0.90)	227,747,828–227,815,647	RNA5SP121, SNRPGP8, LOC100533842, CCL20†,§
4p16.3	rs11724804	971,991	A/G	3.67 × 10−7	4.25 × 10−6	9.01 × 10−12	1.22 (1.12–1.33)	853,681–1,014,424	GAK, TMEM175, DGKQ†, SLC26A1a, IDUA, FGFRL1
5q21.1	rs526231	103,345,680	T/C	3.10 × 10−5	9.39 × 10−5	1.14 × 10−8	0.87 (0.81–0.93)	102,939,698–103,416,571	PAM§, EIF3KP1, GIN1, PPIP5K2, C5orf30‡,§	RA
5q33.3	rs2546890	159,332,892	G/A	1.20 × 10−6	1.89 × 10−5	1.06 × 10−10	0.87 (0.82–0.93)	159,117,927–159,414,310	RNF145, UBLCP1, RNU4ATAC2P, IL12B, LOC285626†	Pso, MS, CD
6q23.3	rs6933404	137,638,098	C/T	9.47 × 10−7	2.84 × 10−5	1.27 × 10−10	1.18 (1.09–1.27)	137,571,557–137,803,754	LOC102723649, LOC442263, OLIG3†, TNFAIP3†	RA, SLE, SjS, CeD, UC, MS

A1, tested allele; CD, Crohn disease; CeD, coeliac disease; CI, confidence interval; MS, multiple sclerosis; OR, odds ratio in validation cohorts; Pso, psoriasis; RA, rheumatoid arthritis; SjS, Sjogren syndrome; SLE, systemic lupus erythematosus; SNP, single-nucleotide polymorphism; UC, ulcerative colitis.

PBC risk loci identified in the current study. SNPs were taken forward for validation based on having a discovery P value <2 × 10 (or, in the case of rs526231 and rs2434360, based on acting as a proxy for a SNP with a P value <2 × 10). Discovery P values were calculated using logistic regression of individual discovery data sets in ProbABEL followed by genomic control correction of individual discovery data sets in R and fixed-effects meta-analysis in META; validation P values were calculated using logistic regression of individual data sets in PLINK followed by fixed-effect meta-analysis in META; joint P values were calculated using fixed-effect meta-analysis of discovery and validation data sets in META; see Methods. Autoimmune overlap refers to overlap between risk loci for PBC and those of other autoimmune conditions.

*Functional annotation.

†Regulatory variants: The index SNP or variants in strong linkage disequilibrium (LD, r2≥0.8) with the index SNP at this locus overlap regulatory elements that are related to the annotated gene (Supplementary Table 3).

‡mQTLs: The index SNP or variants in strong LD are correlated to methylation related to the annotated gene (Supplementary Data 4).

§eQTLs: The index SNP or variants in strong LD are correlated to expression of the annotated gene (see Supplementary Data 3).

Discovery of candidate causal disease variants

In functional annotation of risk loci, we identified 199 candidate variants across 28 non-HLA risk loci with probabilistic identification of causal SNPs (PICS) probability >0.0275 (ref. 9). At each risk locus, the most-likely causal variant was the index variant, with median PICS probability of 0.224 and values up to 0.998 for rs2546890 at 5q33.3 (Supplementary Data 2). Looking at all candidate variants across all risk loci, the majority were intronic, upstream or downstream gene variants with no predicted functional consequence (99/199, 40%). However, a substantial proportion (59/199, 30%) were regulatory region variants, defined as SNPs located within regulatory features, including enhancers, promoters, transcription factor-binding sites and open chromatin regions (Supplementary Data 3). Notably, candidate variants at 18 (64%) of the 28 annotated risk loci included at least one regulatory region variant. In contrast, only 5 of 199 candidates were missense variants (2.5%) (Supplementary Table 3a). However, these included rs2297067 in EXOC3L4 at 14q32.32 and rs2304256 in TYK2 at 19p13.2, both predicted by SIFT and/or PolyPhen to be deleterious or potentially damaging1011. Candidate variants included a single splice region variant, that is, rs17641524 at 1q31.3 that is predicted to affect splicing of DENND1B (Supplementary Table 3b).

We found that candidate variants at several risk loci are methylation quantitative trait loci (mQTLs), including mQTLs for DENNDIB, PLCL2, IRF5 and TNFRSF1A, all genes that are implicated in risk for other autoimmune diseases (Supplementary Data 4). We also found that candidate variants at several risk loci are expression quantitative trait loci (eQTLs) in lymphoblastoid and other cell lineages, including eQTLs for CCL20, IL12A, IRF5 and TYK2 (Supplementary Data 5).

At many risk loci, functional annotation highlighted a single candidate gene (Supplementary Data 2). However, most risk loci contained multiple compelling candidate variants. This complexity is well illustrated by the composite of candidate variants at the PLCL2 gene and MANBA gene loci, which include multiple eQTL and mQTL SNPs. Thus, despite the presence of many candidate variants with regulatory or epigenetic roles within PBC risk loci, more direct biological experimental approaches are required to pinpoint the disease-causal variants at these loci.

We also applied functional GWAS (FGWAS) and its associated annotation file12 to our full set of discovery GWMA results and thereby identified 75 annotations with enrichment (P<0.01 from FGWAS) of GWMA association signals (Supplementary Data 6). After a stepwise selection approach similar to that of Pickrell12, the best-fitting model included six annotations highlighting negative enrichment of repressed chromatin regions in a lymphoblastoid cell line, and positive enrichment of DNase-I-hypersensitive sites in a variety of cell types, in particular CD20+ and Th1 T cells (Supplementary Table 4).

Identification of candidate targetable biological pathways

To identify biological pathways involved in development of PBC, we conducted pathway analysis using GCTA13 followed by i-GSEA4GWAS14. We identified several immunoregulatory pathways associated with PBC, in particular, IL-12 and other cytokines as well as T-cell signalling pathways. To account for bias that might result from the strong HLA association with PBC, we repeated this analysis with SNPs/genes in the HLA region excluded. Notably, IL-12, IL-27 and JAK-STAT signalling pathways were still associated with PBC, even after their HLA contribution had been removed (Table 2).

Table 2

Results from pathway analysis in iGSEA4GWAS.

Gene set	Source	FDR (with HLA)	FDR (without HLA)
NO2-dependent IL-12 pathway in NK cells	Biocarta	6.7 × 10−4
JAK-STAT signalling pathway*,†	KEGG	0.001	0.013
IL-12 mediated signalling events	PID	0.001
IL-12- and Stat4-dependent signalling in Th1 development*,†	Biocarta		<0.001
Interferon signalling	REACTOME	0.001
PD-1 signalling	REACTOME	0.001
Phosphorylation of CD3 and TCR-ζ chains	REACTOME	0.001
IL-27-mediated signalling events*,†	PID	0.001	<0.001
Cytokine–cytokine receptor interaction‡,§,||,¶,#,**,†,††,‡‡	KEGG	0.002	0.010
IFN-γ signalling	REACTOME	0.002
MHC class II antigen presentation	REACTOME	0.003
Cytokine signalling in immune system	REACTOME	0.004
Antigen processing and presentation	KEGG	0.004
Intestinal immune network for IgA production	KEGG	0.004
Co-stimulation by the CD28 family	REACTOME	0.005
IL-2 mediated signalling events	PID	0.008
TCR signalling	REACTOME	0.008
Downstream TCR signalling	REACTOME	0.008
Cell adhesion molecules	KEGG	0.015
Th1, Th2 differentiation†	Biocarta	0.019	0.012
IL-2 receptor beta chain in T-cell activation	Biocarta	0.021
Interferon α/β signalling	REACTOME	0.035
IL-23-mediated signalling events	PID	0.039

FDR, false discovery rate; IFN, interferon; IL, interleukin; PD, programmed cell death; TCR, T cell antigen receptor; NK, natural killer.

Gene sets with FDR <0.05 are listed. Results are shown with or without inclusion in the analysis of SNPs within the HLA region. The top 10 hits from our drug-positioning analysis using a combined pathway from the HLA excluded set are indicated by symbols for the associated pathways that they affect.

*Tofacitinib.

†Glatiramer acetate.

‡Axitinib.

§Pazopanib.

||Vatalanib.

¶Cediranib.

#X-82.

**Telatinib.

††Linifanib.

‡‡Tandutinib.

We identified molecules that targeted these pathways by overlaying the Drug Gene Interaction database15 and calculating a pathway specificity score and Jaccard index of each drug for each of the pathways that remained associated with PBC after the HLA contribution had been removed (Table 2, Supplementary Data 7). This combined analysis identified pathways and immunomodulatory agents that represent promising leads for further study in models of PBC.

Discussion

The current study adds to our knowledge of the genetic architecture of PBC. Notably, our data identify CCL20 as a candidate risk gene for PBC. Chemokine ligand 20 (CCL20) and its chemokine receptor CCR6 contribute to the formation and function of mucosal lymphoid tissues and are notably, in the context of the immune-mediated lymphocytic cholangitis characteristic of PBC, involved in the localization of Th17 cells and CD8 effector T cells to cholangiocytes and the periductal area in portal tracts16. This study also reinforces the importance of IL-12 and JAK-STAT signalling in this disease.

The functional annotation of risk loci has helped to assign priority to the candidate genes at newly identified and established risk loci. Furthermore, the identification of disease-associated regulatory variants at multiple risk loci emphasizes the potential importance of gene regulation in the pathogenesis of PBC (and presumably other complex disorders). This possibility is corroborated by the finding of numerous risk loci wherein the index and/or closely related SNPs that appear to represent regulatory, mQTL and/or eQTLs variants related to the nearby gene. Via the FGWAS analysis, this study also suggests particular importance of CD20+ B cells and Th1 cells in the pathogenesis of PBC. However, both the cell types and the specific gene variants most relevant to PBC require further investigation and in particular exploration of the tissue-specific functional effects of the disease-associated variants.

By looking for drug–gene interactions, we have identified candidate drugs targeting specific, PBC-associated pathways, creating new opportunities to re-purpose available drugs for targeted immune therapy. Despite the speculative nature of this analysis, the data provide a start point in the search for novel therapies that are urgently needed to improve outcomes for PBC patients.

Methods

Study samples and genotyping

The use of human subjects for this study was approved by the University Health Network Research Ethics Board, The Mayo Clinic Institutional Review Board, Etico Indipendente IRCCS Istituto Clinico Humanitas, UC Davis Institutional Review Board and the Oxford Research Ethics Committee.

All PBC cases included in the Canadian–US, Italian and UK discovery and validation cohorts fulfilled the American Association for the Study of Liver Diseases criteria for PBC.

The Canadian–US discovery cohort included 499 PBC cases who were self-reported whites of European descent and 390 healthy Canadian controls, all genotyped using the Illumina HumanHap370 BeadChip. Additional controls included in this cohort were 1,094 control subjects provided from the Prostate Cancer Genetics Markers Susceptibility (CGEMS), 1,142 controls from the Breast CGEMS studies and 1,748 controls from the New York Cancer Project, all of whom who were genotyped on an Illumina 550 K bead array4. Following all quality control (QC) procedures, the final Canadian–US discovery set included 499 PBC cases and 4,374 controls.

The PBC cases included in the Italian discovery cohort were self-reported whites of Italian descent genotyped using the Illumina Human610-Quad BeadChip. Controls in this cohort were healthy Italians genotyped using the Illumina 1M-duo array. Following QC procedures, the final Italian discovery set comprised 449 cases and 940 controls.

The PBC cases included in the UK discovery cohort were self-reported whites of British descent, genotyped using the Illumina Human-660 W Quad array. Controls in this cohort were 5,163 population controls genotyped on the Illumina 1M-Duo array as part of the Wellcome Trust Case Control Consortium 2 project. Following QC procedures, the UK discovery set comprised 1,816 cases and 5,161 controls.

The ‘Canadian' 903 PBC cases and 834 controls included in the validation studies were self-reported whites of European descent recruited from Canada, Europe and the United States to an ongoing PBC genetics study based in Toronto. The 721 ‘US' PBC cases and 294 controls included in the validation studies were self-reported whites of European descent enroled in the Mayo Clinic PBC Genetic Epidemiology registry and biorepository based at the Mayo Clinic in Rochester (https://clinicaltrials.gov/ct2/show/NCT01161953?term=pbc&rank=5).The Italian PBC cases and controls included in the validation studies were self-reported whites of Italian descent recruited to the Italian PBC Genetics study based at Instituto Humanitas in Milan. The Italian controls were obtained from Ospedale Alessandro Manzoni, Lecco, Italy and were unrelated healthy volunteers with no known non-Italian heritage. Cases and controls from the Canadian, Italian and the US cohorts were genotyped at the University Health Network/Mount Sinai Hospital Clinical Genomics Centre using a Sequenom iPLEX Gold assay. Following QC procedures, the final validation set included 903 cases and 834 controls from Canada; 300 cases and 618 controls from Italy; and 721 cases and 294 controls from the United States (Supplementary Table 1).

The ‘UK' PBC cases included in the validation studies were self-reported whites of British descent recruited to the UK-PBC project via the UK-PBC Consortium (http://www.uk-pbc.com/). Cases were genotyped using Sequenom iPLEX Gold assay at the Wellcome Trust Sanger Institute Genotyping Facility (http://www.sanger.ac.uk/). The UK validation control data were obtained from the TwinsUK resource, an adult twin registry comprising 12,000 (predominantly female) British twins. Genotype data for 3,512 twin individuals (genotyped using the Illumina HumanHap610 array) were obtained from the Department of Twin Research and Genetic Epidemiology at King's College London. One twin from each genotyped pair was included in the current study, amounting to 2,603 unrelated individuals. Following QC procedures, the final UK validation set comprised 1,792 PBC cases and 2,515 TwinsUK controls (Supplementary Table 1).

Quality control

We implemented a standard QC pipeline across all three discovery data sets, over-and-above QC procedures carried out in the respective primary analyses245. QC checks were carried out using the software package PLINK17. Within each discovery data set we removed SNPs with a genotype call rate <95%; minor allele frequency <0.05; significant deviation from Hardy Weinberg Equilibrium in controls (P<10−5) or a large difference (>5%) in the proportion of missing genotypes in cases versus controls. We removed samples showing high rates of missing data (>90%); whole-genome heterozygosity >six s.d. from the mean; estimated proportion of identity by descent (IBD) sharing with another sample >0.1, or apparent gender discrepancies (based on X-chromosomal heterozygosity >0.2 for men and <0.2 for women). Principal component analysis (based on a subset of 32,000 highly informative SNPs) was carried out using the ‘smartpca' routine of the EIGENSOFT package18 to identify population outliers for exclusion and to identify principal components that differed between cases and controls; these principal components were used as covariates in subsequent association analyses.

Genome-wide imputation

We used the SNPs and samples passing QC to carry out genome-wide imputation within each of our cohorts using the software package MaCH19 with HapMap3 CEU+TSI samples as reference data sets. Within each cohort we used approximately the same set of genotyped SNPs in cases and controls to ensure similar levels of informativity. Following imputation, we retained only those SNPs displaying minor allele frequency >0.005 and imputation quality score R2>0.5 in all three cohorts.

Statistical analysis of discovery cohorts

Within each cohort we carried out association analysis of the genome-wide imputed data allowing for imputation uncertainty using the software package ProbABEL20. We performed logistic regression of disease phenotype on allele dosage; principal components that differed between cases and controls were included as covariates to help correct for population stratification. Quantile–quantile plots of the genome-wide set of test statistics were examined and genomic control correction was carried out within each cohort by multiplying the standard error of the estimated log odds ratio for each SNP by the square root of the genomic control inflation factor λ (ref. 21). The resulting log odds ratios and adjusted standard errors from all three cohorts were meta-analysed using the software package META to produce the final set of genome-wide discovery results22.

Validation analysis

We selected loci for validation if they achieved suggestive level of significance in the discovery analysis (minimum P<2 × 10−5) and were not already known to be associated with PBC. We also selected loci for validation if they had achieved genome-wide significant association in one previous study but had never been validated in an independent cohort. We selected approximately two validation SNPs per locus; for loci displaying extended patterns of linkage disequilibrium or harbouring several putative independent association signals we attempted to select two validation SNPs within each subregion.

Within each locus chosen for validation we assigned priority to SNPs according to whether they had been genotyped in the TwinsUK cohort (which was used as a validation cohort for the UK validation cases). One SNP selected for validation (rs2297067) did not have genotype data available in TwinsUK and was therefore imputed within TwinsUK based on genotyped SNPs in the surrounding 5-Mb region using the software packages SHAPEIT23 and IMPUTE2 (ref. 24), with 1,000 Genomes (Phase I version 3 integrated data, released on March 2012) used as a reference sample. The TwinsUK cohort was subjected to a variety of additional QC checks as described previously25; the 2,515 controls used here correspond to the 2,520 controls used previously with an additional five exclusions due to discrepant gender25.

Within each validation cohort we carried out case/control association analysis of those SNPs that were successfully genotyped using logistic regression in PLINK. Results from the four validation cohorts (or from the combined discovery and validation cohorts) were combined via meta-analysis in META.

Imputation to 1,000 Genomes within validated loci

Imputation within the discovery cohorts was carried out at the six validated loci using the software packages SHAPEIT23 and IMPUTE2 (ref. 24) with the 1,000 Genomes (Phase I integrated variant set, release December and June 2013) used as a reference panel. The same genotyped SNPs that had been used to inform HapMap3 imputation for the discovery analysis were used for the 1,000 Genomes imputation within these targeted regions. Association analysis of SNPs passing post-imputation QC (‘info' score >0.5) was carried out separately within each cohort, the results were genomic control corrected by multiplying the standard error of the estimated log odds ratio for each SNP by the square root of the previously estimated genomic control inflation factor λ for each cohort, and results were combined across the cohorts via meta-analysis in META. This confirmed the findings from our original (HapMap3) imputation experiment but did not identify any substantially stronger associations or candidate causal variants than we had already found.

Functional annotation of validated loci

Left and right boundaries for each associated region were defined by finding a 0.1-cM interval either side of the most strongly associated SNP where no SNP has P<1 × 10−5. We looked for overlap between PBC risk loci and confirmed risk loci for other autoimmune conditions using ImmunoBase, a web-based resource focused on the genetics and genomics of immunologically related human diseases (http://www.immunobase.org/). To assign priority to candidate genes and candidate variants at risk loci, we used the online PICS (Probabilistic Identification of Causal SNPs) algorithm to identify candidate variants at each risk locus with a PICS probability >0.0275 (http://www.broadinstitute.org/pubs/finemapping/?q=pics)9. We adopted this threshold to be consistent with Farh et al.9 in their paper describing the approach. Given an index SNP corresponding to the most associated SNP in a locus, the PICS algorithm calculates (based on the known linkage disequilibrium pattern in the region, as measured in a large Immunochip or 1000 Genomes reference sample) a score for each SNP in the region, representing the extent to which that SNP could, in fact, be the true causal SNP, allowing for statistical sampling variation.

We then used the Ensembl Variant Effect Predictor web tool to annotate candidate variants for their predicted functional consequences (http://www.ensembl.org/info/docs/tools/vep/index.html). We used Genevar to evaluate the measured effects of these variants on DNA methylation in tissue collected from 856 healthy female twins of the MuTHER resource (http://www.sanger.ac.uk/resources/software/genevar/)2627. We used Genevar26, seeQTL (http://www.bios.unc.edu/research/genomic_software/seeQTL/)28 and the University of Chicago eQTL browser (http://eqtl.uchicago.edu/cgi-bin/gbrowse/eqtl/) to identify eQTLs amongst candidate variants.

We also used the FGWAS software and its associated annotation file (containing 450 genomic annotations of various types), applied to our full set of GWMA results, to investigate the extent to which genetic variants associated with PBC were enriched within specific annotation categories12. Testing each annotation individually, we found 75 annotations that showed enrichment (P<0.01) of GWMA association signals; as many of these annotations are correlated with one another we used a stepwise selection approach followed by cross-validation to mitigate overfitting (similar to the procedure performed by Pickrell12) on these 75 annotations to identify a final best-fitting model that included 6 annotations. Annotation information used by FGWAS was derived from a variety of sources including Maurano et al.29, Thurman et al.30 and Hffman et al.31 (see Appendix of Pickrell12 for details).

Pathway analysis

Using summary results from the GWMA (effect size, standard error and allele frequency) along with SNP linkage disequilibrium estimated from the Italian GWAS individual-level genotype data, we performed approximate conditional analysis using the software GCTA13. Only the independently associated signals with conditional P value and PGWMA both <0.001 were retained for further consideration. We submitted the rsIDs and PGWMA of these SNPs as well as gene sets from BioCarta, KEGG, PID and Reactome curated by MSigDB (as of 26 March 2014) to the i-GSEA4GWAS web server14. This programme identified genes within 20 kb of the SNPs and represented each gene by the greatest –log PGWMA of the SNP(s) mapped to it. Gene sets were then assessed for enrichment with significant genes while SNP label permutations were conducted to correct for bias from variations in gene size and gene set size. False discovery rate was used to correct for multiple testing based on the distributions of enrichment scores generated by permutation.

Drug-pathway analysis

To identify drugs that affected the pathways associated with PBC (when the HLA locus was excluded), we first identified the genes participating in each pathway. We then downloaded drug–gene associations from the Drug Gene Interaction database15 and scored each drug by the proportion of each its targets that were in each pathway, which we termed as the drug's pathway specificity. As a secondary scoring metric, we evaluated the proportion of each pathway affected by the drug using the Jaccard index on the respective sets of pathway genes and targeted genes. To identify promising drug candidates, we ranked drugs first by our primary specificity metric and then by the secondary Jaccard index.

Additional information

How to cite this article: Cordell, H. J. et al. International genome-wide meta-analysis identifies new primary biliary cirrhosis risk loci and targetable pathogenic pathways. Nat. Commun. 6:8019 doi: 10.1038/ncomms9019 (2015).