Tidy LCL data for Jingxin

Identifying cassette exons

identifying GA|GT introns that have a nearby upstream splice acceptor:

MaxCassetteExonLen <- 300

splicing.gagt.donors.granges <- splicing.gagt.df %>%
  dplyr::select(chrom=`#Chrom`, intron.start=start, intron.end=end, strand) %>%
  mutate(start = case_when(
    strand == "+" ~ intron.start,
    strand == "-" ~ intron.end -2
  )) %>%
  mutate(end = start +2) %>%
  dplyr::select(chrom, start, end, strand) %>% 
  distinct() %>%
  add_row(chrom="DUMMY", start=0, end=1, strand="+", .before=1) %>%
  makeGRangesFromDataFrame()

splicing.acceptors.granges <- all.samples.PSI %>%
  dplyr::select(1:6) %>%
  filter(gid %in% splicing.gagt.df$gid) %>%
  filter(!junc %in% splicing.gagt.df$junc) %>%
  dplyr::select(chrom=`#Chrom`, intron.start=start, intron.end=end, strand) %>%
  mutate(start = case_when(
    strand == "+" ~ intron.end -2,
    strand == "-" ~ intron.start -2
  )) %>%
  mutate(end = start +2) %>%
  dplyr::select(chrom, start, end, strand) %>% 
  distinct() %>%
  makeGRangesFromDataFrame()

PrecedingIndexes <- GenomicRanges::follow(splicing.gagt.donors.granges, splicing.acceptors.granges, select="last") %>%
  replace_na(1)

UpstreamAcceptors.df <- cbind(
  splicing.gagt.donors.granges %>% as.data.frame() %>%
    dplyr::rename(chrom=1, start=2, end=3, width=4, strand=5),
  splicing.acceptors.granges[PrecedingIndexes] %>% as.data.frame() %>%
    dplyr::rename(chrom.preceding=1, start.preceding=2, end.preceding=3, width.preceding=4, strand.preceding=5)
) %>%
  filter(as.character(chrom) == as.character(chrom.preceding)) %>%
  filter((strand=="+" & start - end.preceding   < MaxCassetteExonLen) | (strand=="-" & start.preceding - end < MaxCassetteExonLen)) %>%
  mutate(SpliceDonor = case_when(
    strand == "+" ~ paste(chrom, start, strand, sep="."),
    strand == "-" ~ paste(chrom, end, strand, sep=".")
  )) %>%
  mutate(UpstreamSpliceAcceptor = case_when(
    strand == "+" ~ paste(chrom, end.preceding, strand, sep="."),
    strand == "-" ~ paste(chrom, start.preceding, strand, sep=".")
  )) %>%
  dplyr::select(SpliceDonor, UpstreamSpliceAcceptor)

# Add the UpstreamSpliceAcceptor column to the splicing.gagt.df. NA indicates there is no upstream splice acceptor
splicing.gagt.df <- splicing.gagt.df %>%
  mutate(SpliceDonor = case_when(
    strand == "+" ~ paste(`#Chrom`, start, strand, sep="."),
    strand == "-" ~ paste(`#Chrom`, end, strand, sep=".")
    )) %>%
  left_join(UpstreamAcceptors.df, by="SpliceDonor") %>%
  mutate(IntronType = recode(anchor, A="Alt 5'ss", D="Alt 3'ss", DA="Annotated", N="New Intron", NDA="New Splice Site combination"))

So I added a new column called UpstreamSpliceAcceptor that has the coordinates of the nearest upstream 3’ss within 300bp. I have found that this is a good way to filter for the GA|GT introns involved in cassette exon inclusion. If this column has NA, it is unlikely this intron is involved in exon inclusion. Probably, if i wanted to filter for the introns with the highest confidence for being involved in new cassette exons, I would filter for Alt 5’ss with an UpstreamSpliceAcceptor:

splicing.gagt.df %>%
  filter(IntronType == "Alt 5'ss" & !is.na(UpstreamSpliceAcceptor))

# A tibble: 1,033 × 52
   `#Chrom`    start    end junc  gid   strand Branaplam_3160_… Branaplam_1000_…
   <chr>       <dbl>  <dbl> <chr> <chr> <chr>             <dbl>            <dbl>
 1 chr1      1300068 1.30e6 chr1… chr1… -                  0.28             0.53
 2 chr1      1390865 1.39e6 chr1… chr1… -                  7.9              1.2 
 3 chr1      8015692 8.02e6 chr1… chr1… -                  4.4             16   
 4 chr1      8365974 8.40e6 chr1… chr1… -                 49               48   
 5 chr1      9714131 9.72e6 chr1… chr1… +                 33               24   
 6 chr1      9773394 9.81e6 chr1… chr1… -                  1.4              0   
 7 chr1     10178166 1.02e7 chr1… chr1… +                  8.6              0.47
 8 chr1     11243300 1.12e7 chr1… chr1… -                  9.7              6.1 
 9 chr1     11952068 1.20e7 chr1… chr1… +                  8.5              3   
10 chr1     11981853 1.20e7 chr1… chr1… +                 26               21   
# … with 1,023 more rows, and 44 more variables:
#   Branaplam_316_LCL_polyA_1 <dbl>, Branaplam_100_LCL_polyA_1 <dbl>,
#   Branaplam_31.6_LCL_polyA_1 <dbl>, Branaplam_10_LCL_polyA_1 <dbl>,
#   Branaplam_3.16_LCL_polyA_1 <dbl>, Branaplam_1_LCL_polyA_1 <dbl>,
#   C2C5_1000_LCL_polyA_1 <dbl>, C2C5_316_LCL_polyA_1 <dbl>,
#   C2C5_100_LCL_polyA_1 <dbl>, C2C5_31.6_LCL_polyA_1 <dbl>,
#   C2C5_10_LCL_polyA_1 <dbl>, C2C5_3.16_LCL_polyA_1 <dbl>, …

Prepping data for model fitting:

As before, I will include the DMSO data for model-fitting… The drc package is designed to handle data with dose=0, so I will code the DMSO samples as treated with dose=0. There are 3 replicates of DMSO. I will randomly assign 1 DMSO sample for each titration series. I will actually do this on a per-junctino level, so it’s not like branaplam titration series is always paired with the same DMSO replicate… Rather, it will vary randomly by junction.

dmso.data <- splicing.gagt.df %>%
  dplyr::select(junc, contains("LCL")) %>%
  gather("SampleName", "PSI", -junc) %>%
  left_join(sample.list, by="SampleName") %>%
  filter(treatment=="DMSO") %>%
  arrange(junc, rep) %>%
  mutate(dose.nM = 0)
dmso.data$treatment <- replicate(nrow(dmso.data)/3, sample(c("Branaplam", "Risdiplam", "C2C5"), 3)) %>% as.vector()

# Tidy the splicing.gagt.df dataframe
splicing.gagt.df.tidy <- splicing.gagt.df %>%
  dplyr::select(junc, contains("LCL")) %>%
  gather("SampleName", "PSI", contains("LCL")) %>%
  left_join(sample.list, by="SampleName") %>%
  filter(!treatment=="DMSO") %>%
  bind_rows(dmso.data) %>%
  inner_join(
    splicing.gagt.df %>%
        dplyr::select(1:6, seq, Donor.score=score.x, gene_names, gene_ids, SpliceDonor, UpstreamSpliceAcceptor, IntronType),
    by="junc"
  )

Now the data is really tidy for plotting the raw data for a dose-response curve with ggplot2 functions:

splicing.gagt.df.tidy %>%
  filter(junc %in% c("chr22:35267514:35283962:clu_38658_+", "chr17:30183823:30184615:clu_32603_+")) %>%
  ggplot(aes(color=treatment)) +
    geom_line(aes(x=dose.nM, y=PSI)) +
    scale_x_continuous(trans="log1p", limits=c(0, 10000), breaks=c(10000, 3160, 1000, 316, 100, 31.6, 10, 3.16, 1, 0.316, 0)) +
    facet_wrap(~junc, scales = "free_y") +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 45, vjust = 1))  

Note that the x-axis I’m plotting the actual nanomolar dosage, and since the genome-wide effective conentration of these drugs are a bit different, it might be more interpretable to plot them using their within-treatment ranked dosage:

splicing.gagt.df.tidy %>%
  filter(junc %in% c("chr22:35267514:35283962:clu_38658_+", "chr17:30183823:30184615:clu_32603_+")) %>%
  group_by(treatment) %>%
  mutate(DosageRank = dense_rank(dose.nM)) %>%
  ggplot(aes(color=treatment)) +
    geom_line(aes(x=DosageRank, y=PSI)) +
    facet_wrap(~junc, scales = "free_y") +
    theme_bw() +
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1, size=3)) 

In any case, when I do the model fitting, I think it will be more appropriate/useful to use the actual nanomolar dosage as the explanatory variable, so the units of the EC50 are easily interpretable. Keep that in mind when interpreting my EC50 estimate results.

Filtering by spearman correlation coef

As before, I think a reasonable way to pre-filter introns that have a reasonable dose-response curve is to filter for introns that have a spearman correlation coef of at least <0.9 across doses in at least one treatment dose-response series. Let’s still do that filter:

# Add spearman correlation ceof column
splicing.gagt.df.tidy <- splicing.gagt.df.tidy %>%
  nest(-treatment, -junc) %>%
  mutate(cor=map(data,~cor.test(.x$dose.nM, .x$PSI, method = "sp"))) %>%
  mutate(tidied = map(cor, tidy)) %>% 
  unnest(tidied, .drop = T) %>%
  dplyr::select(junc:data, spearman=estimate) %>%
  unnest(data)

# filter the data for model fitting
model.dat.df <- splicing.gagt.df.tidy %>%
  group_by(junc) %>%
  filter(any(spearman > 0.9)) %>%
  ungroup() %>%
  group_by(junc, treatment) %>% 
  nest(-junc)

Model-fitting

this is basically the hackish code i used before. I kept trying to use broom functions to fit things, but I kept running into errors that I don’t feel like figuring out, and the following code seems to get the job done, with a tryCatch statement to print out errors that might occur in the first 100 introns. so you will see, some things error because convergence failed (likely data dosen’t even resemeble a nice dose-response curve), whereas some things error because “0 (non-NA) cases” which is confusing because i previously showed that error pops up even when there are no NAs in my data, but rather when PSI is exactly 0 for a lot of data points… In any case, when those errors pop up, I will just exclude those introns from the modelling results.

Results <- list()
for(i in 1:nrow(model.dat.df)) {
  tryCatch(
    expr = {
      junc <- model.dat.df$junc[i]
      data <- model.dat.df$data[i] %>% as.data.frame()
      fit <- drm(formula = PSI ~ dose.nM,treatment,
                  data = data,
                  # curveid = treatment,
                  fct = LL.4(names=c("Steepness", "LowerLimit", "UpperLimit", "ED50")),
                  pmodels=data.frame(1, 1, 1, treatment)
                  )
      df.out <- fit$coefficients %>%
        as.data.frame() %>%
        rownames_to_column("param") %>%
        rename(`.` = "estimate")
      message("Successfully fitted model.")
      Results[[junc]] <- df.out
      # do stuff with row
    },
    error=function(e){
      if (i < 100){
        cat("ERROR :",conditionMessage(e), junc, "\n")
      }
      })
}

Error in optim(startVec, opfct, hessian = TRUE, method = optMethod, control = list(maxit = maxIt,  : 
  non-finite finite-difference value [5]
ERROR : Convergence failed chr1:12066023:12084464:clu_2210_+ 
ERROR : 0 (non-NA) cases chr1:41028726:41033593:clu_638_- 
Error in optim(startVec, opfct, hessian = TRUE, method = optMethod, control = list(maxit = maxIt,  : 
  non-finite finite-difference value [4]
ERROR : Convergence failed chr1:182874291:182874854:clu_3655_+ 
Error in optim(startVec, opfct, hessian = TRUE, method = optMethod, control = list(maxit = maxIt,  : 
  non-finite finite-difference value [6]
ERROR : Convergence failed chr1:203228562:203273376:clu_1714_- 
ERROR : 0 (non-NA) cases chr10:68446413:68448247:clu_20590_- 
ERROR : 0 (non-NA) cases chr10:124414432:124416736:clu_21009_- 
ERROR : 0 (non-NA) cases chr10:131934767:131934894:clu_20222_+ 
ERROR : 0 (non-NA) cases chr11:33867062:33869346:clu_21326_- 
Error in optim(startVec, opfct, hessian = TRUE, method = optMethod, control = list(maxit = maxIt,  : 
  non-finite finite-difference value [6]
Error in optim(startVec, opfct, hessian = TRUE, method = optMethod, control = list(maxit = maxIt,  : 
  non-finite finite-difference value [6]
Error in optim(startVec, opfct, hessian = TRUE, method = optMethod, control = list(maxit = maxIt,  : 
  non-finite finite-difference value [4]
Error in optim(startVec, opfct, hessian = TRUE, method = optMethod, control = list(maxit = maxIt,  : 
  non-finite finite-difference value [5]
Error in optim(startVec, opfct, hessian = TRUE, method = optMethod, control = list(maxit = maxIt,  : 
  non-finite value supplied by optim
Error in optim(startVec, opfct, hessian = TRUE, method = optMethod, control = list(maxit = maxIt,  : 
  non-finite finite-difference value [4]
Error in optim(startVec, opfct, hessian = TRUE, method = optMethod, control = list(maxit = maxIt,  : 
  non-finite finite-difference value [5]
Error in optim(startVec, opfct, hessian = TRUE, method = optMethod, control = list(maxit = maxIt,  : 
  non-finite finite-difference value [5]
Error in optim(startVec, opfct, hessian = TRUE, method = optMethod, control = list(maxit = maxIt,  : 
  non-finite finite-difference value [6]
Error in optim(startVec, opfct, hessian = TRUE, method = optMethod, control = list(maxit = maxIt,  : 
  non-finite value supplied by optim
Error in optim(startVec, opfct, hessian = TRUE, method = optMethod, control = list(maxit = maxIt,  : 
  non-finite finite-difference value [5]
Error in optim(startVec, opfct, hessian = TRUE, method = optMethod, control = list(maxit = maxIt,  : 
  non-finite finite-difference value [1]
Error in optim(startVec, opfct, hessian = TRUE, method = optMethod, control = list(maxit = maxIt,  : 
  non-finite finite-difference value [6]

ModelFits.Coefficients <- bind_rows(Results, .id="junc") %>%
  separate(param, into=c("param", "treatment"), sep=":") %>%
  pivot_wider(names_from=c("param", "treatment"), values_from="estimate") %>%
  dplyr::rename("LowerLimit" =`LowerLimit_(Intercept)`, "UpperLimit" =`UpperLimit_(Intercept)`, "Steepness" =`Steepness_(Intercept)`)

Now let’s do some merging of the data into a useful format for sharing with Jingxin, some quick exploration, and write out to a file for sharing:

EC50 estimate data exploration

How many errors were thrown, compared to how many introns I attempted to model:

#Number introns with results; no error
nrow(ModelFits.Coefficients)

[1] 598

#Number introns attempted to model
nrow(model.dat.df)

[1] 653

Ok, so most things fit succesfully.

First let’s look at the distribution of the steepness parameter estimates. The way the model parameterized is such that a negative steepness parameter corresponds to a positive dose-response relationship (increase doses means increasing the expected response). So I expect these to be largely negative, since I only modelled GA|GT introns which I expect to have positive slope dose-response relationships. Keep in mind I kept this parameter fixed across treatments for simpler models with less degrees of freedom.

ModelFits.Coefficients %>%
  ggplot(aes(x=Steepness)) +
  stat_ecdf()

Now let’s plot a random sample of dose-response curves, with the EC50 parameter estimates as vertical lines:

sample_n_of <- function(data, size, ...) {
  dots <- quos(...)
  
  group_ids <- data %>% 
    group_by(!!! dots) %>% 
    group_indices()
  
  sampled_groups <- sample(unique(group_ids), size)
  
  data %>% 
    filter(group_ids %in% sampled_groups)
}

set.seed(0)
model.dat.df %>%
  unnest(data) %>%
  ungroup() %>%
  sample_n_of(10, junc) %>%
  inner_join(
    bind_rows(Results, .id="junc") %>%
      separate(param, into=c("param", "treatment"), sep=":") %>%
      filter(param=="ED50"),
    by=c("junc", "treatment")) %>%
ggplot(aes(color=treatment)) +
  geom_vline(
    data = (. %>%
      distinct(junc, treatment, estimate)),
    aes(xintercept=estimate, color=treatment),
    linetype='dashed') +
  geom_label(
    data = (. %>%
      distinct(junc, treatment, estimate) %>%
      group_by(junc) %>%
      mutate(vjust = row_number()) %>%
      ungroup()),
    aes(vjust = vjust, label=signif(estimate,3)),
    y=Inf, x=-Inf, size=2.5, hjust=-0.1, label.padding=unit(0.1, "lines")
  ) +
  geom_line(aes(x=dose.nM, y=PSI)) +
  scale_x_continuous(trans="log1p", limits=c(0, 10000), breaks=c(10000, 3160, 1000, 316, 100, 31.6, 10, 3.16, 1, 0.316, 0)) +
  facet_wrap(~junc, scales = "free_y") +
  labs(title="All modeled introns, with EC50 estimate") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1, size=3))

I think the EC50 estimates are actually pretty reasonable, given the limitations of the noisy data.

Now let’s plot scatter plots of the EC50 estimates across conditions

diag_limitrange <- function(data, mapping, ...) { 
  ggplot(data = data, mapping = mapping, ...) + 
    geom_density(...) + 
    scale_x_continuous(trans="log10") +
    theme_bw()
}

upper_point <- function(data, mapping, ...) { 
  ggplot(data = data, mapping = mapping, ...) + 
    geom_point(..., alpha=0.05) + 
    scale_y_continuous(trans="log10") +
    scale_x_continuous(trans="log10") +
    geom_abline() +
    theme_bw()
}

ModelFits.Coefficients %>%
  dplyr::select(junc, contains("ED")) %>%
    ggpairs(
    title="ED50 across treatments",
    columns = 2:4,
    upper=list(continuous = wrap("cor", method = "spearman", hjust=0.7)),
    lower=list(continuous = upper_point),
    diag=list(continuous = diag_limitrange))

Ok, so keep in mind the offset from the diagonal is somewhat expected, as the ED50 units are in nanomolar and these drugs have different genome-wide effective concentrations. And as expected, the C2C5 looks a lot like risdiplam.

Lastly, since I previously used PC2 loading as a proxy for molecule specific effects, let’s look at how that correlates with the difference in EC50 between risdiplam and branaplam.

ModelFits.Coefficients %>%
  mutate(ED50.Ratio = ED50_Branaplam/ED50_Risdiplam) %>%
  inner_join(
    pca.results.splicing$rotation %>%
    as.data.frame() %>%
    dplyr::select(PC2) %>%
    rownames_to_column("junc"),
    by="junc"
  ) %>%
  ggplot(aes(x=ED50.Ratio, y=PC2)) +
  geom_point() +
  scale_x_continuous(trans="log10", breaks=10**c(-3:3), limits=c(10**-3, 10**3)) +
  theme_bw() +
  labs(x="log2(EC50.Branaplam/EC50.Risdiplam", y="PC2 loading", title="PC2 as a proxy for drug-specific response")

Next let’s plot dose response curves for the top20 juncs by difference in EC50. Before seeing these plots, I’m not sure whether I should expect to be outliers with really questionable dose response curves (keep in mind, the standard errors for these parameter estimates I expect to be biased and I didn’t even report them), or whether they are really believable results

ModelFits.Coefficients %>%
  mutate(ED50.Ratio = ED50_Branaplam/ED50_Risdiplam) %>%
  arrange(ED50.Ratio) %>%
  head(20) %>%
  inner_join(splicing.gagt.df.tidy, by="junc") %>%
ggplot(aes(color=treatment)) +
  geom_line(aes(x=dose.nM, y=PSI)) +
  scale_x_continuous(trans="log1p", limits=c(0, 10000), breaks=c(10000, 3160, 1000, 316, 100, 31.6, 10, 3.16, 1, 0.316, 0)) +
  facet_wrap(~junc, scales = "free_y") +
  labs(title="Top20 GA|GT introns by ED50 ratio") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1, size=3))

ModelFits.Coefficients %>%
  mutate(ED50.Ratio = ED50_Branaplam/ED50_Risdiplam) %>%
  arrange(ED50.Ratio) %>%
  tail(20) %>%
  inner_join(splicing.gagt.df.tidy, by="junc") %>%
ggplot(aes(color=treatment)) +
  geom_line(aes(x=dose.nM, y=PSI)) +
  scale_x_continuous(trans="log1p", limits=c(0, 10000), breaks=c(10000, 3160, 1000, 316, 100, 31.6, 10, 3.16, 1, 0.316, 0)) +
  facet_wrap(~junc, scales = "free_y") +
  labs(title="Bottom20 GA|GT introns by ED50 ratio") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1, size=3))

Hmm, by in large the results are as expected, but you should certainly look at the dose response curve yourself before really believing any particular junction. In the future I will fix this with more careful modelling, and filtering, and estimates of standard error. But for now this I think this will do.

Now I will write out the results in a useful way for sharing with Jingxin.

# Data to recreate dose-response data plots
splicing.gagt.df.tidy %>%
  dplyr::select(junc, SampleName, old.sample.name, treatment, dose.nM, cell.type, PSI) %>%
  write_tsv("../output/PSI.ForGAGTIntrons.tidy.txt.gz")



# Data with intron-level metadata and EC50 estimates
splicing.gagt.df.tidy %>% distinct(junc, .keep_all=T) %>%
  dplyr::select(junc, `#Chrom`:IntronType) %>%
  left_join(ModelFits.Coefficients, by="junc") %>%
  left_join(
    splicing.gagt.df.tidy %>%
      dplyr::select(junc, treatment,spearman) %>%
      distinct() %>%
      pivot_wider(names_from="treatment", values_from="spearman", names_prefix="spearman.coef."),
    by="junc"
  ) %>%
  write_tsv("../output/EC50Estimtes.FromPSI.txt.gz")

• So I put the output files here in output/EC50Estimtes.FromPSI.txt.gz and output/PSI.ForGAGTIntrons.tidy.txt.gz. They should also be tracked on github here

• Email me if it is anything is unclear =)

Session information

sessionInfo()

R version 3.6.1 (2019-07-05)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

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         LC_TIME=C           
 [4] LC_COLLATE=C         LC_MONETARY=C        LC_MESSAGES=C       
 [7] LC_PAPER=C           LC_NAME=C            LC_ADDRESS=C        
[10] LC_TELEPHONE=C       LC_MEASUREMENT=C     LC_IDENTIFICATION=C 

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

other attached packages:
 [1] GGally_1.4.0         broom_1.0.0          drc_3.0-1           
 [4] MASS_7.3-51.4        GenomicRanges_1.36.1 GenomeInfoDb_1.20.0 
 [7] IRanges_2.18.1       S4Vectors_0.22.1     BiocGenerics_0.30.0 
[10] forcats_0.4.0        stringr_1.4.0        dplyr_1.0.9         
[13] purrr_0.3.4          readr_1.3.1          tidyr_1.2.0         
[16] tibble_3.1.7         ggplot2_3.3.6        tidyverse_1.3.0     

loaded via a namespace (and not attached):
 [1] bitops_1.0-6           fs_1.5.2               lubridate_1.7.4       
 [4] RColorBrewer_1.1-2     httr_1.4.4             rprojroot_2.0.2       
 [7] tools_3.6.1            backports_1.4.1        utf8_1.1.4            
[10] R6_2.4.0               DBI_1.1.0              colorspace_1.4-1      
[13] withr_2.5.0            tidyselect_1.1.2       curl_3.3              
[16] compiler_3.6.1         git2r_0.26.1           cli_3.3.0             
[19] rvest_0.3.5            xml2_1.3.2             sandwich_3.0-2        
[22] labeling_0.3           scales_1.1.0           mvtnorm_1.1-3         
[25] digest_0.6.20          foreign_0.8-71         rmarkdown_1.13        
[28] rio_0.5.16             XVector_0.24.0         pkgconfig_2.0.2       
[31] htmltools_0.5.3        plotrix_3.8-2          dbplyr_1.4.2          
[34] fastmap_1.1.0          highr_0.9              rlang_1.0.5           
[37] readxl_1.3.1           rstudioapi_0.14        farver_2.1.0          
[40] generics_0.1.3         zoo_1.8-10             jsonlite_1.6          
[43] gtools_3.9.2.2         zip_2.2.1              car_3.0-5             
[46] RCurl_1.98-1.1         magrittr_1.5           GenomeInfoDbData_1.2.1
[49] Matrix_1.2-18          Rcpp_1.0.5             munsell_0.5.0         
[52] fansi_0.4.0            abind_1.4-5            lifecycle_1.0.1       
[55] multcomp_1.4-19        stringi_1.4.3          yaml_2.2.0            
[58] carData_3.0-5          zlibbioc_1.30.0        plyr_1.8.4            
[61] grid_3.6.1             promises_1.0.1         crayon_1.3.4          
[64] lattice_0.20-38        splines_3.6.1          haven_2.3.1           
[67] hms_0.5.3              knitr_1.39             pillar_1.7.0          
[70] codetools_0.2-18       reprex_0.3.0           glue_1.6.2            
[73] evaluate_0.15          data.table_1.14.2      modelr_0.1.8          
[76] vctrs_0.4.1            httpuv_1.5.1           cellranger_1.1.0      
[79] gtable_0.3.0           reshape_0.8.8          assertthat_0.2.1      
[82] xfun_0.31              openxlsx_4.1.0.1       later_0.8.0           
[85] survival_2.44-1.1      workflowr_1.6.2        TH.data_1.1-1         
[88] ellipsis_0.3.2