Last updated: 2019-03-05

Checks: 5 1

Knit directory: drift-workflow/analysis/

This reproducible R Markdown analysis was created with workflowr (version 1.2.0). The Report tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.


The R Markdown file has unstaged changes. To know which version of the R Markdown file created these results, you’ll want to first commit it to the Git repo. If you’re still working on the analysis, you can ignore this warning. When you’re finished, you can run wflow_publish to commit the R Markdown file and build the HTML.

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

The command set.seed(20190211) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:


Ignored files:
    Ignored:    .Rhistory
    Ignored:    analysis/.Rhistory
    Ignored:    analysis/figure/
    Ignored:    analysis/flash_cache/
    Ignored:    data.tar.gz
    Ignored:    data/datasets/
    Ignored:    data/raw/
    Ignored:    output.tar.gz
    Ignored:    output/admixture/
    Ignored:    output/admixture_benchmark/
    Ignored:    output/benchmark/
    Ignored:    output/flash_greedy/
    Ignored:    output/log/
    Ignored:    output/pca/
    Ignored:    output/sim/

Unstaged changes:
    Modified:   analysis/hoa_global.Rmd

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.


These are the previous versions of the R Markdown and HTML files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view them.

File Version Author Date Message
html fb4c571 jhmarcus 2019-03-05 fixed md formatting;
Rmd 9a69c08 jhmarcus 2019-03-04 updated data hoa
html 9a69c08 jhmarcus 2019-03-04 updated data hoa
html c203e00 jhmarcus 2019-03-03 Build site.
Rmd 7392ef2 jhmarcus 2019-03-03 wflow_publish("analysis/*.Rmd")
Rmd 90a0a02 jhmarcus 2019-02-28 updated to autosomes and added manhatten plot
html 90a0a02 jhmarcus 2019-02-28 updated to autosomes and added manhatten plot
Rmd d01c17c jhmarcus 2019-02-28 added chrom exploration
Rmd 5ee97ed jhmarcus 2019-02-27 updating manhatten plots
Rmd e1b4f85 jhmarcus 2019-02-25 added snakemake rule
html e1b4f85 jhmarcus 2019-02-25 added snakemake rule
html 38a461d jhmarcus 2019-02-24 built all
html b5aafbb jhmarcus 2019-02-24 Build site.
Rmd 63e7173 jhmarcus 2019-02-24 expanded upon mean variance
html 63e7173 jhmarcus 2019-02-24 expanded upon mean variance
Rmd 38b57c5 jhmarcus 2019-02-24 simplified greedy flash global analysis
html 38b57c5 jhmarcus 2019-02-24 simplified greedy flash global analysis
Rmd 403bc6b jhmarcus 2019-02-15 added hide code blocks
html 403bc6b jhmarcus 2019-02-15 added hide code blocks
Rmd b4749ac jhmarcus 2019-02-15 fixed some typos
html b4749ac jhmarcus 2019-02-15 fixed some typos
Rmd 7a2b6c7 jhmarcus 2019-02-15 added backfit
html 7a2b6c7 jhmarcus 2019-02-15 added backfit
Rmd f5ef1af jhmarcus 2019-02-15 added workflows for human origins datasets
html f5ef1af jhmarcus 2019-02-15 added workflows for human origins datasets
Rmd 4afc77e jhmarcus 2019-02-15 init hoa global analysis

This is an analysis of the full Human Origins dataset which includes 2068 sampled from around the world. I filtered out rare variants with global minor allele frequency less than 5%, removed any variants with a missingness level greater than 0.5%, and removed any SNPs on the sex chromosomes, resulting in 343758 SNPs … see Human Origins Array Global Data for details on the data pre-processing.

Imports

Lets import some needed packages:

library(ggplot2)
library(tidyr)
library(dplyr)
library(RColorBrewer)
library(biomaRt)
library(knitr)
source("../code/viz.R")

FLASH-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,
                                         ash.param=list(method='fdr'),
                                         prior.type=c('nonnegative', 'normal.mixture'),
                                         var.type=2,
                                         fixed.factors=flashier::ones.factor(1),
                                         output.lvl=4,
                                         verbose.lvl=2)
          # save the rds
          saveRDS(flash_fit, '{output.rds}')

Lets first read the greedy flashier fit

# read the flash fit output by snakemake
flash_fit = readRDS("../output/flash_greedy/hoa_global/HumanOriginsPublic2068_auto_maf05_geno005.rds")

# read the snp meta data
bim_df = read.table("../data/datasets/hoa_global/HumanOriginsPublic2068_auto_maf05_geno005.bim", header=F)
colnames(bim_df) = c("chrom", "rsid", "cm", "pos", "a1", "a2")

# drift event loadings lfsr
l_lfsr_df = data.frame(flash_fit$lfsr[[1]])
colnames(l_lfsr_df) = 1:21

# drift event lfsr
delta_lfsr_df = data.frame(flash_fit$lfsr[[2]])
colnames(delta_lfsr_df) = 1:21
delta_lfsr_df$chrom = bim_df$chrom
delta_lfsr_df$pos = bim_df$pos
delta_lfsr_df$rsid = bim_df$rsid

# drift events
delta_df = as.data.frame(flash_fit$loadings$normalized.loadings[[2]])
colnames(delta_df) = 1:21
delta_df$chrom = bim_df$chrom
delta_df$pos = bim_df$pos
delta_df$rsid = bim_df$rsid

# number of drift events
K = ncol(flash_fit$loadings$normalized.loadings[[1]]) 

# number of individuals
n = nrow(flash_fit$loadings$normalized.loadings[[1]])

# number of SNPs
p = nrow(flash_fit$loadings$normalized.loadings[[2]])
print(K)
[1] 21
print(n)
[1] 2068
print(p)
[1] 343758

PVEs

Here is a plot of the proportion of variance (PVE) explained by each drift event:

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

Version Author Date
c203e00 jhmarcus 2019-03-03
90a0a02 jhmarcus 2019-02-28
38b57c5 jhmarcus 2019-02-24
print(flash_fit$pve)
 [1] 0.4698466181 0.0195673554 0.0212183760 0.0110132180 0.0026319565
 [6] 0.0047422304 0.0016326107 0.0015178954 0.0012378907 0.0008897922
[11] 0.0007602681 0.0006964420 0.0005662946 0.0002224860 0.0002020599
[16] 0.0003186703 0.0002843464 0.0001670124 0.0002173951 0.0001456609
[21] 0.0001285236

It looks like the PVE drops off at around 14 or so?

Fitted Mean and Variance

I setup the flashier run so it estimates a SNP specific precision term. Here is a histogram of fitted variances:

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

Lets now look the the fitted means:

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

The mean term should roughly be interpreted as the mean minor allele frequency at the SNP:

maf_df = read.table("../data/datasets/hoa_global/HumanOriginsPublic2068_auto_maf05_geno005.frq", 
                    stringsAsFactors=F, 
                    header=T)

maf_df$mu = mu

p_mu_maf = ggplot(maf_df %>% sample_n(10000), aes(x=MAF, y=mu)) + 
           geom_point() +
           xlab("MAF") + 
           ylab("Estimated Mean") + 
           theme_bw() +
           xlim(0.0, .5) +
           ylim(0.0, .5) +
           geom_abline(slope=1, intercept=0)

p_mu_maf

Its looks like the estimated mean is matching the MAF but is biased … it tends to be underestimating the MAF. 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("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

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). I’m not sure why there is a sharp cuttof at ~.35.

Drift Event Distributions

Lets now plot the distribution of drift events:

# 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_de = 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_de

We can see the lower PVE drift events tend to get sparser!

Drift Event Local False Sign Rates

Here histograms of the lfsrs for each drift event:

delta_lfsr_gath_df = delta_lfsr_df %>% 
                     gather(K, value, -chrom, -pos, -rsid) %>%
                     filter(K %in% paste0(2:21))

p_lfsr = ggplot(delta_lfsr_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() + 
         xlab("Local False Sign Rate")
p_lfsr
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Version Author Date
c203e00 jhmarcus 2019-03-03

It is interesting to see the lower PVE factors (around 13 onwards) shift to having many SNPs with high uncertainty of the sign of the drift event.

Drift Event Loadings (2-11)

Lets now take a look at the drift event loadings. First we setup the data:

# 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)

# all unique pop labels
pops = unique(l_df$Simple.Population.ID) 

# 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.344880e-11 0.009957689 0.02220615 3.206359e-07 2.177500e-07
2 0.02198997 3.092709e-11 0.009454017 0.02319302 5.477147e-07 1.330105e-07
3 0.02198997 2.875096e-11 0.006741585 0.02582431 1.623606e-07 3.339426e-07
4 0.02198997 2.939241e-11 0.009538174 0.02405254 1.030657e-07 4.949653e-07
5 0.02198997 2.955584e-11 0.010435541 0.02369975 1.748882e-07 8.876189e-07
6 0.02198997 3.338298e-11 0.010010116 0.02221653 1.824629e-07 3.673969e-07
            7            8            9         10           11         12
1 0.009714513 3.894804e-07 6.901002e-09 0.05180365 4.698266e-07 0.01754951
2 0.010167281 3.028628e-07 1.073222e-08 0.04694396 2.341897e-06 0.01181955
3 0.009951402 1.618942e-07 5.953407e-09 0.05151676 1.361510e-06 0.01284319
4 0.007536562 2.554642e-07 7.551714e-09 0.05886751 4.496300e-07 0.01684531
5 0.005643071 3.048893e-07 6.285200e-09 0.05447787 5.993484e-07 0.01786413
6 0.010378956 1.590768e-07 7.388866e-09 0.05291977 3.552595e-07 0.01761416
            13           14         15           16           17
1 2.003361e-06 2.147565e-07 0.06517419 0.0015261431 0.0016538921
2 4.600137e-03 1.131473e-07 0.04171043 0.0060590808 0.0049734659
3 8.890164e-06 1.486041e-07 0.06591569 0.0049143357 0.0004319419
4 1.135608e-04 9.539378e-08 0.06514533 0.0047344501 0.0002431952
5 6.119798e-06 1.062652e-07 0.07251071 0.0008047695 0.0010858248
6 4.957157e-06 1.422404e-07 0.05927887 0.0004251974 0.0000393936
            18           19           20           21            ID
1 9.852204e-06 1.654588e-06 1.900904e-06 6.540409e-07 Algerian43A22
2 3.790969e-06 1.970198e-06 7.016039e-06 1.386323e-06 Algerian43A21
3 2.049506e-03 1.821203e-06 7.141382e-06 5.548505e-07 Algerian43A34
4 6.631587e-06 1.336589e-06 2.237992e-05 3.818413e-07 Algerian43A13
5 6.361275e-06 2.150546e-06 4.261356e-04 4.396811e-07 Algerian43A24
6 2.049671e-05 2.222535e-06 2.644428e-06 4.196996e-07 Algerian43A32
  Simple.Population.ID Verbose.Population.ID Region Country Latitude
1             Algerian              Algerian Africa Algeria     36.8
2             Algerian              Algerian Africa Algeria     36.8
3             Algerian              Algerian Africa Algeria     36.8
4             Algerian              Algerian Africa Algeria     36.8
5             Algerian              Algerian Africa Algeria     36.8
6             Algerian              Algerian Africa Algeria     36.8
  Longitude Samples Passed.QC Contributor
1         3       7         7 David Comas
2         3       7         7 David Comas
3         3       7         7 David Comas
4         3       7         7 David Comas
5         3       7         7 David Comas
6         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(gath_df=l_gath_df %>% 
                                   filter(Region == "Africa"), 
                                   colset="Set3",
                                   facet_levels=africa_pops,
                                   facet_grp="Simple.Population.ID", 
                                   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(gath_df=l_gath_df %>% 
                                    filter(Region == "America"), 
                                    colset="Set3",
                                    facet_levels=america_pops,
                                    facet_grp="Simple.Population.ID", 
                                    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(gath_df=l_gath_df %>% 
                                                 filter(Region == "CentralAsiaSiberia"), 
                                                 colset="Set3",
                                                 facet_levels=central_asia_siberia_pops,  
                                                 facet_grp="Simple.Population.ID",
                                                 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(gath_df=l_gath_df %>% 
                                      filter(Region == "EastAsia"), 
                                      colset="Set3",
                                      facet_levels=east_asia_pops,  
                                      facet_grp="Simple.Population.ID",
                                      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(gath_df=l_gath_df %>% 
                                      filter(Region == "SouthAsia"),
                                      colset="Set3",
                                      facet_levels=south_asia_pops, 
                                      facet_grp="Simple.Population.ID",
                                      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(gath_df=l_gath_df %>% 
                                         filter(Region == "WestEurasia"), 
                                         colset="Set3",
                                         facet_levels=west_eurasia_pops, 
                                         facet_grp="Simple.Population.ID",
                                         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(gath_df=l_gath_df %>% 
                                    filter(Region == "Oceania"), 
                                    colset="Set3",
                                    facet_levels=oceania_pops, 
                                    facet_grp="Simple.Population.ID",
                                    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

Drift Event Loadings (12-21)

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(gath_df=l_gath_df %>% 
                                   filter(Region == "Africa"), 
                                   colset="Set3",
                                   facet_levels=africa_pops,
                                   facet_grp="Simple.Population.ID", 
                                   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(gath_df=l_gath_df %>% 
                                    filter(Region == "America"), 
                                    colset="Set3",
                                    facet_levels=america_pops,
                                    facet_grp="Simple.Population.ID", 
                                    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(gath_df=l_gath_df %>% 
                                                 filter(Region == "CentralAsiaSiberia"), 
                                                 colset="Set3",
                                                 facet_levels=central_asia_siberia_pops,  
                                                 facet_grp="Simple.Population.ID",
                                                 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(gath_df=l_gath_df %>% 
                                      filter(Region == "EastAsia"), 
                                      colset="Set3",
                                      facet_levels=east_asia_pops,  
                                      facet_grp="Simple.Population.ID",
                                      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(gath_df=l_gath_df %>% 
                                      filter(Region == "SouthAsia"),
                                      colset="Set3",
                                      facet_levels=south_asia_pops, 
                                      facet_grp="Simple.Population.ID",
                                      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(gath_df=l_gath_df %>% 
                                         filter(Region == "WestEurasia"), 
                                         colset="Set3",
                                         facet_levels=west_eurasia_pops, 
                                         facet_grp="Simple.Population.ID",
                                         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(gath_df=l_gath_df %>% 
                                    filter(Region == "Oceania"), 
                                    colset="Set3",
                                    facet_levels=oceania_pops, 
                                    facet_grp="Simple.Population.ID",
                                    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

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.

Drift Event Loadings Local False Sign Rates

l_lfsr_gath_df = l_lfsr_df %>% 
                 gather(K, value) %>%
                 filter(K %in% paste0(2:21))

p_lfsr = ggplot(l_lfsr_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() + 
         xlab("Local False Sign Rate")
p_lfsr
`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Version Author Date
c203e00 jhmarcus 2019-03-03

Its very interesting to see that the lfsr seem very bi-model for the loadings for many of the drift events.

PCA

Here are proportions of variances explained by each PC:

pca_fit = readRDS("../output/pca/hoa_global/HumanOriginsPublic2068_auto_maf05_geno005.rds")
lamb2 = pca_fit$d^2

p_pve = qplot(x=1:length(lamb2), y=lamb2 / sum(lamb2)) + 
        xlab("PC") +
        ylab("PVE") +
        theme_bw()
        
p_pve

Lets plot the PC loadings

l_df = cbind(data.frame(pca_fit$v), meta_df)
colnames(l_df)[1:20] = 1:20
l_df = l_df %>% rowwise() %>% mutate(lab=substr(Simple.Population.ID, 1, 3))

# 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, 
                   -lab) 

pc_df = l_gath_df %>% 
        inner_join(l_gath_df, by="ID") %>% 
        filter((K.x == "1" & K.y == "2") | 
               (K.x == "3" & K.y == "4") |
               (K.x == "5" & K.y == "6") |
               (K.x == "7" & K.y == "8") |
               (K.x == "9" & K.y == "10") |
               (K.x == "11" & K.y == "12") |
               (K.x == "13" & K.y == "14") |
               (K.x == "15" & K.y == "16") |
               (K.x == "17" & K.y == "18") |
               (K.x == "19" & K.y == "20")) %>%
        mutate(Region=Region.x)

p_pca = ggplot(pc_df, aes(x=value.x, y=value.y, label=lab.x, color=Region)) +
        geom_text(size=2) +
        scale_color_brewer(palette = "Set1") + 
        facet_wrap(factor(K.x, levels=seq(1,20,2))~factor(K.y, levels=seq(2,20,2)), scales="free") + 
        theme_bw() + 
        xlab("PCx") +
        ylab("PCy")
        
p_pca

The classic horseshoe pattern is obtained that we often see for PCA on global populations … then different regions and populations have high loadings in subsequent PCs i.e. PC3~Ocenica, PC4~Americas, etc. I now see the advantage of the bar plot for visualizing lower PVE factors. Its “OK” if we only focus on PC1 / PC2 but the bar plot has a massive advantage when trying to compare many 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/HumanOriginsPublic2068_auto_maf05_geno005.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(gath_df=l_gath_df %>% 
                                   filter(Region == "Africa"), 
                                   colset="Set3",
                                   facet_levels=africa_pops,
                                   facet_grp="Simple.Population.ID", 
                                   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(gath_df=l_gath_df %>% 
                                    filter(Region == "America"), 
                                    colset="Set3",
                                    facet_levels=america_pops,
                                    facet_grp="Simple.Population.ID", 
                                    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(gath_df=l_gath_df %>% 
                                                 filter(Region == "CentralAsiaSiberia"), 
                                                 colset="Set3",
                                                 facet_levels=central_asia_siberia_pops,  
                                                 facet_grp="Simple.Population.ID",
                                                 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(gath_df=l_gath_df %>% 
                                      filter(Region == "EastAsia"), 
                                      colset="Set3",
                                      facet_levels=east_asia_pops,  
                                      facet_grp="Simple.Population.ID",
                                      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(gath_df=l_gath_df %>% 
                                      filter(Region == "SouthAsia"),
                                      colset="Set3",
                                      facet_levels=south_asia_pops, 
                                      facet_grp="Simple.Population.ID",
                                      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(gath_df=l_gath_df %>% 
                                         filter(Region == "WestEurasia"), 
                                         colset="Set3",
                                         facet_levels=west_eurasia_pops, 
                                         facet_grp="Simple.Population.ID",
                                         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(gath_df=l_gath_df %>% 
                                    filter(Region == "Oceania"), 
                                    colset="Set3",
                                    facet_levels=oceania_pops, 
                                    facet_grp="Simple.Population.ID",
                                    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?

Outlier SNPs

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:

manh_df = delta_lfsr_df %>% 
          group_by(chrom) %>% 
          summarise(chr_len=max(pos)) %>% 
          mutate(tot=cumsum(chr_len)-chr_len) %>%
          dplyr::select(-chr_len) %>%
          left_join(delta_lfsr_df, ., by=c("chrom"="chrom")) %>%
          arrange(chrom, pos) %>%
          mutate(BPcum=pos+tot)

manh_axis_df = manh_df %>% 
               group_by(chrom) %>% 
               summarize(center=(max(BPcum) + min(BPcum)) / 2)

manh_gath_df = manh_df %>% gather(K, value, -chrom, -pos, -rsid, -tot, -BPcum) %>% 
               filter(K %in% paste0(2:11)) %>% 
               filter(value < 1e-10) %>%
               filter(chrom %in% 1:22)
          
p = ggplot(manh_gath_df, aes(x=BPcum, y=-log10(value))) +
    geom_point(aes(color=as.factor(chrom)), alpha=.7, size=.5) +
    scale_color_manual(values=rep(c("grey", "orange"), 22)) +
    scale_x_continuous(label=manh_axis_df$chrom, breaks=manh_axis_df$center) +
    scale_y_continuous(expand=c(0, 0)) +     
    theme_bw() +
    theme(axis.title.x=element_blank(), axis.ticks.x=element_blank(),
          panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
    facet_grid(factor(K, levels=2:11)~., scales = "free") + 
    ylab("-log10(lfsr)") + 
    guides(alpha=F, color=F)
p
Warning in FUN(X[[i]], ...): NaNs produced
Warning: Removed 5767 rows containing missing values (geom_point).

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:

# SNP bioMart database
grch37_snp = useMart(biomart="ENSEMBL_MART_SNP", 
                     host="grch37.ensembl.org", 
                     path="/biomart/martservice",
                     dataset="hsapiens_snp")

# GENE bioMart database
grch37 = useMart(biomart="ENSEMBL_MART_ENSEMBL", 
                 host="grch37.ensembl.org", 
                 path="/biomart/martservice", 
                 dataset="hsapiens_gene_ensembl")

# top rank SNPs per each drift event
delta_tophit_df = delta_df %>%
                  gather(K, value, -chrom, -pos, -rsid) %>%
                  filter(K %in% paste0(2:21)) %>% 
                  group_by(K) %>%
                  top_n(5, value)

# SNP meta data
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

# GENE meta data
table2 = getBM(attributes = c("ensembl_gene_id", "external_gene_name", "description"),
               filters = "ensembl_gene_id", 
               values =  table1$ensembl_gene_stable_id, 
               mart = grch37)

# annotation data
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
# unique genes
print(unique(anno_df$external_gene_name))
 [1] NA              "ANXA7"         "NAV2"          "SLC22A10"     
 [5] "APPL2"         "LINC00523"     "HERC2"         "RP11-109D20.1"
 [9] "SORD"          "FBXL22"        "USP3-AS1"      "PRKCB"        
[13] "EPN2"          "CACNB1"        "RP11-515E23.1" "RP11-1055B8.4"
[17] "ARHGAP28"      "KIF3C"         "EXOC6B"        "IL1R1"        
[21] "EDAR"          "CCNT2"         "UBXN4"         "PLA2R1"       
[25] "TLK1"          "KIF16B"        "LINC00478"     "NHP2L1"       
[29] "MEI1"          "KAT2B"         "FRMD4B"        "ROBO2"        
[33] "RP11-103J17.1" "TBC1D1"        "RP11-696N14.1" "ADH4"         
[37] "ADH1B"         "SLC45A2"       "PDE4D"         "STK32A"       
[41] "DPYSL3"        "RP11-267A15.1" "RP1-179N16.6"  "FOXP4"        
[45] "AKAP12"        "GS1-124K5.12"  "TMEM248"       "MAGI2"        
[49] "CYP7B1"        "RP11-1D12.2"   "SMARCA2"       "GLIS3"        
[53] "RCL1"         
# formatted table
d = anno_df %>% 
    distinct(external_gene_name, .keep_all = T) %>% 
    dplyr::select(external_gene_name, K, rsid, consequence_type_tv) %>% 
    arrange(consequence_type_tv)
kable(d)
external_gene_name K rsid consequence_type_tv
NA 9 rs4382718
ADH1B 7 rs1042026 3_prime_UTR_variant
ANXA7 11 rs56967412 intron_variant
NAV2 3 rs10741783 intron_variant
HERC2 12 rs12913832 intron_variant
SORD 5 rs2470686 intron_variant
FBXL22 2 rs4247077 intron_variant
PRKCB 13 rs4787651 intron_variant
EPN2 4 rs11652784 intron_variant
CACNB1 15 rs16531 intron_variant
ARHGAP28 19 rs11663962 intron_variant
KIF3C 2 rs11126303 intron_variant
EXOC6B 3 rs7420667 intron_variant
IL1R1 13 rs3755295 intron_variant
EDAR 2 rs922452 intron_variant
UBXN4 10 rs6732236 intron_variant
PLA2R1 8 rs17830904 intron_variant
TLK1 11 rs78813632 intron_variant
KIF16B 20 rs6044020 intron_variant
NHP2L1 6 rs5758406 intron_variant
MEI1 6 rs126092 intron_variant
KAT2B 14 rs2929402 intron_variant
FRMD4B 13 rs34266487 intron_variant
ROBO2 17 rs77908466 intron_variant
TBC1D1 11 rs62298552 intron_variant
SLC45A2 4 rs35397 intron_variant
PDE4D 5 rs173591 intron_variant
FOXP4 15 rs2148342 intron_variant
AKAP12 19 rs17081009 intron_variant
TMEM248 8 rs10234322 intron_variant
MAGI2 14 rs2364341 intron_variant
CYP7B1 17 rs4236936 intron_variant
SMARCA2 20 rs7035991 intron_variant
GLIS3 14 rs10974315 intron_variant
RCL1 3 rs172447 intron_variant
CCNT2 10 rs41521049 NMD_transcript_variant
RP11-1055B8.4 5 rs3744146 non_coding_transcript_exon_variant
SLC22A10 13 rs12804561 non_coding_transcript_variant
LINC00523 4 rs12881700 non_coding_transcript_variant
RP11-109D20.1 5 rs2470686 non_coding_transcript_variant
USP3-AS1 2 rs4247077 non_coding_transcript_variant
RP11-515E23.1 5 rs9908046 non_coding_transcript_variant
LINC00478 3 rs2823880 non_coding_transcript_variant
RP11-103J17.1 9 rs11733992 non_coding_transcript_variant
RP11-696N14.1 7 rs17028609 non_coding_transcript_variant
ADH4 7 rs17028609 non_coding_transcript_variant
STK32A 16 rs114085738 non_coding_transcript_variant
RP11-267A15.1 19 rs2936935 non_coding_transcript_variant
RP1-179N16.6 9 rs10447375 non_coding_transcript_variant
GS1-124K5.12 8 rs3757419 non_coding_transcript_variant
RP11-1D12.2 17 rs6472154 non_coding_transcript_variant
APPL2 20 rs12303948 synonymous_variant
DPYSL3 16 rs17106709 synonymous_variant
We can see many of these top outliers are in genes (but we have to consider the array design to know if this is unusual). Its cool to see a couple very interesting genes pop up including HERC2 (eye color) and SLC45A2 (skin color), EDAR (hair), ADH4 (alcohol dehydrogenase), ADH1B (alcohol dehydrogenase), which have been previously studied for their selection signatures. Its also fun to look at some of the rsids on https://popgen.uchicago.edu/ggv/ to get a sense of what kind of allele frequency distributions define each factor.

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