CheckFinemapSNPEnrichments

Intro

One way to explore my QTL results is to look for enrichments in QTL SNPs. There are many ways to do this. I propose something simple: Start with the finemapping results from hyprcoloc (the posterior probability for each SNP for each hyprcoloc cluster) and ask to what degree that posterior probability mass is enriched in annotations, compared to some other set of SNPs (eg sQTLs versus eQTLs). The annotations could be anything, but here I am just going to explore this method and make some initial plots using GM12878 ChromHMM annotations from ENCODE (lifted over from hg19->hg38).

library(tidyverse)
library(data.table)
library(boot)

dat <- fread("../code/QTL_SNP_Enrichment/FinemapIntersections.bed.gz", col.names=c("SNP_chrom", "SNP_start", "SNP_end", "SNP", "PosteriorPr", "Annotation_chrom", "Annotation_start", "Annotation_stop", "Annotation", "Overlap")) %>%
  separate(SNP, into=c("SNP", "HyprcolocCluster", "GeneLocus"), sep="_")

hyprcoloc_results <- read_tsv("../code/hyprcoloc/Results/ForColoc/MolColocStandard/tidy_results_OnlyColocalized.txt.gz")



hyprcoloc_results_finemap <- read_tsv("../code/hyprcoloc/Results/ForColoc/MolColocStandard/snpscores.txt.gz")

First I want to verify that I can match up the right TopSNPs by the clusterNumber. If I matching these TopSNPs correctly, the finemapping probability should be equal from both tables…

hyprcoloc_results_finemap %>%
  group_by(ColocalizedCluster, Locus) %>%
  filter(FinemapPr == max(FinemapPr)) %>%
  ungroup() %>%
  full_join(
    hyprcoloc_results %>% distinct(snp, iteration, Locus, .keep_all=T),
    by = c("Locus", "snp", "ColocalizedCluster"="iteration")
  ) %>%
  ggplot(aes(x=FinemapPr, y=TopSNPFinemapPr)) +
  geom_point()

That looks right…

Now let’s pick clusters with only splicing phenotypes, and compare to clusters with an eQTL

hyprcoloc_results %>%
  separate(phenotype_full, into=c("PC", "P"), sep=";") %>%
  pull(PC) %>% unique()

 [1] "polyA.Splicing"                      "chRNA.IR"                           
 [3] "polyA.IR"                            "MetabolicLabelled.60min.IR"         
 [5] "MetabolicLabelled.30min.IR"          "MetabolicLabelled.30min"            
 [7] "MetabolicLabelled.60min"             "Expression.Splicing"                
 [9] "H3K27AC"                             "H3K4ME3"                            
[11] "H3K4ME1"                             "chRNA.Expression_eRNA"              
[13] "polyA.IR.Subset_YRI"                 "chRNA.IER"                          
[15] "chRNA.IRjunctions"                   "ProCap"                             
[17] "polyA.Splicing.Subset_YRI"           "CTCF"                               
[19] "chRNA.Expression_lncRNA"             "Expression.Splicing.Subset_YRI"     
[21] "chRNA.Splicing"                      "MetabolicLabelled.30min.Splicing"   
[23] "MetabolicLabelled.60min.Splicing"    "chRNA.Expression.Splicing"          
[25] "H3K36ME3"                            "polyA.IER"                          
[27] "polyA.IRjunctions"                   "chRNA.Expression_cheRNA"            
[29] "chRNA.Slopes"                        "MetabolicLabelled.30min.IER"        
[31] "MetabolicLabelled.60min.IER"         "MetabolicLabelled.60min.IRjunctions"
[33] "MetabolicLabelled.30min.IRjunctions" "chRNA.Expression_snoRNA"            

SplicingClusters <- hyprcoloc_results %>%
  separate(phenotype_full, into=c("PC", "P"), sep=";") %>%
  rowwise() %>%
  mutate(PC_new = if_else(
    any(str_detect(PC, c("polyA.Splicing","chRNA.IR", "polyA.IR", "polyA.IR.Subset_YRI", "chRNA.IER", "chRNA.IRjunctions", "polyA.Splicing.Subset_YRI", "chRNA.Splicing"))),
    "Splicing",
    "NotSplicing"
  )) %>%
  ungroup() %>%
  group_by(Locus, iteration) %>%
  filter(all(PC_new == "Splicing")) %>%
  ungroup() %>%
  distinct(Locus,iteration) %>%
  unite(Locus.iteration, Locus,iteration) %>%
  mutate(ClusterGroup = "SplicingClusters")



H3K27AC_AndEqtlClusters <- hyprcoloc_results %>%
  separate(phenotype_full, into=c("PC", "P"), sep=";") %>%
  rowwise() %>%
  mutate(PC_new = if_else(
    any(str_detect(PC, c("H3K27AC","H3K4ME1", "H3K36ME3", "H3K4ME3", "Expression.Splicing", "Expression.Splicing.Subset_YRI", "chRNA.Expression.Splicing", "MetabolicLabelled.30min", "MetabolicLabelled.60min"))),
    "EnhancerOrExpression",
    "NotEnhancerOrExpression"
  )) %>%
  ungroup() %>%
  group_by(Locus, iteration) %>%
  filter(all(PC_new == "EnhancerOrExpression")) %>%
  ungroup() %>%
  distinct(Locus,iteration) %>%
  unite(Locus.iteration, Locus,iteration) %>%
  mutate(ClusterGroup = "EnhancerOrExpression")

intersect(H3K27AC_AndEqtlClusters$Locus.iteration, SplicingClusters$Locus.iteration)

character(0)

FractionPosteriorInEachAnnotation <- bind_rows(H3K27AC_AndEqtlClusters, SplicingClusters) %>%
  left_join(
    dat %>%
    mutate(Locus.iteration = paste(GeneLocus, HyprcolocCluster, sep="_")),
    by = "Locus.iteration"
  ) %>%
  group_by(Annotation, ClusterGroup) %>%
  summarise(Sum = sum(PosteriorPr)) %>%
  group_by(ClusterGroup) %>%
  mutate(FractionAnnotatedPosteriorInGroup = Sum/sum(Sum, na.rm=T)) %>%
  ungroup()

FractionPosteriorInEachAnnotation %>%
  pivot_wider(names_from = "ClusterGroup", values_from = c("FractionAnnotatedPosteriorInGroup", "Sum"), id_cols="Annotation") %>%
  mutate(Enrichment = FractionAnnotatedPosteriorInGroup_SplicingClusters/FractionAnnotatedPosteriorInGroup_EnhancerOrExpression) %>%
  mutate(Annotation = str_replace_all(Annotation, "_0$", "_UnannotatedButObserved")) %>%
  mutate(Annotation = str_replace_all(Annotation, "_1$", "_Annotated")) %>%
ggplot(aes(x=Annotation, y=Enrichment)) +
  geom_col() +
  theme_bw() +
  scale_y_continuous(trans='log2') +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  coord_flip() +
  labs(x = "Annotation", y="Finemapping posterior mass fold enrichment\n(ChromatinQTL|eQTL <-- --> sQTL)")

Let’s try to add error bars by bootstrapping. Perhaps randomly sample the QTL clusters. Let’s try 100 bootstrap replicates and make the plot twice to see how stable the confidence interval is.

set.seed(0)

dat$Annotation %>% unique()

 [1] "SpliceBranchpointRegion_0" "10_Txn_Elongation"        
 [3] "6_Weak_Enhancer"           "SpliceDonor_0"            
 [5] "2_Weak_Promoter"           "11_Weak_Txn"              
 [7] "1_Active_Promoter"         "12_Repressed"             
 [9] "4_Strong_Enhancer"         "13_Heterochrom/lo"        
[11] "SpliceAcceptor_0"          "SpliceBranchpointRegion_1"
[13] "7_Weak_Enhancer"           "9_Txn_Transition"         
[15] "SpliceAcceptor_1"          "SpliceDonor_1"            
[17] "."                         "5_Strong_Enhancer"        
[19] "8_Insulator"               "3_Poised_Promoter"        
[21] "14_Repetitive/CNV"         "15_Repetitive/CNV"        

results <- list()

for (i in 1:100){
  results[[i]] <-
  bind_rows(
  sample_frac(H3K27AC_AndEqtlClusters, replace=T),
  sample_frac(SplicingClusters, replace=T)) %>%
  left_join(
    dat %>%
    mutate(Locus.iteration = paste(GeneLocus, HyprcolocCluster, sep="_")),
    by = "Locus.iteration"
  ) %>%
  mutate(Annotation = str_replace_all(Annotation, "_0$", "_UnannotatedButObserved")) %>%
  mutate(Annotation = str_replace_all(Annotation, "_1$", "_Annotated")) %>%
  mutate(Annotation = str_replace_all(Annotation, "^\\d+?_", "ChromHMM:")) %>%
  filter(!is.na(Annotation)) %>%
  mutate(Annotation = if_else(Annotation==".", "No annotation", Annotation)) %>%
  group_by(Annotation, ClusterGroup) %>%
  summarise(Sum = sum(PosteriorPr)) %>%
  group_by(ClusterGroup) %>%
  mutate(FractionAnnotatedPosteriorInGroup = Sum/sum(Sum, na.rm=T)) %>%
  ungroup() %>%
  pivot_wider(names_from = "ClusterGroup", values_from = c("FractionAnnotatedPosteriorInGroup", "Sum"), id_cols="Annotation") %>%
  mutate(Enrichment = FractionAnnotatedPosteriorInGroup_SplicingClusters/FractionAnnotatedPosteriorInGroup_EnhancerOrExpression) %>%
  mutate(i=i)
}

bind_rows(results) %>%
  group_by(Annotation) %>%
  summarise(mean = exp(mean(log(Enrichment), na.rm=T)),
            Percentile10 = quantile(Enrichment, probs = 0.05, na.rm=T),
            Percentile90 = quantile(Enrichment, probs = 0.95, na.rm=T)) %>%
  ggplot(aes(x=Annotation, y=mean)) +
    # geom_vline(xintercept = 1) +
    geom_col() +
    geom_errorbar(aes(ymin=Percentile10, ymax=Percentile90)) +
    theme_bw() +
    scale_y_continuous(trans='log2') +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
    coord_flip() +
    labs(x = "Annotation", y="Finemapping posterior mass fold enrichment\n(ChromatinQTL|eQTL <-- --> sQTL)",
         caption = "Bootsrapped 90% confidence interval")

set.seed(100)
results <- list()

for (i in 1:100){
  results[[i]] <-
  bind_rows(
  sample_frac(H3K27AC_AndEqtlClusters, replace=T),
  sample_frac(SplicingClusters, replace=T)) %>%
  left_join(
    dat %>%
    mutate(Locus.iteration = paste(GeneLocus, HyprcolocCluster, sep="_")),
    by = "Locus.iteration"
  ) %>%
  mutate(Annotation = str_replace_all(Annotation, "_0$", "_UnannotatedButObserved")) %>%
  mutate(Annotation = str_replace_all(Annotation, "_1$", "_Annotated")) %>%
  mutate(Annotation = str_replace_all(Annotation, "^\\d+?_", "ChromHMM:")) %>%
  filter(!is.na(Annotation)) %>%
  mutate(Annotation = if_else(Annotation==".", "No annotation", Annotation)) %>%
  group_by(Annotation, ClusterGroup) %>%
  summarise(Sum = sum(PosteriorPr)) %>%
  group_by(ClusterGroup) %>%
  mutate(FractionAnnotatedPosteriorInGroup = Sum/sum(Sum, na.rm=T)) %>%
  ungroup() %>%
  pivot_wider(names_from = "ClusterGroup", values_from = c("FractionAnnotatedPosteriorInGroup", "Sum"), id_cols="Annotation") %>%
  mutate(Enrichment = FractionAnnotatedPosteriorInGroup_SplicingClusters/FractionAnnotatedPosteriorInGroup_EnhancerOrExpression) %>%
  mutate(i=i)
}

bind_rows(results) %>%
  group_by(Annotation) %>%
  summarise(mean = exp(mean(log(Enrichment), na.rm=T)),
            Percentile10 = quantile(Enrichment, probs = 0.05, na.rm=T),
            Percentile90 = quantile(Enrichment, probs = 0.95, na.rm=T)) %>%
  ggplot(aes(x=Annotation, y=mean)) +
    # geom_vline(xintercept = 1) +
    geom_col() +
    geom_errorbar(aes(ymin=Percentile10, ymax=Percentile90)) +
    theme_bw() +
    scale_y_continuous(trans='log2') +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
    coord_flip() +
    labs(x = "Annotation", y="Finemapping posterior mass fold enrichment\n(ChromatinQTL|eQTL <-- --> sQTL)",
         caption = "Bootsrapped 90% confidence interval")

Ok those two separate graphs look similar enough for me, that I think this is enough bootstrap replicates. Let’s wrap this process into a function and try making similar plots comparing different sets of QTLs…

Actually, first, before considering “relative enrichment” in comparing two sets of QTLs, I want to get a sense of the total mass under each annotation in both sets…

FractionPosteriorInEachAnnotation %>%
  ggplot(aes(x=Annotation, y=FractionAnnotatedPosteriorInGroup*100)) +
  geom_col() +
  facet_wrap(~ClusterGroup) +
  coord_flip() +
  theme_bw()

FractionPosteriorInEachAnnotation %>%
  ggplot(aes(x=ClusterGroup, y=FractionAnnotatedPosteriorInGroup*100, fill=Annotation)) +
  geom_col(position='fill') +
  theme_bw()

Ok, so perhaps only a 1/3 of the posterior mass for sQTLs lies in splice site regions.

Now, let’s go back to wrapping the enrichment plot between two sets of QTLs into a function, so I can reuse it to comapre different sets of QTLs

iterations = 100
GeneLocus_Cluster_QTLSet1 = H3K27AC_AndEqtlClusters
GeneLocus_Cluster_QTLSet2 = SplicingClusters

#Define the function
PlotDifferentialEnrichmentBetweenQTLSets <- function(dat, GeneLocus_Cluster_QTLSet1, GeneLocus_Cluster_QTLSet2, seed=0, iterations=100){
  output <- list()
  results <- list()
  for (i in 1:iterations){
    results[[i]] <-
    bind_rows(
    sample_frac(GeneLocus_Cluster_QTLSet1, replace=T),
    sample_frac(GeneLocus_Cluster_QTLSet2, replace=T),
    .id="QTL.Set") %>%
    left_join(
      dat %>%
      mutate(Locus.iteration = paste(GeneLocus, HyprcolocCluster, sep="_")),
      by = "Locus.iteration"
    ) %>%
    mutate(Annotation = str_replace_all(Annotation, "_0$", "_UnannotatedButObserved")) %>%
    mutate(Annotation = str_replace_all(Annotation, "_1$", "_Annotated")) %>%
    mutate(Annotation = str_replace_all(Annotation, "^\\d+?_", "ChromHMM:")) %>%
    filter(!is.na(Annotation)) %>%
    mutate(Annotation = if_else(Annotation==".", "No annotation", Annotation)) %>%
    group_by(Annotation, ClusterGroup, QTL.Set) %>%
    summarise(Sum = sum(PosteriorPr)) %>%
    group_by(ClusterGroup, QTL.Set) %>%
    mutate(FractionAnnotatedPosteriorInGroup = Sum/sum(Sum, na.rm=T)) %>%
    ungroup() %>%
    # select(-ClusterGroup) %>%
    pivot_wider(names_from = "QTL.Set", values_from = c("FractionAnnotatedPosteriorInGroup", "Sum"), id_cols="Annotation") %>%
    mutate(Enrichment = FractionAnnotatedPosteriorInGroup_1/FractionAnnotatedPosteriorInGroup_2) %>%
    mutate(i=i)
  }
  output[["data"]] <- bind_rows(results)
  output[["plot"]] <- output[["data"]] %>%
  group_by(Annotation) %>%
  summarise(mean = exp(mean(log(Enrichment), na.rm=T)),
            Percentile10 = quantile(Enrichment, probs = 0.05, na.rm=T),
            Percentile90 = quantile(Enrichment, probs = 0.95, na.rm=T)) %>%
  ggplot(aes(x=Annotation, y=mean)) +
    # geom_vline(xintercept = 1) +
    geom_col() +
    geom_errorbar(aes(ymin=Percentile10, ymax=Percentile90)) +
    theme_bw() +
    scale_y_continuous(trans='log2') +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
    coord_flip() +
    labs(x = "Annotation", y="Finemapping posterior mass relative fold enrichment\nQTL.Set2 <----> QTL.Set1",
         caption = "Bootsrapped 90% confidence interval")
  return(output)
}

#Test the function
test.results <- PlotDifferentialEnrichmentBetweenQTLSets(dat, H3K27AC_AndEqtlClusters, SplicingClusters)

test.results$data %>% head()

# A tibble: 6 × 7
  Annotation    FractionAnnotat… FractionAnnotat…  Sum_1  Sum_2 Enrichment     i
  <chr>                    <dbl>            <dbl>  <dbl>  <dbl>      <dbl> <int>
1 ChromHMM:Act…          0.0585           0.0438  3.58e2 142.        1.33      1
2 ChromHMM:Het…          0.207            0.144   1.27e3 465.        1.45      1
3 ChromHMM:Ins…          0.00775          0.00438 4.74e1  14.2       1.77      1
4 ChromHMM:Poi…          0.00399          0.00161 2.44e1   5.21      2.48      1
5 ChromHMM:Rep…          0.00124          0.00151 7.58e0   4.88      0.822     1
6 ChromHMM:Rep…          0.0177           0.00780 1.08e2  25.3       2.27      1

test.results$plot +
  labs(y="FoldEnrichment\nSplicingQTLs <----> Chromatin|Expression QTLs")

Another hypothesis I have is that the chromatin-specific sQTLs will more reflect splicing and therefore be more enriched to splice sites, whereas the poly-A specific sQTLs will more reflect isoform stability. Let’s see…

hyprcoloc_results %>%
  separate(phenotype_full, into=c("PC", "P"), sep=";") %>%
  pull(PC) %>% unique()

 [1] "polyA.Splicing"                      "chRNA.IR"                           
 [3] "polyA.IR"                            "MetabolicLabelled.60min.IR"         
 [5] "MetabolicLabelled.30min.IR"          "MetabolicLabelled.30min"            
 [7] "MetabolicLabelled.60min"             "Expression.Splicing"                
 [9] "H3K27AC"                             "H3K4ME3"                            
[11] "H3K4ME1"                             "chRNA.Expression_eRNA"              
[13] "polyA.IR.Subset_YRI"                 "chRNA.IER"                          
[15] "chRNA.IRjunctions"                   "ProCap"                             
[17] "polyA.Splicing.Subset_YRI"           "CTCF"                               
[19] "chRNA.Expression_lncRNA"             "Expression.Splicing.Subset_YRI"     
[21] "chRNA.Splicing"                      "MetabolicLabelled.30min.Splicing"   
[23] "MetabolicLabelled.60min.Splicing"    "chRNA.Expression.Splicing"          
[25] "H3K36ME3"                            "polyA.IER"                          
[27] "polyA.IRjunctions"                   "chRNA.Expression_cheRNA"            
[29] "chRNA.Slopes"                        "MetabolicLabelled.30min.IER"        
[31] "MetabolicLabelled.60min.IER"         "MetabolicLabelled.60min.IRjunctions"
[33] "MetabolicLabelled.30min.IRjunctions" "chRNA.Expression_snoRNA"            

# chRNA.Splicing.PCs <- c("chRNA.Splicing", "chRNA.IER")
# polyA.Splicing.PCs <- c("polyA.Splicing.Subset_YRI", "polyA.IER")

chRNA.Splicing.PCs <- c("chRNA.Splicing")
polyA.Splicing.PCs <- c("polyA.Splicing.Subset_YRI")


chRNA.Specific.SplicingClusters <- hyprcoloc_results %>%
  separate(phenotype_full, into=c("PC", "P"), sep=";") %>%
  group_by(Locus, iteration) %>%
  filter(any(PC %in% chRNA.Splicing.PCs)) %>%
  filter(any(PC %in% polyA.Splicing.PCs)) %>%
  distinct(Locus,iteration) %>%
  unite(Locus.iteration, Locus,iteration) %>%
  mutate(ClusterGroup = "chRNA.Specific.SplicingClusters")

polyA.Specific.SplicingClusters <- hyprcoloc_results %>%
  separate(phenotype_full, into=c("PC", "P"), sep=";") %>%
  group_by(Locus, iteration) %>%
  filter(any(PC %in% polyA.Splicing.PCs)) %>%
  filter(!any(PC %in% chRNA.Splicing.PCs)) %>%
  distinct(Locus,iteration) %>%
  unite(Locus.iteration, Locus,iteration) %>%
  mutate(ClusterGroup = "polyA.Specific.SplicingClusters")

results <- PlotDifferentialEnrichmentBetweenQTLSets(dat, chRNA.Specific.SplicingClusters, polyA.Specific.SplicingClusters, iterations=400)


results$plot +
  labs(y = "FoldEnrichment\npolyA-specific sQTLs <----> chromatin-specific sQTLs")

Now let’s compare eQTLs with colocalizing H3K27AC, versus eQTLs that colocalize with something but not H3K27AC

chRNA.Splicing.PCs <- c("chRNA.Splicing")
polyA.Splicing.PCs <- c("polyA.Splicing.Subset_YRI")


eQTL.Chromatin.Clusters <- hyprcoloc_results %>%
  separate(phenotype_full, into=c("PC", "P"), sep=";") %>%
  group_by(Locus, iteration) %>%
  filter(any(PC == "Expression.Splicing.Subset_YRI")) %>%
  filter(any(PC %in% c("H3K27AC", "H3K4ME3", "H3K4ME1"))) %>%
  distinct(Locus,iteration) %>%
  unite(Locus.iteration, Locus,iteration) %>%
  mutate(ClusterGroup = "eQTL.Chromatin.Clusters")

eQTL.NoChromatin.Clusters <- hyprcoloc_results %>%
  separate(phenotype_full, into=c("PC", "P"), sep=";") %>%
  group_by(Locus, iteration) %>%
  filter(any(PC == "Expression.Splicing.Subset_YRI")) %>%
  filter(!any(PC %in% c("H3K27AC", "H3K4ME3", "H3K4ME1"))) %>%
  distinct(Locus,iteration) %>%
  unite(Locus.iteration, Locus,iteration) %>%
  mutate(ClusterGroup = "eQTL.NoChromatin.Clusters")

intersect(eQTL.Chromatin.Clusters$Locus.iteration, eQTL.NoChromatin.Clusters$Locus.iteration)

character(0)

results <- PlotDifferentialEnrichmentBetweenQTLSets(dat, eQTL.Chromatin.Clusters, eQTL.NoChromatin.Clusters, iterations=400)


results$plot +
  labs(y = "FoldEnrichment\neQTL.NoChromatin.Clusters <----> eQTL.Chromatin.Clusters")

Session information

sessionInfo()

R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so

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       

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

other attached packages:
 [1] boot_1.3-23       data.table_1.14.2 forcats_0.4.0     stringr_1.4.0    
 [5] dplyr_1.0.9       purrr_0.3.4       readr_1.3.1       tidyr_1.2.0      
 [9] tibble_3.1.7      ggplot2_3.3.6     tidyverse_1.3.0  

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.5        lubridate_1.7.4   assertthat_0.2.1  rprojroot_2.0.2  
 [5] digest_0.6.20     utf8_1.1.4        R6_2.4.0          cellranger_1.1.0 
 [9] backports_1.4.1   reprex_0.3.0      evaluate_0.15     highr_0.9        
[13] httr_1.4.1        pillar_1.7.0      rlang_1.0.3       readxl_1.3.1     
[17] rstudioapi_0.10   R.utils_2.9.0     R.oo_1.22.0       rmarkdown_1.13   
[21] labeling_0.3      munsell_0.5.0     broom_1.0.0       compiler_3.6.1   
[25] httpuv_1.5.1      modelr_0.1.8      xfun_0.31         pkgconfig_2.0.2  
[29] htmltools_0.3.6   tidyselect_1.1.2  workflowr_1.6.2   fansi_0.4.0      
[33] crayon_1.3.4      dbplyr_1.4.2      withr_2.4.1       later_0.8.0      
[37] R.methodsS3_1.7.1 grid_3.6.1        jsonlite_1.6      gtable_0.3.0     
[41] lifecycle_1.0.1   DBI_1.1.0         git2r_0.26.1      magrittr_1.5     
[45] scales_1.1.0      cli_3.3.0         stringi_1.4.3     farver_2.1.0     
[49] fs_1.3.1          promises_1.0.1    xml2_1.3.2        ellipsis_0.3.2   
[53] generics_0.1.3    vctrs_0.4.1       tools_3.6.1       glue_1.6.2       
[57] hms_0.5.3         yaml_2.2.0        colorspace_1.4-1  rvest_0.3.5      
[61] knitr_1.39        haven_2.3.1