## Load Relevant Packages
library(gdsfmt)
library(SeqArray)
library(SeqVarTools)
library(dplyr)
library(STAAR)
library(TxDb.Hsapiens.UCSC.hg38.knownGene)
library(Matrix)
library(SCANG)
library(STAARpipeline)
library(glmnet)
library(caret)
library(RISCA)
library(RICE)
RICE performs risk prediction for complex traits and diseases by incorporating both common (RICE-CV) and rare variants (RICE-RV). RICE framework consists of three main steps: (1) RICE-CV combines multiple existing common variant PRSs using ensemble learning to generate a single robust PRS; (2) RICE-RV identifies significant rare variants sets conditioned on RICE-CV, models them using statistical methods, and combines the results through ensemble learning to produce a rare variant PRS; (3) The PRS from RICE-CV and RICE-RV are then evaluated together in a regression model, adjusting for covariates. RICE requires three independent datasets: (1) a training dataset to generate GWAS summary statistics, rare variant association p-values, and train model for both common and rare variant models; (2) a tuning dataset to optimize model parameters and train the ensemble learning step; (3) a validation dataset to assess the final prediction performance.
This vignette illustrates an example of constructing a rare variant PRS using chromosome 22 of 1000 Genomes data and provides details on the ensemble regression algorithm used both RICE-RV. Data included in this vignette is located at https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/Y7AXYA and can be either directly downloaded from this link or using wget.
Note, this vignette shows how to create a rare variant PRS using burden scores extracted from an AGDS file. However, if a user has burden scores from another source, a rare variant PRS can still be constructed following the RICE-RV structure. If this is the case, start at the Penalized Regression section below.
RICE-RV performs rare variant association analysis to identify significant rare variant sets using the training dataset with STAARpipeline. With these significant rare variant sets, burden scores are constructed for each set. Then penalized regression models are trained using the burden scores and predicted onto the tuning and validation datasets. Ensemble regression optimizes the predictions to form a single PRS. Lastly, this PRS is projected onto the validation set to assess the final prediction performance.
The first step of RICE-RV involves STAARpipeline. RICE-RV can be seen as an extension of the STAARpipeline, as it relies on the results of STAARpipeline and operates with identical data structure. Therefore, before attempting to implement RICE-RV, we suggest reading through the STAARpipeline tutorial (https://github.com/xihaoli/STAARpipeline-Tutorial). This vignette will only provide brief summaries of conducting rare variant analysis with STAARpipeline.
The first step of both RICE-RV is creating annotated genomic data structure files (aGDS) for each chromosome. In the manuscript, aGDS files included the train, tune, and validation sets, as we can extract each set by subsetting. These are highly compressed files that include all variants and the functional annotations of each variant. These are commonly created using VCF files. Full details of creating aGDS files can be found in the STAARpipeline tutorial (https://github.com/xihaoli/STAARpipeline-Tutorial?tab=readme-ov-file#generate-genomic-data-structure-gds-file).
The second step of both RICE-RV is creating a baseline null model to conduct rare variant association analysis with the training dataset. RICE-RV is conducted after RICE-CV therefore the baseline model for rare variant association analysis includes the common variant PRS from RICE-CV. In the manuscript, RICE-CV was constructed with CT, LDpred2, and Lassosum2 in the case of the training dataset consisting of a single ancestry and CT-SLEB, PROSPER, and JointPRS in the case of the training dataset consisting of multiple ancestries. Full details on the implementation of this can be found at … and https://github.com/jwilliams10/RareVariantPRS. However, in practice RICE-RV can be constructed conditional on any common variant PRS. Therefore, here we fit the baseline model with a simulated common variant PRS.
## Create Train/Tune/Validation Partitions
set.seed(1330)
### Get Sample IDs
genofile <- seqOpen("Chr22_1000G.gds")
sampleids <- seqGetData(genofile,"sample.id")
seqClose(genofile)
### 70% Train, 15% Tune, 15% Validation
train_set <- sample(sampleids,round(length(sampleids)*0.7))
remaining_set <- sampleids[!(sampleids %in% train_set)]
tune_set <- sample(remaining_set,round(length(remaining_set)*0.5))
validation_set <- sampleids[!(sampleids %in% c(train_set,tune_set))]
## Fit Baseline Model for Training Data
### Load Phenotype and Common Variant PRS
load("Phenotype.RData")
load("Common_Variant_PRS.RData")
### Merge Phenotype and Common Variant PRS
Phenotype_Data <- inner_join(Y,Common_Variant_PRS)
### Subset to the Training Set
Phenotype_Data_Train <- Phenotype_Data[Phenotype_Data$IDs %in% train_set,]
### Fit Null Model without kinship, specifying continuous outcome, and column name of ID column
Null_Model <- fit_nullmodel(Y~PRS, data = Phenotype_Data_Train,id = "IDs",kins = NULL,family = gaussian(link = "identity"))
#> [1] "kins is NULL, fit generalized linear model."
The third step of RICE-RV is to perform rare variant association analysis using STAARpipeline. Rare variant association analysis within STAARpipeline can be divided into two analyses; gene-centric coding analysis and gene-centric noncoding analysis. Gene-centric coding analysis conducts association analysis with rare variants located in protein-coding genes in the coding region and gene-centric noncoding analysis in the noncoding region. Whole exome sequencing data is largely comprised of rare variants located in the coding region while whole genome sequencing data has rare variants in both the coding and noncoding regions. This vignette only provides details for an example gene-centric coding analysis. Full details on how to implement gene-centric noncoding analysis and implementing both coding and noncoding analysis within a HPC is available at https://github.com/xihaoli/STAARpipeline-Tutorial?tab=readme-ov-file#step-31-gene-centric-coding-analysis.
Here we implement gene-centric coding analysis on the training dataset for a subset of chromosome 22. We implement for only a subset as this is a non-optimal way of implementing this and takes a while to run, if you have access to a HPC or cloud computing we recommend submitting this as a batch job as shown in the STAARpipeline tutorial.
### Gene-Centric Coding Analysis
## QC_label
QC_label <- "annotation/info/QC_label"
## variant_type
variant_type <- "SNV"
## geno_missing_imputation
geno_missing_imputation <- "mean"
## Annotation_dir
Annotation_dir <- "annotation/info/FunctionalAnnotation"
## Annotation channel
Annotation_name_catalog <- read.csv("Annotation_name_catalog.csv")
## Use_annotation_weights
Use_annotation_weights <- TRUE
## Annotation name
Annotation_name <- c("CADD","LINSIGHT","FATHMM.XF","aPC.EpigeneticActive","aPC.EpigeneticRepressed","aPC.EpigeneticTranscription",
"aPC.Conservation","aPC.LocalDiversity","aPC.Mappability","aPC.TF","aPC.Protein")
## Protein-Coding Genes in Chromosome 22
genes_info_chr <- genes_info[genes_info[,2]==22,]
Genes <- c("OR11H1","POTEH","CCT8L2","XKR3","GAB4","IL17RA","RANBP1","PI4KA","SERPIND1","ZNF280B","GGTLC2","PIK3IP1")
genes_info_chr_subset <- genes_info_chr[genes_info_chr$hgnc_symbol %in% Genes,]
### Gene-Centric Coding Analysis
genofile <- seqOpen("Chr22_1000G.gds")
results_coding <- c()
for(kk in 1:nrow(genes_info_chr_subset)){
print(kk)
gene_name <- genes_info_chr_subset[kk,1]
results <- Gene_Centric_Coding(chr=22,gene_name=gene_name,genofile=genofile,obj_nullmodel=Null_Model,
rare_maf_cutoff=0.01,rv_num_cutoff=2,
QC_label=QC_label,variant_type=variant_type,geno_missing_imputation=geno_missing_imputation,
Annotation_dir=Annotation_dir,Annotation_name_catalog=Annotation_name_catalog,
Use_annotation_weights=Use_annotation_weights,Annotation_name=Annotation_name,silent = TRUE)
results_coding <- append(results_coding,results)
}
seqClose(genofile)
The fourth step of is extracting the burden scores of the significant rare variant sets. In our manuscript, significance is defined as having the STAAR-Burden p-value less than \(10^{-3}\), however, an individual can construct burden scores for any set of rare variant sets. To identify significant rare variant sets, we first compile the results from the gene-centric coding analysis.
### Compile Gene-Centric Coding Analysis
coding_sig <- NULL
for(j in 1:length(results_coding)){
coding_sig <- rbind(coding_sig,unlist(results_coding[[j]][,c(1:4,61)]))
}
coding_sig <- as.data.frame(coding_sig)
colnames(coding_sig) <- c("Gene","Chr","Category","Number_SNV","STAARB")
coding_sig$Chr <- as.numeric(coding_sig$Chr)
coding_sig$Number_SNV <- as.numeric(coding_sig$Number_SNV)
coding_sig$STAARB <- as.numeric(coding_sig$STAARB)
coding_sig <- coding_sig[coding_sig$Number_SNV < 2000,]
Again this code is written for this vignette, when using a HPC or cloud computing refer to the STAARpipeline tutorial to optimally compile the results. Once compiled the structure of the results should look like this:
### Head of Significant Rare Variant Set Results
coding_sig <- coding_sig[coding_sig$STAARB < 1e-3,]
head(coding_sig)
#> Gene Chr Category Number_SNV STAARB
#> 18 RANBP1 22 missense 5 4.566458e-05
#> 23 PI4KA 22 synonymous 60 5.117859e-04
#> 25 SERPIND1 22 missense 15 6.412555e-04
#> 28 ZNF280B 22 missense 8 8.451143e-08
#> 30 GGTLC2 22 missense 26 5.101352e-06
#> 33 PIK3IP1 22 synonymous 3 5.468301e-07
Note, the results include the protein-coding gene, the chromosome, and the functional category. These are essential elements to define a rare variant set to build and extract a burden score. Below is how to compute the burden scores for all significant rare variant sets for all individuals. Burden scores are independent across individuals, therefore they can be computed for all individuals then separated into train, tune, and validation sets. Further note, that the STAAR-Burden(1,1) p-value is the only p-value, this is because this the p-value used to determine significance in the manuscript. However, the raw output of the rare variant association analysis includes much more information outlined in the STAARpipeline tutorial. With the significant rare variant sets, we can extract and build the burden scores of the significant rare variant sets.
### Significant Rare Variant Set's Burden Scores
genofile <- seqOpen("Chr22_1000G.gds")
Sig_Burdens <- NULL
### Loop over significant rare variant sets
### Using IDs = sampleids to build burden scores for all individuals.
for(i in 1:nrow(coding_sig)){
chr <- coding_sig$Chr[i]
gene_name <- coding_sig$Gene[i]
category <- coding_sig$Category[i]
a <- Burden_Scores(region = "Coding",chr = chr, gene_name = gene_name,category = category,
genofile = genofile,IDs = sampleids,rare_maf_cutoff=0.01,rv_num_cutoff=2,
QC_label=QC_label,variant_type=variant_type,geno_missing_imputation=geno_missing_imputation,
Annotation_dir=Annotation_dir,Annotation_name_catalog=Annotation_name_catalog,silent = TRUE)
### Combine them column-wise
Sig_Burdens <- cbind(Sig_Burdens,a)
}
seqClose(genofile)
### The burden scores align with the order of the sampleids
IDs_Sig_Burdens <- sampleids
The fifth step of RICE-RV trains penalized regression models with the trait as the outcome and the significant burden scores as an input. Once these are trained, we predict the outcome for the tuning and validation. Note we do not choose optimal hyper-parameters for the penalized regression models, we simply use the entire sequence of hyper-parameters.
## Penalized Regression
### Align with phenotype data
X <- data.frame(IDs = IDs_Sig_Burdens,Sig_Burdens)
Phenotype_Data <- inner_join(Phenotype_Data,X)
#> Joining with `by = join_by(IDs)`
### Extract Aligned Burdens, change structure
X_Aligned <- as.matrix(Phenotype_Data[,paste0("X",1:ncol(Sig_Burdens))])
### Subset the Burden Scores to Train, Tune, Validation
X_train <- X_Aligned[Phenotype_Data$IDs %in% train_set,]
X_tune <- X_Aligned[Phenotype_Data$IDs %in% tune_set,]
X_valid <- X_Aligned[Phenotype_Data$IDs %in% validation_set,]
### Build a null model using control covariates to train penalized regression models
train_null_model <- lm(Y~1,data = Phenotype_Data[Phenotype_Data$IDs %in% train_set,])
### Extract residuals to use as a response
Residuals_Train <- train_null_model$residuals
### Train penalized regression models on the training set
lasso_train <- glmnet(X_train,Residuals_Train,family = "gaussian",alpha = 1)
ridge_train <- glmnet(X_train,Residuals_Train,family = "gaussian",alpha = 0)
### Project onto tuning and validation sets
lasso_tune <- predict(lasso_train,as.matrix(X_tune))
ridge_tune <- predict(ridge_train,X_tune)
lasso_valid <- predict(lasso_train,X_valid)
ridge_valid <- predict(ridge_train,X_valid)
### Combine
all_tune <- cbind(lasso_tune,ridge_tune)
colnames(all_tune) <- c(paste0("lasso",1:ncol(lasso_tune)),paste0("ridge",1:ncol(ridge_tune)))
all_valid <- cbind(lasso_valid,ridge_valid)
colnames(all_valid) <- c(paste0("lasso",1:ncol(lasso_valid)),paste0("ridge",1:ncol(ridge_valid)))
The sixth and last step of RICE-RV, trains an ensemble regression algorithm on the predicted penalized regression outcomes with the tuning dataset. The ensemble regression algorithm includes both Lasso and ridge regression as base learners.
## Ensemble Regression
### Build a null model using control covariates to train the Ensembler Regression Algorithm
tune_null_model <- lm(Y~1,data = Phenotype_Data[Phenotype_Data$IDs %in% tune_set,])
y_tune <- tune_null_model$residuals
Results <- Ensemble_Function(PRSs = all_tune,Y = y_tune,family = "continuous")
## Set coefficients of rank deficient columns to 0
Results$Coefficients[is.na(Results$Coefficients)] <- 0
### Save final PRS for tuning and validation sets
BestRareVariantPRS_Tune <- data.frame(IDs = tune_set,PRS = as.matrix(all_tune[,names(Results$Coefficients)[-1]]) %*% matrix(Results$Coefficients[-1],ncol = 1))
BestRareVariantPRS_Valid <- data.frame(IDs = validation_set,PRS = as.matrix(all_valid[,names(Results$Coefficients)[-1]]) %*% matrix(Results$Coefficients[-1],ncol = 1))
We can evaluate the final prediction performance on the validation dataset.
valid_null_model <- lm(Y~1,data = Phenotype_Data[Phenotype_Data$IDs %in% validation_set,])
y_valid <- valid_null_model$residuals
### R2 of RICE-RV on Validation Data
summary(lm(y~x,data = data.frame(y = y_valid,x = BestRareVariantPRS_Valid$PRS)))$r.squared
#> [1] 0.04935384
### Beta of PRS per SD of RICE-RV on Validation Data
coef(lm(y~x - 1,data = data.frame(y = scale(y_valid),x = scale(BestRareVariantPRS_Valid$PRS))))
#> x
#> 0.2221572
Commonly, there is curiosity in knowing the relationship between the variants and the final PRS. As all of the ensemble learning is done using linear algorithms, the final predicted PRS is then a linear projection of the significant burden scores. Using the tuning dataset, we can estimate the final coefficients of the significant rare variant sets with the following code.
## Estimate Final Coefficients
Linear_Projection <- data.frame(y = BestRareVariantPRS_Tune$PRS,X_tune)
### The tuning PRS is a perfect linear projection of the significant burden scores
summary(lm(y~.,data =Linear_Projection))$r.squared
#> [1] 1
### Extract the coefficients
coding_sig$Final_Coefficients <- unname(coef(lm(y~.,data = Linear_Projection))[-1])
head(coding_sig)
#> Gene Chr Category Number_SNV STAARB Final_Coefficients
#> 18 RANBP1 22 missense 5 4.566458e-05 1.08275160
#> 23 PI4KA 22 synonymous 60 5.117859e-04 0.12985495
#> 25 SERPIND1 22 missense 15 6.412555e-04 -0.06010972
#> 28 ZNF280B 22 missense 8 8.451143e-08 1.01263845
#> 30 GGTLC2 22 missense 26 5.101352e-06 -0.20726166
#> 33 PIK3IP1 22 synonymous 3 5.468301e-07 0.97698168
Here we presented how to construct a rare variant PRS using rare
variants located in coding region of chromosome 22. This example is
almost identical to the simulation study in our manuscript but is far
from the complexities in real world data analyses in whole exome and
whole genome sequencing data. To extend this to real data analyses,
common extensions are the inclusion of all chromosomes, including both
rare variants located in coding and noncoding regions, and both binary
and continuous traits. To extend to all chromosomes, perform rare
variant association analysis on all chromosomes using the STAARpipeline,
extract burden scores for each chromosome using the code shown above,
compile the significant burden scores from all chromosomes together, and
perform the identical ensemble regression algorithm above. To extend to
noncoding regions, perform a separate gene-centric noncoding rare
variant association analysis with STAARpipeline, extract significant
burden scores using the same function but changing the region argument,
combine significant coding and noncoding burden scores, and perform the
identical ensemble regression algorithm. Lastly to extend to binary
traits, fit a binary null model with STAARpipeline, perform the
respective rare variant association tests, extract the significant
burden scores, implement the penalized regression models,and set
family = "binary"
in the Ensemble_Function
.
The implementation of this algorithm for binary traits for the UK
Biobank whole exome sequencing analysis can be found at https://github.com/jwilliams10/RareVariantPRS/blob/main/WES/Binary/RareVariant_PRS/Single_RareVariant_PRS_All.R.
sessionInfo()
#> R version 4.4.1 (2024-06-14)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS 15.5
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
#>
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> time zone: America/New_York
#> tzcode source: internal
#>
#> attached base packages:
#> [1] splines stats4 stats graphics grDevices utils datasets
#> [8] methods base
#>
#> other attached packages:
#> [1] RICE_0.1.0
#> [2] RISCA_1.0.7
#> [3] mosaicData_0.20.4
#> [4] ggformula_0.12.0
#> [5] tune_1.3.0
#> [6] reticulate_1.42.0
#> [7] relsurv_2.3-2
#> [8] survival_3.7-0
#> [9] caret_7.0-1
#> [10] lattice_0.22-6
#> [11] ggplot2_3.5.2
#> [12] glmnet_4.1-9
#> [13] STAARpipeline_0.9.7
#> [14] SCANG_1.0.4
#> [15] Matrix_1.7-0
#> [16] TxDb.Hsapiens.UCSC.hg38.knownGene_3.18.0
#> [17] GenomicFeatures_1.56.0
#> [18] AnnotationDbi_1.66.0
#> [19] Biobase_2.64.0
#> [20] GenomicRanges_1.56.1
#> [21] GenomeInfoDb_1.40.1
#> [22] IRanges_2.38.1
#> [23] S4Vectors_0.42.1
#> [24] BiocGenerics_0.50.0
#> [25] STAAR_0.9.7
#> [26] dplyr_1.1.4
#> [27] SeqVarTools_1.42.0
#> [28] SeqArray_1.44.0
#> [29] gdsfmt_1.40.0
#>
#> loaded via a namespace (and not attached):
#> [1] matrixStats_1.4.1 bitops_1.0-9
#> [3] lubridate_1.9.4 DiceDesign_1.10
#> [5] httr_1.4.7 RColorBrewer_1.1-3
#> [7] doParallel_1.0.17 tools_4.4.1
#> [9] backports_1.5.0 R6_2.6.1
#> [11] mgcv_1.9-1 yardstick_1.3.2
#> [13] jomo_2.7-6 SNPRelate_1.38.0
#> [15] withr_3.0.2 quantreg_5.98
#> [17] cli_3.6.5 sandwich_3.1-1
#> [19] sass_0.4.9 nnls_1.6
#> [21] mvtnorm_1.3-1 quantsmooth_1.70.0
#> [23] ggridges_0.5.6 Rsamtools_2.20.0
#> [25] parallelly_1.45.0 labelled_2.14.1
#> [27] rstudioapi_0.16.0 RSQLite_2.3.7
#> [29] generics_0.1.4 shape_1.4.6.1
#> [31] BiocIO_1.14.0 abind_1.4-8
#> [33] lifecycle_1.0.4 multcomp_1.4-28
#> [35] yaml_2.3.10 CompQuadForm_1.4.3
#> [37] GWASTools_1.50.0 SummarizedExperiment_1.34.0
#> [39] recipes_1.3.1 SparseArray_1.4.8
#> [41] grid_4.4.1 blob_1.2.4
#> [43] crayon_1.5.3 mitml_0.4-5
#> [45] haven_2.5.4 KEGGREST_1.44.1
#> [47] pillar_1.10.2 knitr_1.48
#> [49] rjson_0.2.23 boot_1.3-31
#> [51] logistf_1.26.0 estimability_1.5.1
#> [53] future.apply_1.11.3 codetools_0.2-20
#> [55] pan_1.9 glue_1.8.0
#> [57] rsample_1.3.0 data.table_1.17.4
#> [59] vctrs_0.6.5 png_0.1-8
#> [61] gtable_0.3.6 kernlab_0.9-33
#> [63] cachem_1.1.0 gower_1.0.2
#> [65] xfun_0.47 S4Arrays_1.4.1
#> [67] prodlim_2025.04.28 timeDate_4041.110
#> [69] iterators_1.0.14 hardhat_1.4.1
#> [71] lava_1.8.1 gam_1.22-5
#> [73] statmod_1.5.0 TH.data_1.1-3
#> [75] ipred_0.9-15 nlme_3.1-166
#> [77] bit64_4.5.2 kinship2_1.9.6.1
#> [79] bslib_0.8.0 rpart_4.1.23
#> [81] DBI_1.2.3 nnet_7.3-19
#> [83] DNAcopy_1.78.0 tidyselect_1.2.1
#> [85] emmeans_1.10.4 bit_4.5.0
#> [87] compiler_4.4.1 curl_5.2.3
#> [89] GENESIS_2.34.0 mice_3.16.0
#> [91] SparseM_1.84-2 DelayedArray_0.30.1
#> [93] rtracklayer_1.64.0 mosaicCore_0.9.4.0
#> [95] scales_1.4.0 lmtest_0.9-40
#> [97] MultiSTAAR_0.9.7 quadprog_1.5-8
#> [99] stringr_1.5.1 digest_0.6.37
#> [101] minqa_1.2.8 GWASExactHW_1.2
#> [103] rmarkdown_2.28 XVector_0.44.0
#> [105] htmltools_0.5.8.1 pkgconfig_2.0.3
#> [107] lme4_1.1-35.5 lhs_1.2.0
#> [109] MatrixGenerics_1.16.0 fastmap_1.2.0
#> [111] rlang_1.1.6 UCSC.utils_1.0.0
#> [113] farver_2.1.2 jquerylib_0.1.4
#> [115] zoo_1.8-12 jsonlite_1.8.9
#> [117] BiocParallel_1.38.0 ModelMetrics_1.2.2.2
#> [119] RCurl_1.98-1.16 magrittr_2.0.3
#> [121] GenomeInfoDbData_1.2.12 GPfit_1.0-9
#> [123] Rcpp_1.0.14 furrr_0.3.1
#> [125] stringi_1.8.7 pROC_1.18.5
#> [127] zlibbioc_1.50.0 MASS_7.3-61
#> [129] plyr_1.8.9 formula.tools_1.7.1
#> [131] parallel_4.4.1 listenv_0.9.1
#> [133] forcats_1.0.0 Biostrings_2.72.1
#> [135] hms_1.1.3 cubature_2.1.4
#> [137] dials_1.4.0 reshape2_1.4.4
#> [139] parsnip_1.3.2 XML_3.99-0.17
#> [141] evaluate_1.0.0 operator.tools_1.6.3
#> [143] nloptr_2.1.1 foreach_1.5.2
#> [145] MatrixModels_0.5-3 tidyr_1.3.1
#> [147] purrr_1.0.4 future_1.49.0
#> [149] SuperLearner_2.0-29 broom_1.0.7
#> [151] xtable_1.8-4 restfulr_0.0.15
#> [153] class_7.3-22 GMMAT_1.4.2
#> [155] tibble_3.2.1 memoise_2.0.1
#> [157] GenomicAlignments_1.40.0 mosaic_1.9.1
#> [159] workflows_1.2.0 timechange_0.3.0
#> [161] globals_0.18.0