--- title: "HIBAG -- an R Package for HLA Genotype Imputation with Attribute Bagging" author: "Xiuwen Zheng" date: "Jan 3, 2016" output: html_document: theme: default highlight: tango toc: yes pdf_document: toc: yes toc_depth: 3 bibliography: HIBAG_Ref.bib vignette: > %\VignetteIndexEntry{HIBAG vignette html} %\VignetteKeywords{HLA, MHC, Imputation, SNP, GWAS} %\VignetteEngine{knitr::rmarkdown} --- # Overview The human leukocyte antigen (HLA) system, located in the major histocompatibility complex (MHC) on chromosome 6p21.3, is highly polymorphic. This region has been shown to be important in human disease, adverse drug reactions and organ transplantation [@Shiina09]. HLA genes play a role in the immune system and autoimmunity as they are central to the presentation of antigens for recognition by T cells. Since they have to provide defense against a great diversity of environmental microbes, HLA genes must be able to present a wide range of peptides. Evolutionary pressure at these loci have given rise to a great deal of functional diversity. For example, the *HLA--B* locus has 1,898 four-digit alleles listed in the April 2012 release of the IMGT-HLA Database [@Robinson13] ([http://www.ebi.ac.uk/imgt/hla/](http://www.ebi.ac.uk/imgt/hla/)). Classical HLA genotyping methodologies have been predominantly developed for tissue typing purposes, with sequence based typing (SBT) approaches currently considered the gold standard. While there is widespread availability of vendors offering HLA genotyping services, the complexities involved in performing this to the standard required for diagnostic purposes make using a SBT approach time-consuming and cost-prohibitive for most research studies wishing to look in detail at the involvement of classical HLA genes in disease. Here we introduce a new prediction method for **H**LA **I**mputation using attribute **BAG**ging, HIBAG, that is highly accurate, computationally tractable, and can be used with published parameter estimates, eliminating the need to access large training samples [@Zheng13]. It relies on a training set with known HLA and SNP genotypes, and combines the concepts of attribute bagging with haplotype inference from unphased SNPs and HLA types. Attribute bagging is a technique for improving the accuracy and stability of classifier ensembles using bootstrap aggregating and random variable selection [@Breiman96; @Breiman01; @Bryll03]. In this case, individual classifiers are created which utilize a subset of SNPs to predict HLA types and haplotype frequencies estimated from a training data set of SNPs and HLA types. Each of the classifiers employs a variable selection algorithm with a random component to select a subset of the SNPs. HLA type predictions are determined by maximizing the average posterior probabilities from all classifiers. # Features * HIBAG can be used by researchers with published parameter estimates ([http://www.biostat.washington.edu/~bsweir/HIBAG/](http://www.biostat.washington.edu/~bsweir/HIBAG/)) instead of requiring access to large training sample datasets. * A typical HIBAG parameter file contains only haplotype frequencies at different SNP subsets rather than individual training genotypes. * SNPs within the xMHC region (chromosome 6) are used for imputation. * HIBAG employs unphased genotypes of unrelated individuals as a training set. * HIBAG supports parallel computing with R. # Examples ```{r echo=FALSE} options(width=100) ``` ```{r} library(HIBAG) ``` ## Pre-fit HIBAG Models for HLA Imputation ```{r fig.width=10, fig.height=5, fig.align='center'} # load the published parameter estimates from European ancestry # e.g., filename <- "HumanOmniExpressExome-European-HLA4-hg19.RData" # here, we use example data in the package filename <- system.file("extdata", "ModelList.RData", package="HIBAG") model.list <- get(load(filename)) # HLA imputation at HLA-A hla.id <- "A" model <- hlaModelFromObj(model.list[[hla.id]]) summary(model) # SNPs in the model head(model$snp.id) head(model$snp.position) ``` ```{r eval=FALSE} ######################################################################### # Import your PLINK BED file # yourgeno <- hlaBED2Geno(bed.fn=".bed", fam.fn=".fam", bim.fn=".bim") summary(yourgeno) # best-guess genotypes and all posterior probabilities pred.guess <- predict(model, yourgeno, type="response+prob") summary(pred.guess) pred.guess$value pred.guess$postprob ``` ## Build a HIBAG Model for HLA Genotype Imputation ```{r eval=FALSE} library(HIBAG) # load HLA types and SNP genotypes in the package data(HLA_Type_Table, package="HIBAG") data(HapMap_CEU_Geno, package="HIBAG") # a list of HLA genotypes # e.g., train.HLA <- hlaAllele(sample.id, H1=c("01:02", "05:01", ...), # H2=c("02:01", "03:01", ...), locus="A") # the HLA type of the first individual is 01:02/02:01, # the second is 05:01/03:01, ... hla.id <- "A" train.HLA <- hlaAllele(HLA_Type_Table$sample.id, H1 = HLA_Type_Table[, paste(hla.id, ".1", sep="")], H2 = HLA_Type_Table[, paste(hla.id, ".2", sep="")], locus=hla.id, assembly="hg19") # selected SNPs: # the flanking region of 500kb on each side, # or an appropriate flanking size without sacrificing predictive accuracy region <- 500 # kb snpid <- hlaFlankingSNP(HapMap_CEU_Geno$snp.id, HapMap_CEU_Geno$snp.position, hla.id, region*1000, assembly="hg19") length(snpid) # training genotypes # import your PLINK BED file, # e.g., train.geno <- hlaBED2Geno(bed.fn=".bed", fam.fn=".fam", bim.fn=".bim") # or, train.geno <- hlaGenoSubset(HapMap_CEU_Geno, snp.sel = match(snpid, HapMap_CEU_Geno$snp.id)) summary(train.geno) # train a HIBAG model set.seed(100) model <- hlaAttrBagging(train.HLA, train.geno, nclassifier=100) ``` ```{r echo=FALSE} mobj <- get(load(system.file("extdata", "ModelList.RData", package="HIBAG"))) model <- hlaModelFromObj(mobj$A) ``` ```{r} summary(model) ``` ## Build and Predict in Parallel ```{r eval=FALSE} library(HIBAG) library(parallel) # load HLA types and SNP genotypes in the package data(HLA_Type_Table, package="HIBAG") data(HapMap_CEU_Geno, package="HIBAG") # a list of HLA genotypes # e.g., train.HLA <- hlaAllele(sample.id, H1=c("01:02", "05:01", ...), # H2=c("02:01", "03:01", ...), locus="A") # the HLA type of the first individual is 01:02/02:01, # the second is 05:01/03:01, ... hla.id <- "A" train.HLA <- hlaAllele(HLA_Type_Table$sample.id, H1 = HLA_Type_Table[, paste(hla.id, ".1", sep="")], H2 = HLA_Type_Table[, paste(hla.id, ".2", sep="")], locus=hla.id, assembly="hg19") # selected SNPs: # the flanking region of 500kb on each side, # or an appropriate flanking size without sacrificing predictive accuracy region <- 500 # kb snpid <- hlaFlankingSNP(HapMap_CEU_Geno$snp.id, HapMap_CEU_Geno$snp.position, hla.id, region*1000, assembly="hg19") length(snpid) # training genotypes # import your PLINK BED file, # e.g., train.geno <- hlaBED2Geno(bed.fn=".bed", fam.fn=".fam", bim.fn=".bim") # or, train.geno <- hlaGenoSubset(HapMap_CEU_Geno, snp.sel = match(snpid, HapMap_CEU_Geno$snp.id)) summary(train.geno) # Create an environment with an appropriate cluster size, # e.g., 2 -- # of (CPU) cores cl <- makeCluster(2) # Building ... set.seed(1000) hlaParallelAttrBagging(cl, train.HLA, train.geno, nclassifier=100, auto.save="AutoSaveModel.RData", stop.cluster=TRUE) model.obj <- get(load("AutoSaveModel.RData")) model <- hlaModelFromObj(model.obj) summary(model) ``` ```R # best-guess genotypes and all posterior probabilities pred.guess <- predict(model, yourgeno, type="response+prob", cl=cl) ``` ## Evaluate Overall Accuracy, Sensitivity, Specificity, etc The function `hlaReport()` can be used to automatically generate a tex or HTML report when a validation dataset is available. ```{r} library(HIBAG) # load HLA types and SNP genotypes in the package data(HLA_Type_Table, package="HIBAG") data(HapMap_CEU_Geno, package="HIBAG") # make a list of HLA types hla.id <- "A" hla <- hlaAllele(HLA_Type_Table$sample.id, H1 = HLA_Type_Table[, paste(hla.id, ".1", sep="")], H2 = HLA_Type_Table[, paste(hla.id, ".2", sep="")], locus=hla.id, assembly="hg19") # divide HLA types randomly set.seed(100) hlatab <- hlaSplitAllele(hla, train.prop=0.5) names(hlatab) summary(hlatab$training) summary(hlatab$validation) # SNP predictors within the flanking region on each side region <- 500 # kb snpid <- hlaFlankingSNP(HapMap_CEU_Geno$snp.id, HapMap_CEU_Geno$snp.position, hla.id, region*1000, assembly="hg19") length(snpid) # training and validation genotypes train.geno <- hlaGenoSubset(HapMap_CEU_Geno, snp.sel = match(snpid, HapMap_CEU_Geno$snp.id), samp.sel = match(hlatab$training$value$sample.id, HapMap_CEU_Geno$sample.id)) summary(train.geno) test.geno <- hlaGenoSubset(HapMap_CEU_Geno, samp.sel=match( hlatab$validation$value$sample.id, HapMap_CEU_Geno$sample.id)) ``` ```{r eval=FALSE} # train a HIBAG model set.seed(100) model <- hlaAttrBagging(hlatab$training, train.geno, nclassifier=100) ``` ```{r echo=FALSE} mobj <- get(load(system.file("extdata", "OutOfBag.RData", package="HIBAG"))) model <- hlaModelFromObj(mobj) ``` ```{r} summary(model) # validation pred <- predict(model, test.geno) summary(pred) # compare comp <- hlaCompareAllele(hlatab$validation, pred, allele.limit=model, call.threshold=0) comp$overall ``` ### Evaluation in Figures ```{r fig.align="center", fig.height=3.3, fig.width=5} # violin plot hlaReportPlot(pred, model=model, fig="matching") ``` Matching proportion is a measure or proportion describing how the SNP profile matches the SNP haplotypes observed in the training set, i.e., the likelihood of SNP profile in a random-mating population consisting of training haplotypes. It is not directly related to confidence score, but a very low value of matching indicates that it is underrepresented in the training set. ```{r fig.align="center", fig.height=3.3, fig.width=5} hlaReportPlot(pred, hlatab$validation, fig="call.rate") hlaReportPlot(pred, hlatab$validation, fig="call.threshold") ``` ### Report in Text Output to plain text format: ```{r} # report overall accuracy, per-allele sensitivity, specificity, etc hlaReport(comp, type="txt") ``` ### Report in Markdown Output to a markdown file: ```{r results='asis'} # report overall accuracy, per-allele sensitivity, specificity, etc hlaReport(comp, type="markdown") ``` ### Report in LaTeX Output to a tex file, and please add `\usepackage{longtable}` to your tex file: ```{r} # report overall accuracy, per-allele sensitivity, specificity, etc hlaReport(comp, type="tex", header=FALSE) ``` ## Release HIBAG Models without Confidential Information ```{r eval=FALSE} library(HIBAG) # make a list of HLA types hla.id <- "DQA1" hla <- hlaAllele(HLA_Type_Table$sample.id, H1 = HLA_Type_Table[, paste(hla.id, ".1", sep="")], H2 = HLA_Type_Table[, paste(hla.id, ".2", sep="")], locus = hla.id, assembly = "hg19") # training genotypes region <- 500 # kb snpid <- hlaFlankingSNP(HapMap_CEU_Geno$snp.id, HapMap_CEU_Geno$snp.position, hla.id, region*1000, assembly="hg19") train.geno <- hlaGenoSubset(HapMap_CEU_Geno, snp.sel = match(snpid, HapMap_CEU_Geno$snp.id), samp.sel = match(hla$value$sample.id, HapMap_CEU_Geno$sample.id)) set.seed(1000) model <- hlaAttrBagging(hla, train.geno, nclassifier=100) summary(model) # remove unused SNPs and sample IDs from the model mobj <- hlaPublish(model, platform = "Illumina 1M Duo", information = "Training set -- HapMap Phase II", warning = NULL, rm.unused.snp=TRUE, anonymize=TRUE) save(mobj, file="Your_HIBAG_Model.RData") ``` ## Release a Collection of HIBAG Models ```{r eval=FALSE} # assume the HIBAG models are stored in R objects: mobj.A, mobj.B, ... ModelList <- list() ModelList[["A"]] <- mobj.A ModelList[["B"]] <- mobj.B ... # save to an R data file save(ModelList, file="HIBAG_Model_List.RData") ``` ## Association Tests In the association tests of HLA allele, four models (dominant, additive, recessive and genotype) are allowed, and linear/logistic regressions can be conducted according to the dependent variable. | Model | Description (given a specific HLA allele h) | |:----------|:--------------------------------------------| | dominant | [-/-] vs. [-/h,h/h] (0 vs. 1 in design matrix) | | additive | [-] vs. [h] in Chi-squared and Fisher's exact test, the allele dosage in regressions (0: -/-, 1: -/h, 2: h/h in design matrix) | | recessive | [-/-,-/h] vs. [h/h] (0 vs. 1 in design matrix) | | genotype | three categories: [-/-], [-/h], [h/h] | ```{r echo=FALSE} fn <- system.file("doc", "case_control.txt", package="HIBAG") ``` ```R # prepare data fn <- "case_control.txt" ``` The example text file: [case_control.txt](case_control.txt). ```{r} dat <- read.table(fn, header=TRUE, stringsAsFactors=FALSE) head(dat) # make an object for hlaAssocTest hla <- hlaAllele(dat$sample.id, H1=dat$A, H2=dat$A.1, locus="A", assembly="hg19", prob=dat$prob) summary(hla) ``` Or the best-guess HLA genotypes from `predict()` or `hlaPredict()`: ```R hla <- predict(model, yourgeno) ``` ### Allelic Association `h` in the formula is denoted for HLA genotypes. Pearson's Chi-squared and Fisher's exact test will be performed if the dependent variable is categorial, while two-sample t test or ANOVA F test will be used for continuous dependent variables. ```{r} hlaAssocTest(hla, disease ~ h, data=dat) # 95% confidence interval (h.2.5%, h.97.5%) # show details print(hlaAssocTest(hla, disease ~ h, data=dat, verbose=FALSE)) hlaAssocTest(hla, disease ~ h, data=dat, prob.threshold=0.5) # regression with a threshold hlaAssocTest(hla, disease ~ h, data=dat, showOR=TRUE) # report odd ratio instead of log odd ratio hlaAssocTest(hla, disease ~ h + pc1, data=dat) # confounding variable pc1 hlaAssocTest(hla, disease ~ h, data=dat, model="additive") # use an additive model hlaAssocTest(hla, trait ~ h, data=dat) # continuous outcome ``` ### Amino Acid Association We convert [P-coded](http://hla.alleles.org/alleles/p_groups.html) alleles to amino acid sequences. `hlaConvSequence(..., code="P.code.merge")` returns the protein sequence in the 'antigen binding domains' (exons 2 and 3 for HLA Class I genes, exon 2 for HLA Class II genes). ```{r} aa <- hlaConvSequence(hla, code="P.code.merge") ``` 1. the sequence is displayed as a hyphen (-) where it is identical to the reference. 2. an insertion or deletion is represented by a period (.). 3. an unknown or ambiguous position in the alignment is represented by an asterisk (*). 4. a capital X is used for the 'stop' codons in protein alignments. ```{r} head(c(aa$value$allele1, aa$value$allele2)) # show cross tabulation at each amino acid position summary(aa) # association tests hlaAssocTest(aa, disease ~ h, data=dat, model="dominant") # try dominant models hlaAssocTest(aa, disease ~ h, data=dat, model="dominant", prob.threshold=0.5) # try dominant models hlaAssocTest(aa, disease ~ h, data=dat, model="recessive") # try recessive models ``` # Resources * Allele Frequency Net Database (AFND): [http://www.allelefrequencies.net](http://www.allelefrequencies.net). * IMGT/HLA Database: [http://www.ebi.ac.uk/imgt/hla](http://www.ebi.ac.uk/imgt/hla). * HLA Nomenclature: [G Codes](http://hla.alleles.org/alleles/g_groups.html) and [P Codes](http://hla.alleles.org/alleles/p_groups.html) for reporting of ambiguous allele typings. # Session Info ```{r} sessionInfo() ```