Human Origins Array Global Results

Greedy

Here is the snakemake rule I used for running flashier:

    run:
        R("""
          # read the genotype matrix
          Y = t(lfa:::read.bed('{params.bed}'))
          # number of individuals
          n = nrow(Y)
          
          # run greedy flash
          flash_fit = flashier::flashier(Y,
                                         greedy.Kmax=20,
                                         prior.type=c('nonnegative', 'point.normal'),
                                         var.type=2,
                                         fix.dim=list(1),
                                         fix.idx=list(1:n),
                                         fix.vals=list(rep(1, n)))
          # save the rds
          saveRDS(flash_fit, '{output.rds}')
          """)

Lets first read the greedy flashier fit

flash_fit = readRDS("../output/flash_greedy/hoa_global_ld/HumanOriginsPublic2068_maf_geno_auto_ldprune.rds")
K = ncol(flash_fit$loadings$normalized.loadings[[1]]) 
n = nrow(flash_fit$loadings$normalized.loadings[[1]])
p = nrow(flash_fit$loadings$normalized.loadings[[2]])
print(K)

[1] 21

print(n)

[1] 2068

print(p)

[1] 165468

Lets now plot the distribution of factors for each drift event

# read factors
delta_df = as.data.frame(flash_fit$loadings$normalized.loadings[[2]])
colnames(delta_df)[1:K] = 1:K 

# add snp meta data
bim_df = read.table("../data/raw/NearEastPublic/HumanOriginsPublic2068_maf_geno_auto_ldprune.bim", header=F)
colnames(bim_df) = c("chrom", "rsid", "cm", "pos", "a1", "a2")
delta_df$chrom = bim_df$chrom
delta_df$pos = bim_df$pos
delta_df$rsid = bim_df$rsid

# gather the data.frame for plotting
delta_gath_df = delta_df %>% 
                gather(K, value, -chrom, -pos, -rsid) %>%
                filter(K!=1)

# plot the factors
K_ = K
p_fct = ggplot(delta_gath_df, aes(x=value)) + 
        geom_histogram() + 
        facet_wrap(~factor(K, levels=2:K_), scales = "free") + 
        labs(fill="K") + 
        scale_x_continuous(breaks = scales::pretty_breaks(n = 3)) +
        scale_y_continuous(breaks = scales::pretty_breaks(n = 3)) + 
        theme_bw()
p_fct

Past versions of flash-greedy-ld-viz-factors-1.png

We can see the later factors tend to get sparser! Here is a plot of the “proportion of variance explained” of each factor:

qplot(2:K, flash_fit$pve[2:K]) + ylab("Proportion of Varaince Explained") + xlab("K") + theme_bw()

Past versions of flash-greedy-ld-viz-pve-1.png

print(flash_fit$pve)

 [1] 0.4763895014 0.0180272272 0.0221803113 0.0102965471 0.0025143344
 [6] 0.0045277929 0.0015659276 0.0014218575 0.0011760893 0.0008608887
[11] 0.0007633150 0.0006690835 0.0005002149 0.0002133497 0.0001920124
[16] 0.0002281114 0.0001847405 0.0002903768 0.0002700529 0.0001573717
[21] 0.0001176524

It looks like the PVE drops off at around 14 or so? I setup the flashier run so it estimates a SNP specific precision term. Here is a histogram of fitted variances:

qplot(1/flash_fit$fit$est.tau) + xlab("Estimated Variance") + ylab("Count") + theme_bw()

Past versions of flash-greedy-ld-viz-var-1.png

Lets now look the the fitted means:

mu = sqrt(flash_fit$loadings$scale.constant[1]) * delta_df$`1`
qplot(mu) + xlab("Estimated Mean") + ylab("Count") + theme_bw()

Past versions of flash-greedy-ld-viz-mean-1.png

These plots looks about reasonable as each of the SNP variances should roughly be interpreted as average heterozygosity \(\approx 2p(1-p)\)? The mean term should roughly be interpreted as the mean minor allele frequency at the SNP and thus we should see a quadratic relationship with the estimated variance:

d1 = flash_fit$loadings$scale.constant[1]
mv_df = data.frame(var=1/flash_fit$fit$est.tau, mu=mu, chrom=bim_df$chrom)
p_mv = ggplot(mv_df, aes(x=mu, y=var)) + 
       geom_point() + 
       xlab("Rescaled Estimated Mean") + ylab("Estimated Variance") + 
       scale_alpha(guide = "none") + 
       stat_function(fun = function(x){return(2*x*(1-x))}, color="red") + 
       xlim(0, .5) + 
       theme_bw()
p_mv

Past versions of flash-greedy-ld-viz-mv-1.png

Most of the SNPs have a mean-variance relationship expected under a simple Binomial model for the genotypes i.e. \(y_{ij} \sim Binomial(2, p_{ij})\). I wonder if there is anything “special” going on with the high variance SNP (will explore this later). Lets now take a look at the loadings. First we setup a data.frame that we can work with:

# read the meta data
meta_df = read.table("../data/meta/HumanOriginsPublic2068.meta", sep="\t", header=T)

# setup loadings data.frame
l_df = as.data.frame(flash_fit$loadings$normalized.loadings[[1]])
K = ncol(l_df)
l_df = cbind(l_df, meta_df)
pops = unique(l_df$Simple.Population.ID) # all unique pop labels

# join with the meta data
l_df = l_df %>% arrange(Region, Simple.Population.ID) # sort by region then by population
l_df$ID = factor(l_df$ID, levels = l_df$ID) # make sure the ids are sorted
colnames(l_df)[1:K] = 1:K

head(l_df)

           1            2           3          4            5            6
1 0.02198997 3.983040e-06 0.009742430 0.02223268 2.480331e-05 1.933025e-05
2 0.02198997 3.907523e-06 0.009734339 0.02268403 2.577005e-05 1.894343e-05
3 0.02198997 3.849485e-06 0.006489710 0.02495816 2.398744e-05 2.021861e-05
4 0.02198997 3.869767e-06 0.009428810 0.02313380 2.421315e-05 1.971725e-05
5 0.02198997 3.889654e-06 0.009863181 0.02289416 2.525070e-05 2.054919e-05
6 0.02198997 3.977201e-06 0.009502596 0.02229676 2.420936e-05 2.011595e-05
            7            8            9         10           11
1 0.011770312 3.967352e-05 2.832321e-05 0.05538038 3.960260e-05
2 0.010829446 4.011844e-05 2.957381e-05 0.04942202 4.264933e-05
3 0.010805872 3.752266e-05 2.607744e-05 0.05509969 4.206464e-05
4 0.010939192 3.790906e-05 2.729509e-05 0.05984811 3.902960e-05
5 0.008271378 3.933508e-05 2.670218e-05 0.05788346 3.936151e-05
6 0.011583697 3.722886e-05 2.807376e-05 0.05629556 3.791417e-05
           12           13           14         15           16
1 0.014925313 6.531321e-05 7.771835e-05 0.06157679 7.687870e-05
2 0.009683119 8.181238e-05 7.817497e-05 0.04760745 7.947480e-05
3 0.011611347 7.275697e-05 7.761959e-05 0.06522022 7.812401e-05
4 0.016421998 6.915482e-05 7.413616e-05 0.07187187 8.095212e-05
5 0.015993345 6.391435e-05 7.548730e-05 0.07753460 8.029139e-05
6 0.015481074 6.697399e-05 7.687836e-05 0.06055225 7.798898e-05
            17           18           19           20           21
1 0.0001113319 0.0004044934 8.802684e-05 8.731966e-05 8.756886e-05
2 0.0001220043 0.0002632577 3.821495e-04 1.139101e-04 8.716143e-05
3 0.0022902140 0.0002782860 9.505465e-05 2.280344e-04 8.567540e-05
4 0.0001352564 0.0001225972 8.289963e-05 7.181312e-04 8.324215e-05
5 0.0001010870 0.0001548211 1.957986e-04 7.474813e-03 8.673768e-05
6 0.0001327867 0.0005840805 7.520526e-05 9.241005e-05 8.510847e-05
             ID Simple.Population.ID Verbose.Population.ID Region Country
1 Algerian43A22             Algerian              Algerian Africa Algeria
2 Algerian43A21             Algerian              Algerian Africa Algeria
3 Algerian43A34             Algerian              Algerian Africa Algeria
4 Algerian43A13             Algerian              Algerian Africa Algeria
5 Algerian43A24             Algerian              Algerian Africa Algeria
6 Algerian43A32             Algerian              Algerian Africa Algeria
  Latitude Longitude Samples Passed.QC Contributor
1     36.8         3       7         7 David Comas
2     36.8         3       7         7 David Comas
3     36.8         3       7         7 David Comas
4     36.8         3       7         7 David Comas
5     36.8         3       7         7 David Comas
6     36.8         3       7         7 David Comas

Its hard to find a color scale that can sufficiently visualize all of the loadings in a single plot. Instead I just split the loadings up into two plots (K=2,…,11) and (K=12,…,21). Lets first visualize loadings 2 through 12:

# gather the data.frame for plotting
l_gath_df = l_df %>% 
            gather(K, value, -ID, -Verbose.Population.ID, -Simple.Population.ID, 
                   -Region, -Country, -Latitude, -Longitude, -Samples,
                   -Passed.QC,  -Contributor) %>% 
            filter(K %in% paste0(2:11))

# Africa
africa_pops = get_pops(meta_df, "Africa")
p_africa = positive_structure_plot(l_gath_df %>% filter(Region == "Africa"), 
                                   africa_pops, colset="Set3", label_size=5) +
           ggtitle("Africa") + theme(plot.title = element_text(size=6))

# America
america_pops = get_pops(meta_df, "America")
p_america = positive_structure_plot(l_gath_df %>% filter(Region == "America"), 
                                    america_pops, colset="Set3", label_size=5) + 
            ggtitle("America") + theme(plot.title = element_text(size=6))

# Central Asia Siberia
central_asia_siberia_pops = get_pops(meta_df, "CentralAsiaSiberia")
p_central_asia_siberia = positive_structure_plot(l_gath_df %>% filter(Region == "CentralAsiaSiberia"), 
                                                 central_asia_siberia_pops, colset="Set3", label_size=5) + 
                         ggtitle("CentralAsiaSiberia") + theme(plot.title = element_text(size=6))

# East Asia
east_asia_pops = get_pops(meta_df, "EastAsia")
p_east_asia = positive_structure_plot(l_gath_df %>% filter(Region == "EastAsia"), 
                                      east_asia_pops, colset="Set3", label_size=5) + 
              ggtitle("EastAsia") + theme(plot.title = element_text(size=6))

# South Asia
south_asia_pops = get_pops(meta_df, "SouthAsia")
p_south_asia= positive_structure_plot(l_gath_df %>% filter(Region == "SouthAsia"),
                                      south_asia_pops, colset="Set3", label_size=5) + 
              ggtitle("SouthAsia") + theme(plot.title = element_text(size=6))

# West Eurasia
west_eurasia_pops = get_pops(meta_df, "WestEurasia")
p_west_eurasia = positive_structure_plot(l_gath_df %>% filter(Region == "WestEurasia"), 
                                         west_eurasia_pops, colset="Set3", label_size=5) + 
                 ggtitle("WestEurasia") + theme(plot.title = element_text(size=6))

# Oceania
oceania_pops = get_pops(meta_df, "Oceania")
p_oceania = positive_structure_plot(l_gath_df %>% filter(Region == "Oceania"), 
                                    oceania_pops, colset="Set3", label_size=5) + 
            ggtitle("Oceania") + theme(plot.title = element_text(size=6))

# Global
p = cowplot::plot_grid(p_africa, p_west_eurasia, p_central_asia_siberia, p_america, p_east_asia, p_south_asia, p_oceania, 
                       rel_heights = c(1.2, 1.3, 1, 1, 1, 1, 1.1),
                       nrow = 7, align = "v") 
p

Past versions of flash-greedy-ld-viz-loadings-2-11-1.png

Lets now visualize loadings 12 to 21 (be careful: there is no connection to the colors in the last plot):

# gather the data.frame for plotting
l_gath_df = l_df %>% 
            gather(K, value, -ID, -Verbose.Population.ID, -Simple.Population.ID, 
                   -Region, -Country, -Latitude, -Longitude, -Samples,
                   -Passed.QC,  -Contributor) %>% 
            filter(K %in% paste0(12:21))

# Africa
africa_pops = get_pops(meta_df, "Africa")
p_africa = positive_structure_plot(l_gath_df %>% filter(Region == "Africa"), 
                                   africa_pops, colset="Set3", label_size=5) +
           ggtitle("Africa") + theme(plot.title = element_text(size=6))

# America
america_pops = get_pops(meta_df, "America")
p_america = positive_structure_plot(l_gath_df %>% filter(Region == "America"), 
                                    america_pops, colset="Set3", label_size=5) + 
            ggtitle("America") + theme(plot.title = element_text(size=6))

# Central Asia Siberia
central_asia_siberia_pops = get_pops(meta_df, "CentralAsiaSiberia")
p_central_asia_siberia = positive_structure_plot(l_gath_df %>% filter(Region == "CentralAsiaSiberia"), 
                                                 central_asia_siberia_pops, colset="Set3", label_size=5) + 
                         ggtitle("CentralAsiaSiberia") + theme(plot.title = element_text(size=6))

# East Asia
east_asia_pops = get_pops(meta_df, "EastAsia")
p_east_asia = positive_structure_plot(l_gath_df %>% filter(Region == "EastAsia"), 
                                      east_asia_pops, colset="Set3", label_size=5) + 
              ggtitle("EastAsia") + theme(plot.title = element_text(size=6))

# South Asia
south_asia_pops = get_pops(meta_df, "SouthAsia")
p_south_asia= positive_structure_plot(l_gath_df %>% filter(Region == "SouthAsia"),
                                      south_asia_pops, colset="Set3", label_size=5) + 
              ggtitle("SouthAsia") + theme(plot.title = element_text(size=6))

# West Eurasia
west_eurasia_pops = get_pops(meta_df, "WestEurasia")
p_west_eurasia = positive_structure_plot(l_gath_df %>% filter(Region == "WestEurasia"), 
                                         west_eurasia_pops, colset="Set3", label_size=5) + 
                 ggtitle("WestEurasia") + theme(plot.title = element_text(size=6))

# Oceania
oceania_pops = get_pops(meta_df, "Oceania")
p_oceania = positive_structure_plot(l_gath_df %>% filter(Region == "Oceania"), 
                                    oceania_pops, colset="Set3", label_size=5) + 
            ggtitle("Oceania") + theme(plot.title = element_text(size=6))

# Global
p = cowplot::plot_grid(p_africa, p_west_eurasia, p_central_asia_siberia, p_america, p_east_asia, p_south_asia, p_oceania, 
                       rel_heights = c(1.2, 1.3, 1, 1, 1, 1, 1.1),
                       nrow = 7, align = "v") 
p

Past versions of flash-greedy-ld-viz-loadings-12-21-1.png

Its kinda interesting to see that some populations have zero loading on later factors. Its also interesting to see a lot of population specific factors arising. This would be difficult to visualize see if using a single plot for all the factors.

ADMIXTURE

Lets visualize ADMIXTURE with 9 factors which should roughly align to the first plot i.e. FLASH with 2,…,11 (be careful: there is no connection to the colors in the last plot):

l_df = read.table("../output/admixture/hoa_global_ld/HumanOriginsPublic2068_maf_geno_auto_ldprune.K9r1.Q", sep=" ", header=F)
K = ncol(l_df)
l_df = cbind(l_df, meta_df)
pops = unique(l_df$Simple.Population.ID) # all unique pop labels
l_df = l_df %>% arrange(Region, Simple.Population.ID) # sort by region then by population
l_df$ID = factor(l_df$ID, levels = l_df$ID) # make sure the ids are sorted
colnames(l_df)[1:K] = 1:K

# gather the data.frame for plotting
l_gath_df = l_df %>% 
            gather(K, value, -ID, -Verbose.Population.ID, -Simple.Population.ID, 
                   -Region, -Country, -Latitude, -Longitude, -Samples,
                   -Passed.QC,  -Contributor) 

# Africa
africa_pops = get_pops(meta_df, "Africa")
p_africa = positive_structure_plot(l_gath_df %>% filter(Region == "Africa"), 
                                   africa_pops, colset="Set3", label_size=5) +
           ggtitle("Africa") + theme(plot.title = element_text(size=6))

# America
america_pops = get_pops(meta_df, "America")
p_america = positive_structure_plot(l_gath_df %>% filter(Region == "America"), 
                                    america_pops, colset="Set3", label_size=5) + 
            ggtitle("America") + theme(plot.title = element_text(size=6))

# Central Asia Siberia
central_asia_siberia_pops = get_pops(meta_df, "CentralAsiaSiberia")
p_central_asia_siberia = positive_structure_plot(l_gath_df %>% filter(Region == "CentralAsiaSiberia"), 
                                                 central_asia_siberia_pops, colset="Set3", label_size=5) + 
                         ggtitle("CentralAsiaSiberia") + theme(plot.title = element_text(size=6))

# East Asia
east_asia_pops = get_pops(meta_df, "EastAsia")
p_east_asia = positive_structure_plot(l_gath_df %>% filter(Region == "EastAsia"), 
                                      east_asia_pops, colset="Set3", label_size=5) + 
              ggtitle("EastAsia") + theme(plot.title = element_text(size=6))

# South Asia
south_asia_pops = get_pops(meta_df, "SouthAsia")
p_south_asia= positive_structure_plot(l_gath_df %>% filter(Region == "SouthAsia"),
                                      south_asia_pops, colset="Set3", label_size=5) + 
              ggtitle("SouthAsia") + theme(plot.title = element_text(size=6))

# West Eurasia
west_eurasia_pops = get_pops(meta_df, "WestEurasia")
p_west_eurasia = positive_structure_plot(l_gath_df %>% filter(Region == "WestEurasia"), 
                                         west_eurasia_pops, colset="Set3", label_size=5) + 
                 ggtitle("WestEurasia") + theme(plot.title = element_text(size=6))

# Oceania
oceania_pops = get_pops(meta_df, "Oceania")
p_oceania = positive_structure_plot(l_gath_df %>% filter(Region == "Oceania"), 
                                    oceania_pops, colset="Set3", label_size=5) + 
            ggtitle("Oceania") + theme(plot.title = element_text(size=6))
# Global
p = cowplot::plot_grid(p_africa, p_west_eurasia, p_central_asia_siberia, p_america, p_east_asia, p_south_asia, p_oceania, 
                       rel_heights = c(1.2, 1.3, 1, 1, 1, 1, 1.1),
                       nrow = 7, align = "v") 
p

There is a lot that one can compare between the ADMIXTURE and FLASH results. A high level observation seems that the ADMIXTURE results look a bit more clustered i.e. the Americas and East Asia look like they are explained mostly by 1 or 2 factors whereas FLASH uses 3-4. Its hard to tell be it seems that this is true in many of the super regions … ADMIXTURE tends use fewer factors to explain population structure in each region, leading to a more clustered result?

SNP Weights

Lets take a closer look at the drift factors to see if they are clustering in particular regions of the genome. As a first pass I take the top 5% of SNPs weighted on each drift event (to be clear I ignore the sign of each SNP). I would like to use the lfsr here but I will return to later. I then made a Manhatten plot for each chromosome and factor:

delta_thresh_df = delta_df %>%
                  gather(K, value, -chrom, -pos, -rsid) %>%
                  filter(K %in% paste0(2:11)) %>% 
                  group_by(K) %>%
                  summarise(t_95 = quantile(abs(value), probs = .95)[1],
                            t_99 = quantile(abs(value), probs = .99)[1])

# gather the data.frame for plotting
delta_gath_df = delta_df %>% 
                gather(K, value, -chrom, -pos, -rsid) %>%
                filter(K %in% paste0(2:11)) %>% 
                inner_join(delta_thresh_df, on="K") %>%
                filter(abs(value) > t_95) %>%
                filter(chrom %in% 1:22)

Joining, by = "K"

p = ggplot(delta_gath_df, aes(x=pos, y=abs(value), color=factor(K))) + geom_point(size=.5, alpha=.7) + 
    scale_color_brewer(palette = "Set3") + 
    facet_grid(factor(K, levels=2:11)~factor(chrom, levels = 1:22), scales="free") +
    geom_hline(data=delta_gath_df %>% distinct(chrom, K, t_99), aes(yintercept = t_99), linetype="dashed") +
    theme(axis.title.x=element_blank(), axis.text.x=element_blank(), axis.ticks.x=element_blank(),
          panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
    guides(color=FALSE) +
    xlab("Position") + 
    ylab("|Factor|")
p

We can see there are some regions on the chromosomes that are peaky as well as some regions that have no “outliers” at all. I then took the top 5 outliers in each factor and annotated them with some functional information:

grch37.snp = useMart(biomart="ENSEMBL_MART_SNP", host="grch37.ensembl.org", path="/biomart/martservice",dataset="hsapiens_snp")
grch37 = useMart(biomart="ENSEMBL_MART_ENSEMBL", host="grch37.ensembl.org", path="/biomart/martservice", dataset="hsapiens_gene_ensembl")

delta_tophit_df = delta_df %>%
                  gather(K, value, -chrom, -pos, -rsid) %>%
                  filter(K %in% paste0(2:21)) %>% 
                  group_by(K) %>%
                  top_n(5, abs(value))

table1 <- getBM(attributes=c('refsnp_id', 'chrom_start', 'minor_allele_freq', 'ensembl_gene_stable_id',
                             'consequence_type_tv', "associated_gene"), 
                filters = "snp_filter", 
                values = delta_tophit_df$rsid, 
                mart = grch37.snp)

table1$ensembl_gene_id = table1$ensembl_gene_stable_id
table2 <- getBM(attributes = c("ensembl_gene_id", "external_gene_name", "description"),
                filters = "ensembl_gene_id", 
                values =  table1$ensembl_gene_stable_id, 
                mart = grch37)

anno_df = table1 %>% left_join(table2, on="ensembl_gene_id") %>% mutate(rsid=refsnp_id) %>% 
          inner_join(delta_tophit_df, on="rsid") 

Joining, by = "ensembl_gene_id"

Joining, by = "rsid"

Warning: Column `rsid` joining character vector and factor, coercing into
character vector

print(unique(anno_df$external_gene_name))

 [1] "CCNL2"          "CYP4B1"         NA               "RP11-280O1.2"  
 [5] "KIAA1614"       "C1orf132"       "CD34"           "FRMD4A"        
 [9] "SGPL1"          "CCSER2"         "HPSE2"          "NAV2"          
[13] "SLC22A10"       "STT3A-AS1"      "STT3A"          "HOXC13-AS"     
[17] "CNOT2"          "RP11-114H23.1"  "APPL2"          "LINC00641"     
[21] "HERC2"          "RP11-109D20.1"  "SORD"           "PRKCB"         
[25] "PKD1L2"         "RP11-1407O15.2" "CACNB1"         "RP11-515E23.1" 
[29] "RP11-1055B8.4"  "MFSD12"         "IL1R1"          "TMEM87B"       
[33] "RAB3GAP1"       "ZRANB3"         "ITGB6"          "TLK1"          
[37] "PARD3B"         "SF3A1"          "DRG1"           "TEF"           
[41] "KAT2B"          "SEMA3F"         "FRMD4B"         "MYH15"         
[45] "RP11-103J17.1"  "RP11-696N14.1"  "SLC45A2"        "NDUFAF2"       
[49] "DPYSL3"         "UNC5A"          "AKAP12"         "TIAM2"         
[53] "TMEM248"        "SMARCA2"        "GLIS3"          "RCL1"          
[57] "PTPN3"          "TNFSF8"        

kable(anno_df %>% dplyr::select(rsid, chrom, pos, external_gene_name, K, value) %>% arrange(desc(K, value)))

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: macOS  10.14.2

Matrix products: default
BLAS/LAPACK: /Users/jhmarcus/miniconda3/lib/R/lib/libRblas.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

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

other attached packages:
[1] knitr_1.21         biomaRt_2.38.0     RColorBrewer_1.1-2
[4] dplyr_0.8.0.1      tidyr_0.8.2        ggplot2_3.1.0     

loaded via a namespace (and not attached):
 [1] progress_1.2.0       tidyselect_0.2.5     xfun_0.4            
 [4] reshape2_1.4.3       purrr_0.3.0          colorspace_1.4-0    
 [7] htmltools_0.3.6      stats4_3.5.1         yaml_2.2.0          
[10] blob_1.1.1           XML_3.98-1.12        rlang_0.3.1         
[13] pillar_1.3.1         glue_1.3.0           withr_2.1.2         
[16] DBI_1.0.0            BiocGenerics_0.28.0  bit64_0.9-7         
[19] plyr_1.8.4           stringr_1.4.0        munsell_0.5.0       
[22] gtable_0.2.0         workflowr_1.2.0      flashier_0.1.0      
[25] evaluate_0.12        memoise_1.1.0        labeling_0.3        
[28] Biobase_2.42.0       IRanges_2.16.0       curl_3.3            
[31] parallel_3.5.1       AnnotationDbi_1.44.0 highr_0.7           
[34] Rcpp_1.0.0           scales_1.0.0         backports_1.1.3     
[37] S4Vectors_0.20.1     fs_1.2.6             bit_1.1-14          
[40] hms_0.4.2            digest_0.6.18        stringi_1.2.4       
[43] cowplot_0.9.4        grid_3.5.1           rprojroot_1.3-2     
[46] tools_3.5.1          bitops_1.0-6         magrittr_1.5        
[49] lazyeval_0.2.1       RCurl_1.95-4.11      tibble_2.0.1        
[52] RSQLite_2.1.1        crayon_1.3.4         whisker_0.3-2       
[55] pkgconfig_2.0.2      prettyunits_1.0.2    assertthat_0.2.0    
[58] rmarkdown_1.11       httr_1.4.0           R6_2.4.0            
[61] git2r_0.23.0         compiler_3.5.1