Package Installation
To install this package, start R (version “4.3”) and run the following commands:
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("Summix")Summix2 Functionalities
Summix2 has three primary methods to detect and leverage substructure in genetic summary data. The summix() function is used to detect substructure proportions given AFs for the observed and reference groups, summix_local() is the used to detect local substructure proportions, and adjAF() adjusts AFs for an observed sample to match the substructure proportions of a target sample or individual.
Example of Summix2 Workflow
Before diving into running Summix2 functions and interpreting output, let’s start by considering a sample of African (AFR) individuals we are interested in studying. We see in the pie chart below that this sample has 100% African-like genetic substructure.
Lets assume we need access to allele frequencies (AFs) from an appropriate control group (a group with 100% African genetic substructure) to complete the study. We can access publicly available AFs from the Genome Aggregation Database (gnomAD). However, gnomAD only offers AFs across 20,744 African/African American (AFR/AFRAM) individuals. In order for these AFs to be appropriate controls for a study of African individuals, the population substructure within gnomAD’s AFR/AFRAM AFs needs to be accounted for.
Here is where Summix2 comes in.
We begin by using summix() to capture the population substructure in the AFR/AFRAM AFs. We prepare a data frame containing the observed gnomAD AFR/AFRAM AFs, and reference AFs from homogenous groups. Here, we use AFs from homogenous continental-level reference groups- sourced from the Human Genome Diversity Project and 1000 Genomes Project data released with gnomAD v3.1.2.
library(Summix)
print(head(ancestryData))
#> POS REF ALT CHROM reference_AF_afr reference_AF_eas reference_AF_eur
#> 1 31652001 T A chr22 0.040925268 0.0000000 0.000000000
#> 2 34509945 C G chr22 0.217971527 0.0000000 0.000000000
#> 3 34636589 CAA C chr22 0.181117576 0.0000000 0.001149425
#> 4 38889885 A AAG chr22 0.007117446 0.0000000 0.000000000
#> 5 49160931 G T chr22 0.064056997 0.0000000 0.000000000
#> 6 17604199 C CAGGA chr22 0.219750879 0.1654624 0.070476185
#> reference_AF_iam reference_AF_sas gnomad_AF_afr
#> 1 0.00000000 0.00000000 0.04171490
#> 2 0.00000000 0.00000000 0.18774500
#> 3 0.00000000 0.00000000 0.15198300
#> 4 0.00000000 0.00000000 0.00422064
#> 5 0.00000000 0.00000000 0.05445710
#> 6 0.09523803 0.07459678 0.19200400
Next, we apply summix() to the data frame with the observed group set as gnomAD’s AFR/AFRAM AF vector (observed=“gnomad_AF_afr”) the reference groupings set as all homogenous continental-level reference groups (reference=c(“reference_AF_afr”, “reference_AF_eas”, “reference_AF_eur”, “reference_AF_iam”, “reference_AF_sas”).
summix(data = ancestryData,
reference=c("reference_AF_afr",
"reference_AF_eas",
"reference_AF_eur",
"reference_AF_iam",
"reference_AF_sas"),
observed="gnomad_AF_afr")
#> goodness.of.fit iterations time filtered reference_AF_afr
#> 1 0.4853597 20 0.8998544 secs 0 0.812142
#> reference_AF_eur reference_AF_iam
#> 1 0.169953 0.017905
In this output we see summix() estimates the observed gnomAD AFR/AFRAM AFs to contain approximately 81% African-like, 17% EUR-like, and 1.8% IAM-like mixing proportions.
Now that we have captured the substructure within the genetic summary data, we use Summix2’s adjAF() to adjust the substructure to match that of the 100% African target sample.
When adjusting the genetic substructure, we use only the reference groups with non-zero mixing proportions in the observed and target samples, and ensure the order of mixing proportions in pi.target and pi.observed are in the same order as the reference groups.
We also include the sample sizes for each of the reference groups (N_reference = c(704,741, 47)) and the observed group (N_observed = 20744) to get the effective sample size calculation for the adjusted AF vector.
adjusted_data<-adjAF(data = ancestryData,
reference = c("reference_AF_afr", "reference_AF_eur", "reference_AF_iam"),
observed = "gnomad_AF_afr",
pi.target = c(1, 0, 0),
pi.observed = c(0.812142, 0.169953, 0.017905),
adj_method = 'average',
N_reference = c(704,741, 47),
N_observed = 20744,
filter = TRUE)
#>
#>
#> [1] "Note: In this AF adjustment, 0 SNPs (with adjusted AF > -.005 & < 0) were rounded to 0. 0 SNPs (with adjusted AF > 1) were rounded to 1, and 0 SNPs (with adjusted AF <= -.005) were removed from the final results."
#>
#> [1] $pi
#> ref.group pi.observed pi.target
#> 1 reference_AF_afr 0.812142 1
#> 2 reference_AF_eur 0.169953 0
#> 3 reference_AF_iam 0.017905 0
#>
#> [1] $observed.data
#> [1] "observed AF data to update: 'gnomad_AF_afr'"
#>
#> [1] $Nsnps
#> [1] 1000
#>
#>
#> [1] $effective.sample.size
#> [1] 17551
#>
#>
#> [1] "use $adjusted.AF$adjustedAF to see adjusted AF data"
#>
#>
#> [1] "Note: The accuracy of the AF adjustment is likely lower for rare variants (< .5%)."We see that no SNPs were removed from the AF adjusted data frame (only SNPs with AFs that were less than -.005 are removed). Importantly, the effective sample size of the adjusted AF is 17,551 - which can be important for downstream analyses.
We can take a look at the adjusted AF data frame, and the appended adjusted AF vector (adjustedAF).
print(adjusted_data$adjusted.AF[1:5,])
#> POS REF ALT CHROM reference_AF_afr reference_AF_eas reference_AF_eur
#> 1 31652001 T A chr22 0.040925268 0 0.000000000
#> 2 34509945 C G chr22 0.217971527 0 0.000000000
#> 3 34636589 CAA C chr22 0.181117576 0 0.001149425
#> 4 38889885 A AAG chr22 0.007117446 0 0.000000000
#> 5 49160931 G T chr22 0.064056997 0 0.000000000
#> reference_AF_iam reference_AF_sas gnomad_AF_afr adjustedAF
#> 1 0 0 0.04171490 0.044404861
#> 2 0 0 0.18774500 0.222371894
#> 3 0 0 0.15198300 0.183044358
#> 4 0 0 0.00422064 0.006477272
#> 5 0 0 0.05445710 0.065055887
Using Summix2, we can harmonize the genetic substructure across multiple data sets for secondary analyses. 