The UPDhmm’s User Guide

list(name = “Marta Sevilla Porras”, affiliation = c(“Universitat Pompeu Fabra (UPF)”, “Centro de Investigación Biomédica en Red (CIBERER)”), email = “marta.sevilla@upf.edu”) list(name = “Carlos Ruiz Arenas”, affiliation = “Universidad de Navarra (UNAV)”, email = “cruizarenas@unav.es”)

2024-05-01

Introduction

Background

Uniparental disomy (UPD) is a rare genetic condition where an individual inherits both copies of a chromosome from one parent, as opposed to the typical inheritance of one copy from each parent. The extent of its contribution as a causative factor in rare genetic diseases remains largely unexplored. UPDs can lead to disease either by inheriting a carrier pathogenic mutation as homozygous from a carrier parent or by causing errors in genetic imprinting. Nevertheless, there are currently no standardized methods available for the detection and characterization of these events.

We have developed UPDhmm R/Bioconductor package. The package provides a tool method to detect, classify and stablish the location of uniparental disomy events.

Method overview

UPDhmm relies on a Hidden Markov Model (HMM) to identify regions with UPD.
In our HMM model, observations are the combination of the genotypes from the father, mother and proband for every genomic variant in the input data. The hidden states of the model represent five inheritance patterns: normal (mendelian inheritance), maternal isodisomy, paternal isodisomy, maternal heterodisomy and paternal heterodisomy. Emission probabilities were defined based on the inheritance patterns. Viterbi algorithm was employed to infer the most likely combination of hidden states underlying a sequence of observations, thus defining the most likely inheritance pattern for each genetic variant. UPDhmm reports segments in the genome having an inheritance pattern different than normal, and thus, likely harbouring a UPD.”

Set-up the packages

Loading libraries

library(UPDhmm)
library(dplyr)

Quick start

The input of the package is a multisample vcf file with the following requirements:

Datasets

The UPDhmm package includes one example dataset, adapted from GIB This dataset serves as a practical illustration and can be utilized for testing and familiarizing users with the functionality of the package.

Preprocessing

After reading the VCF file, the vcfCheck() function is employed for preprocessing the input data. This function facilitates reading the VCF in the suitable format for the UPDhmm package.

library(UPDhmm)
library(VariantAnnotation)
file <- system.file(package = "UPDhmm", "extdata", "test_het_mat.vcf.gz")
vcf <- readVcf(file)
processedVcf <- vcfCheck(
    vcf,
    proband = "NA19675", mother = "NA19678", father = "NA19679"
)

Uniparental disomy detection

The principal function of UPDhmm package, calculateEvents(), is the central function for identifying Uniparental Disomy (UPD) events. It takes the output from the previous vcfCheck() function and systematically analyzes genomic data, splitting the VCF into chromosomes and applying the Viterbi algorithm.

calculateEvents() function encompasses a serie of subprocesses for identifying Uniparental disomies (UPDs):

(1) Split vcf into chromosomes (2) Apply Viterbi algorithm to every chromosome (3) Transform the inferred hidden states to coordinates data.frame (4) Create blocks of contiguous variants with the same state (5) Calculate the statistics (log-likelihood and p-value).

updEvents <- calculateEvents(largeCollapsedVcf = processedVcf)
head(updEvents)
  seqnames     start       end   group n_snps log_likelihood      p_value
1        3 100374740 197574936 het_mat     47       75.21933 4.212224e-18
2        6  32489853  32490000 iso_mat      6       38.63453 5.110683e-10
3        6  32609207  32609341 iso_mat     11       95.23543 1.690378e-22
4       11  55322539  55587117 iso_mat      7       40.02082 2.512703e-10
5       14 105408955 105945891 het_fat     30       20.84986 4.967274e-06
6       15  22382655  23000272 iso_mat     12       48.33859 3.586255e-12
  n_mendelian_error
1                 2
2                 5
3                 7
4                 3
5                 2
6                 3

Results description

The calculateEvents function returns a data.frame containing all detected events in the provided trio. If no events are found, the function will return an empty data.frame.

Column name Description
seqnames Chromosome
start Start position of the block
end End position of the block
n_snps Number of variants within the event
group Predicted state
log_likelihood log likelihood ratio
p_value p-value
n_mendelian_error Count of Mendelian errors (genotypic combinations that are inconsistent with Mendelian inheritance principles)

Results Visualization

To visualize the results, the karyoploteR package can be used. Here, a custom function is provided to facilitate easy implementation with the output format results. This function enables the visualization of the detected blocks by the calculateEvents function.

In this example, we can observe how the package detects the simulated event on chromosome 3, as well as the specific type of event. In the plot, the autosomes chromosomes are represented, and we can see that the entire q arm of chromosome 3 is colored. With the legend, we can identify that the event is a maternal heterodisomy.

library(karyoploteR)
library(regioneR)
plotUPDKp <- function(updEvents) {
    updEvents$seqnames <- paste0("chr", updEvents$seqnames)
    het_fat <- toGRanges(subset(
        updEvents,
        group == "het_fat"
    )[, c("seqnames", "start", "end")])
    het_mat <- toGRanges(subset(
        updEvents,
        group == "het_mat"
    )[, c("seqnames", "start", "end")])
    iso_fat <- toGRanges(subset(
        updEvents,
        group == "iso_fat"
    )[, c("seqnames", "start", "end")])
    iso_mat <- toGRanges(subset(
        updEvents,
        group == "iso_mat"
    )[, c("seqnames", "start", "end")])

    kp <- plotKaryotype(genome = "hg19")
    kpPlotRegions(kp, het_fat, col = "#AAF593")
    kpPlotRegions(kp, het_mat, col = "#FFB6C1")
    kpPlotRegions(kp, iso_fat, col = "#A6E5FC")
    kpPlotRegions(kp, iso_mat, col = "#E9B864")

    colors <- c("#AAF593", "#FFB6C1", "#A6E5FC", "#E9B864")
    legend("topright",
        legend = c("Het_Fat", "Het_Mat", "Iso_Fat", "Iso_Mat"),
        fill = colors
    )
}

plotUPDKp(updEvents)

plot of chunk unnamed-chunk-1

Session Info

sessionInfo()
R version 4.4.0 beta (2024-04-15 r86425)
Platform: x86_64-pc-linux-gnu
Running under: Ubuntu 22.04.4 LTS

Matrix products: default
BLAS:   /home/biocbuild/bbs-3.19-bioc/R/lib/libRblas.so 
LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

time zone: America/New_York
tzcode source: system (glibc)

attached base packages:
[1] stats4    stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] karyoploteR_1.30.0          regioneR_1.36.0            
 [3] VariantAnnotation_1.50.0    Rsamtools_2.20.0           
 [5] Biostrings_2.72.0           XVector_0.44.0             
 [7] SummarizedExperiment_1.34.0 Biobase_2.64.0             
 [9] GenomicRanges_1.56.0        GenomeInfoDb_1.40.0        
[11] IRanges_2.38.0              S4Vectors_0.42.0           
[13] MatrixGenerics_1.16.0       matrixStats_1.3.0          
[15] BiocGenerics_0.50.0         dplyr_1.1.4                
[17] UPDhmm_1.0.0               

loaded via a namespace (and not attached):
 [1] DBI_1.2.2                bitops_1.0-7             gridExtra_2.3           
 [4] rlang_1.1.3              magrittr_2.0.3           biovizBase_1.52.0       
 [7] compiler_4.4.0           RSQLite_2.3.6            GenomicFeatures_1.56.0  
[10] png_0.1-8                vctrs_0.6.5              ProtGenerics_1.36.0     
[13] stringr_1.5.1            pkgconfig_2.0.3          crayon_1.5.2            
[16] fastmap_1.1.1            backports_1.4.1          utf8_1.2.4              
[19] rmarkdown_2.26           UCSC.utils_1.0.0         bit_4.0.5               
[22] xfun_0.43                zlibbioc_1.50.0          cachem_1.0.8            
[25] jsonlite_1.8.8           blob_1.2.4               highr_0.10              
[28] DelayedArray_0.30.0      BiocParallel_1.38.0      parallel_4.4.0          
[31] cluster_2.1.6            R6_2.5.1                 stringi_1.8.3           
[34] RColorBrewer_1.1-3       bezier_1.1.2             rtracklayer_1.64.0      
[37] rpart_4.1.23             Rcpp_1.0.12              knitr_1.46              
[40] base64enc_0.1-3          Matrix_1.7-0             nnet_7.3-19             
[43] tidyselect_1.2.1         rstudioapi_0.16.0        dichromat_2.0-0.1       
[46] abind_1.4-5              yaml_2.3.8               codetools_0.2-20        
[49] curl_5.2.1               lattice_0.22-6           tibble_3.2.1            
[52] KEGGREST_1.44.0          evaluate_0.23            foreign_0.8-86          
[55] pillar_1.9.0             BiocManager_1.30.22      checkmate_2.3.1         
[58] generics_0.1.3           RCurl_1.98-1.14          ensembldb_2.28.0        
[61] ggplot2_3.5.1            munsell_0.5.1            scales_1.3.0            
[64] BiocStyle_2.32.0         glue_1.7.0               lazyeval_0.2.2          
[67] Hmisc_5.1-2              tools_4.4.0              BiocIO_1.14.0           
[70] data.table_1.15.4        BSgenome_1.72.0          GenomicAlignments_1.40.0
[73] XML_3.99-0.16.1          grid_4.4.0               AnnotationDbi_1.66.0    
[76] colorspace_2.1-0         GenomeInfoDbData_1.2.12  htmlTable_2.4.2         
[79] restfulr_0.0.15          Formula_1.2-5            cli_3.6.2               
[82] HMM_1.0.1                fansi_1.0.6              S4Arrays_1.4.0          
[85] AnnotationFilter_1.28.0  gtable_0.3.5             digest_0.6.35           
[88] SparseArray_1.4.0        rjson_0.2.21             htmlwidgets_1.6.4       
[91] memoise_2.0.1            htmltools_0.5.8.1        lifecycle_1.0.4         
[94] httr_1.4.7               bit64_4.0.5              bamsignals_1.36.0       

References