MiDAS tutorial

Are classical HLA alleles associated with disease status?

The following function prepares our data for analysis, combining HLA and phenotypic data into one object. Here, we are interested in analyzing our data on the level of classical HLA alleles.

HLA <- prepareMiDAS(
  hla_calls = dat_HLA,
  colData = dat_pheno,
  experiment = "hla_alleles"
)

We can now test our HLA data for deviations from Hardy-Weinberg-Equilibrium (HWE) to filter out alleles that strongly deviate from HWE expectations (imputation or genotyping errors, …). Here, let’s remove alleles with a HWE p-value below 0.05 divided by the number of alleles tested / present in our data (N=447). We can create a filtered MiDAS object right away (as.MiDAS = TRUE), as done in this example, or output actual HWE test results.

HLA <- HWETest(
  object = HLA,
  experiment = "hla_alleles",
  HWE_cutoff = 0.05 / 447,
  as.MiDAS = TRUE
)

Now, we define our statistical model and run the analysis. Since we want to test for association with disease status, we use a logistic regression approach. The term is necessary as placeholder for the tested HLA alleles and needs to be included. It will become handy when for example defining more complex interaction models.

In the runMiDAS function, we then refer to this model, choose our analysis type and define a inheritance model. Here we use dominant model, meaning that individuals will be defined as non-carriers (0) vs. carriers (1) for a given allele. Alternatively, it is also possible to choose recessive (0 = non-carrier or heterozygous carrier, 1 = homozygous carrier), overdominant (assuming heterozygous (dis)advantage: 0 = non-carrier or homozygous carrier, 1 = heterozygous carrier), or additive (N of alleles) inheritance models. Moreover, we define allele inclusion criteria, such that we only consider alleles frequencies above 2% and below 98%. We also apply the Bonferroni method to not only get nominal P values, but also such adjusted for multiple testing. For alternative multiple testing correction methods, as well as the option to choose a custom number of tests, please refer to the documentation. exponentiate = TRUE means that the effect estimate will already be shown as odds ratio, since we use a logistic regression model.

HLA_model <- glm(disease ~ term, data = HLA, family = binomial())
HLA_results <- runMiDAS(
  object = HLA_model, 
  experiment = "hla_alleles", 
  inheritance_model = "dominant",
  lower_frequency_cutoff = 0.02, 
  upper_frequency_cutoff = 0.98, 
  correction = "bonferroni", 
  exponentiate = TRUE
)

kableResults(HLA_results)

Three HLA alleles show significant association with the disease after multiple testing adjustment. Due to the complex linkage disequilibrium structure in the MHC, it is likely that associations are not statistically independent. The two alleles HLA-DRB1*15:01 and HLA-DQB1*06:02 are a common class II haplotype. We can therefore test if there are associations that are statistically independent of our top-associated allele, by setting the conditional flag to TRUE. MiDAS will now perform stepwise conditional testing until the top associated allele does not reach a defined significance threshold (here th = 0.05, based on adjusted p value).

HLA_results_cond <- runMiDAS(
  object = HLA_model, 
  experiment = "hla_alleles", 
  inheritance_model = "dominant", 
  conditional = TRUE,
  lower_frequency_cutoff = 0.02, 
  upper_frequency_cutoff = 0.98, 
  correction = "bonferroni", 
  exponentiate = TRUE
)

kableResults(HLA_results_cond, scroll_box_height = "200px")

The results for conditional testing are displayed in a way that for each step the top associated allele is shown, along with a list of alleles conditioned on.

As we can see, HLA-DRB1*15:01 was not independently associated with the disease when correcting for our top-associated allele HLA-DQB1*06:02. However, HLA-B*57:01 can be considered an independent association signal.

Can we find evidence for a role of HLA variation related to NK cell interactions?

Let’s assume outcome differences in our disease of interest have been shown to be related to NK cell biology. HLA class I alleles can be grouped according to how they interact with KIR, expressed on NK cells. We here now prepare a new MiDAS object, grouping HLA alleles into categories defining their interactions with KIRs (Bw4 vs. Bw6, C1 vs. C2). Then we test these variables for association with disease outcome:

NKlig <- prepareMiDAS(
  hla_calls = dat_HLA,
  colData = dat_pheno,
  experiment = c("hla_alleles", "hla_NK_ligands")
)

NKlig_model <- glm(disease ~ term, data = NKlig, family = binomial())
NKlig_results <- runMiDAS(
  NKlig_model,
  experiment = "hla_NK_ligands",
  inheritance_model = "dominant",
  correction = "bonferroni",
  exponentiate = TRUE
)

kableResults(NKlig_results)

We find that HLA-Bw6 and HLA-Bw4 carrier status are associated with decreased and increased disease risk, respectively. Of course, this is interesting enough for us to invest in some KIR typing efforts.

Do we see association on the level of KIR genes, and when considering defined HLA-KIR interactions?

Now that we have performed KIR genotyping, or e.g. inferred KIR types from available whole-genome sequencing data, we can import this information, and check the gene frequencies. In our example, we could successfully infer KIR gene presence status for 935 out of the 1,000 individuals in our data set.

dat_KIR <- readKirCalls(
  file = system.file("extdata", "MiDAS_tut_KIR.txt", package = "midasHLA")
)