Performing a Genetic Mixture Analysis
This is done with the infer_mixture function. In the example data chinook_mix our data consist of fish caught in three different fisheries, rec1, rec2, and rec3 as denoted in the collection column. Each of those collections is treated as a separate sample, getting its own mixing proportion estimate. This is how it is run with the default options:
mix_est <- infer_mixture(reference = chinook,
mixture = chinook_mix,
gen_start_col = 5)## Collating data; compiling reference allele frequencies, etc. time: 1.95 seconds
## Computing reference locus specific means and variances for computing mixture z-scores time: 0.29 seconds
## Working on mixture collection: rec2 with 772 individuals
## calculating log-likelihoods of the mixture individuals. time: 0.12 seconds
## performing 2000 total sweeps, 100 of which are burn-in and will not be used in computing averages in method "MCMC" time: 0.64 seconds
## tidying output into a tibble. time: 0.14 seconds
## Working on mixture collection: rec1 with 743 individuals
## calculating log-likelihoods of the mixture individuals. time: 0.11 seconds
## performing 2000 total sweeps, 100 of which are burn-in and will not be used in computing averages in method "MCMC" time: 0.61 seconds
## tidying output into a tibble. time: 0.13 seconds
## Working on mixture collection: rec3 with 741 individuals
## calculating log-likelihoods of the mixture individuals. time: 0.11 seconds
## performing 2000 total sweeps, 100 of which are burn-in and will not be used in computing averages in method "MCMC" time: 0.61 seconds
## tidying output into a tibble. time: 0.14 secondsThe result comes back as a list of four tidy data frames:
mixing_proportions: the mixing proportions. The columnpiholds the estimated mixing proportion for each collection.indiv_posteriors: this holds, for each individual, the posterior means of group membership in each collection. ColumnPofZholds those values. Columnlog_likelihoodholds the log of the probability of the individuals genotype given it is from the collection. Also included aren_non_miss_lociandn_miss_lociwhich are the number of observed loci and the number of missing loci at the individual. A list columnmissing_locicontains vectors with the indices (and the names) of the loci that are missing in that individual. It also includes a columnz_scorewhich can be used to diagnose fish that don’t belong to any samples in the reference data base (see below).mix_prop_traces:MCMC traces of the mixing proportions for each collection. You will use these if you want to make density estimates of the posterior distribution of the mixing proportions or if you want to compute credible intervals.bootstrapped_proportions: This is NULL in the above example, but if we had chosenmethod = "PB"then this would be a tibble of bootstrap-corrected reporting unit mixing proportions.
These data frames look like this:
lapply(mix_est, head)## $mixing_proportions
## # A tibble: 6 × 4
## mixture_collection repunit collection pi
## <chr> <chr> <chr> <dbl>
## 1 rec2 CentralValleyfa Feather_H_sp 0.0759
## 2 rec2 CentralValleysp Butte_Cr_Sp 0.0000338
## 3 rec2 CentralValleysp Mill_Cr_sp 0.0000383
## 4 rec2 CentralValleysp Deer_Cr_sp 0.0000583
## 5 rec2 CentralValleysp UpperSacramento_R_sp 0.000470
## 6 rec2 CentralValleyfa Feather_H_fa 0.157
##
## $indiv_posteriors
## # A tibble: 6 × 10
## mixture_collection indiv repunit collection PofZ log_likelihood z_score
## <chr> <chr> <chr> <chr> <dbl> <dbl> <dbl>
## 1 rec2 T124711 Central… Feather_H… 1.75e-28 -137. -13.1
## 2 rec2 T124711 Central… Feather_H… 1.04e-27 -136. -12.6
## 3 rec2 T124711 Central… Butte_Cr_… 1.54e-24 -130. -10.5
## 4 rec2 T124711 Central… Mill_Cr_fa 9.50e-31 -135. -11.8
## 5 rec2 T124711 Central… Deer_Cr_fa 3.39e-28 -134. -11.6
## 6 rec2 T124711 Central… Mokelumne… 1.84e-27 -134. -12.3
## # … with 3 more variables: n_non_miss_loci <int>, n_miss_loci <int>,
## # missing_loci <list>
##
## $mix_prop_traces
## # A tibble: 6 × 5
## mixture_collection sweep repunit collection pi
## <chr> <int> <chr> <chr> <dbl>
## 1 rec2 0 CentralValleyfa Feather_H_sp 0.0145
## 2 rec2 0 CentralValleysp Butte_Cr_Sp 0.0145
## 3 rec2 0 CentralValleysp Mill_Cr_sp 0.0145
## 4 rec2 0 CentralValleysp Deer_Cr_sp 0.0145
## 5 rec2 0 CentralValleysp UpperSacramento_R_sp 0.0145
## 6 rec2 0 CentralValleyfa Feather_H_fa 0.0145
##
## $bootstrapped_proportions
## # A tibble: 0 × 0Setting the prior for the mixing proportions
In some cases there might be a reason to explicitly set the parameters of the Dirichlet prior on the mixing proportions of the collections. For a contrived example, we could imagine that we wanted a Dirichlet prior with all parameters equal to 1/(# of collections), except for the parameters for all the Central Valley Fall Run populations, to which we would like to assign Dirichlet parameters of 2. That can be accomplished with the pi_prior argument to the infer_mixture() function, which will let you pass in a tibble in which one column named “collection” gives the collection, and the other column, named “pi_param” gives the desired parameter.
Here we construct that kind of input:
prior_tibble <- chinook %>%
count(repunit, collection) %>%
filter(repunit == "CentralValleyfa") %>%
select(collection) %>%
mutate(pi_param = 2)
# see what it looks like:
prior_tibble## # A tibble: 8 × 2
## collection pi_param
## <chr> <dbl>
## 1 Battle_Cr 2
## 2 Butte_Cr_fa 2
## 3 Deer_Cr_fa 2
## 4 Feather_H_fa 2
## 5 Feather_H_sp 2
## 6 Mill_Cr_fa 2
## 7 Mokelumne_R_fa 2
## 8 Sacramento_R_lf 2Then we can run that in infer_mixture():
set.seed(12)
mix_est_with_prior <- infer_mixture(reference = chinook,
mixture = chinook_mix,
gen_start_col = 5,
pi_prior = prior_tibble)## Collating data; compiling reference allele frequencies, etc. time: 1.55 seconds
## Computing reference locus specific means and variances for computing mixture z-scores time: 0.25 seconds
## Working on mixture collection: rec2 with 772 individuals
## Joining, by = "collection"
## calculating log-likelihoods of the mixture individuals. time: 0.12 seconds
## performing 2000 total sweeps, 100 of which are burn-in and will not be used in computing averages in method "MCMC" time: 0.63 seconds
## tidying output into a tibble. time: 0.14 seconds
## Working on mixture collection: rec1 with 743 individuals
## Joining, by = "collection"
## calculating log-likelihoods of the mixture individuals. time: 0.11 seconds
## performing 2000 total sweeps, 100 of which are burn-in and will not be used in computing averages in method "MCMC" time: 0.61 seconds
## tidying output into a tibble. time: 0.13 seconds
## Working on mixture collection: rec3 with 741 individuals
## Joining, by = "collection"
## calculating log-likelihoods of the mixture individuals. time: 0.11 seconds
## performing 2000 total sweeps, 100 of which are burn-in and will not be used in computing averages in method "MCMC" time: 0.60 seconds
## tidying output into a tibble. time: 0.14 secondsand now, for fun, we can compare the results for the mixing proportions of different collections there with and without the prior for the mixture collection rec1:
comp_mix_ests <- list(
`pi (default prior)` = mix_est$mixing_proportions,
`pi (cv fall gets 2s prior)` = mix_est_with_prior$mixing_proportions
) %>%
bind_rows(.id = "prior_type") %>%
filter(mixture_collection == "rec1") %>%
select(prior_type, repunit, collection, pi) %>%
spread(prior_type, pi) %>%
mutate(coll_group = ifelse(repunit == "CentralValleyfa", "CV_fall", "Not_CV_fall"))
ggplot(comp_mix_ests,
aes(x = `pi (default prior)`,
y = `pi (cv fall gets 2s prior)`,
colour = coll_group
)) +
geom_point() +
geom_abline(slope = 1, intercept = 0, linetype = "dashed")
Yep, slightly different than before. Let’s look at the sums of everything:
comp_mix_ests %>%
group_by(coll_group) %>%
summarise(with_explicit_prior = sum(`pi (cv fall gets 2s prior)`),
with_default_prior = sum(`pi (default prior)`))## # A tibble: 2 × 3
## coll_group with_explicit_prior with_default_prior
## <chr> <dbl> <dbl>
## 1 CV_fall 0.824 0.821
## 2 Not_CV_fall 0.176 0.179We see that for the most part this change to the prior changed the distribution of fish into different collections within the Central Valley Fall reporting unit. This is not suprising—it is very hard to tell apart fish from those different collections. However, it did not greatly change the estimated proportion of the whole reporting unit. This also turns out to make sense if you consider the effect that the extra weight in the prior will have.
Aggregating collections into reporting units
This is a simple operation in the tidyverse:
# for mixing proportions
rep_mix_ests <- mix_est$mixing_proportions %>%
group_by(mixture_collection, repunit) %>%
summarise(repprop = sum(pi)) # adding mixing proportions over collections in the repunit## `summarise()` has grouped output by 'mixture_collection'. You can override using the `.groups` argument.# for individuals posteriors
rep_indiv_ests <- mix_est$indiv_posteriors %>%
group_by(mixture_collection, indiv, repunit) %>%
summarise(rep_pofz = sum(PofZ))## `summarise()` has grouped output by 'mixture_collection', 'indiv'. You can override using the `.groups` argument.Creating posterior density curves from the traces
The full MCMC output for the mixing proportions is available by default in the field $mix_prop_traces. This can be used to obtain an estimate of the posterior density of the mixing proportions.
Here we plot kernel density estimates for the 6 most abundant repunits from the rec1 fishery:
# find the top 6 most abundant:
top6 <- rep_mix_ests %>%
filter(mixture_collection == "rec1") %>%
arrange(desc(repprop)) %>%
slice(1:6)
# check how many MCMC sweeps were done:
nsweeps <- max(mix_est$mix_prop_traces$sweep)
# keep only rec1, then discard the first 200 sweeps as burn-in,
# and then aggregate over reporting units
# and then keep only the top6 from above
trace_subset <- mix_est$mix_prop_traces %>%
filter(mixture_collection == "rec1", sweep > 200) %>%
group_by(sweep, repunit) %>%
summarise(repprop = sum(pi)) %>%
filter(repunit %in% top6$repunit)## `summarise()` has grouped output by 'sweep'. You can override using the `.groups` argument.# now we can plot those:
ggplot(trace_subset, aes(x = repprop, colour = repunit)) +
geom_density()
Computing Credible Intervals from the Traces
Following on from the above example, we will use trace_subset to compute the equal-tail 95% credible intervals for the 6 most abundant reporting units in the rec1 fishery:
top6_cis <- trace_subset %>%
group_by(repunit) %>%
summarise(loCI = quantile(repprop, probs = 0.025),
hiCI = quantile(repprop, probs = 0.975))
top6_cis## # A tibble: 6 × 3
## repunit loCI hiCI
## <chr> <dbl> <dbl>
## 1 CaliforniaCoast 1.85e- 2 0.0439
## 2 CentralValleyfa 7.91e- 1 0.848
## 3 KlamathR 5.00e- 2 0.0877
## 4 NCaliforniaSOregonCoast 3.04e- 3 0.0182
## 5 RogueR 4.32e- 2 0.0809
## 6 UColumbiaRsufa 5.84e-21 0.0111Assessing whether individuals are not from any of the reference populations
Sometimes totally unexpected things happen. One situation we saw in the California Chinook fishery was samples coming to us that were actually coho salmon. Before we included coho salmon in the reference sample, these coho always assigned quite strongly to Alaska populations of Chinook, even though they don’t really look like Chinook at all.
In this case, it is useful to look at the raw log-likelihood values computed for the individual, rather than the scaled posterior probabilities. Because aberrantly low values of the genotype log-likelihood can indicate that there is something wrong. However, the raw likelihood that you get will depend on the number of missing loci, etc. rubias deals with this by computing a z-score for each fish. The Z-score is the Z statistic obtained from the fish’s log-likelihood (by subtracting from it the expected log-likelihood and dividing by the expected standard deviation). rubias’s implementation of the z-score accounts for the pattern of missing data, but it does this without all the simulation that gsi_sim does. This makes it much, much, faster—fast enough that we can compute it be default for every fish and every population.
Here, we will look at the z-score computed for each fish to the population with the highest posterior. (It is worth noting that you would never want to use the z-score to assign fish to different populations—it is only there to decide whether it looks like it might not have actually come from the population that it was assigned to, or any other population in the reference data set.)
# get the maximum-a-posteriori population for each individual
map_rows <- mix_est$indiv_posteriors %>%
group_by(indiv) %>%
top_n(1, PofZ) %>%
ungroup()If everything is kosher, then we expect that the z-scores we see will be roughly normally distributed. We can compare the distribution of z-scores we see with a bunch of simulated normal random variables.
normo <- tibble(z_score = rnorm(1e06))
ggplot(map_rows, aes(x = z_score)) +
geom_density(colour = "blue") +
geom_density(data = normo, colour = "black")
The normal density is in black and the distribution of our observed z_scores is in blue. They fit reasonably well, suggesting that there is not too much weird stuff going on overall. (That is good!)
The z_score statistic is most useful as a check for individuals. It is intended to be a quick way to identify aberrant individuals. If you see a z-score to the maximum-a-posteriori population for an individual in your mixture sample that is considerably less than z_scores you saw in the reference, then you might infer that the individual doesn’t actually fit any of the populations in the reference well.
Individuals of known origin in the mixture
Here I include a small, contrived example. We use the small_chinook data set so that it goes fast.
First, we analyze the data with no fish in the mixture of known collection
no_kc <- infer_mixture(small_chinook_ref, small_chinook_mix, gen_start_col = 5)## Collating data; compiling reference allele frequencies, etc. time: 0.19 seconds
## Computing reference locus specific means and variances for computing mixture z-scores time: 0.02 seconds
## Working on mixture collection: rec3 with 29 individuals
## calculating log-likelihoods of the mixture individuals. time: 0.00 seconds
## performing 2000 total sweeps, 100 of which are burn-in and will not be used in computing averages in method "MCMC" time: 0.03 seconds
## tidying output into a tibble. time: 0.03 seconds
## Working on mixture collection: rec1 with 36 individuals
## calculating log-likelihoods of the mixture individuals. time: 0.00 seconds
## performing 2000 total sweeps, 100 of which are burn-in and will not be used in computing averages in method "MCMC" time: 0.03 seconds
## tidying output into a tibble. time: 0.04 seconds
## Working on mixture collection: rec2 with 35 individuals
## calculating log-likelihoods of the mixture individuals. time: 0.00 seconds
## performing 2000 total sweeps, 100 of which are burn-in and will not be used in computing averages in method "MCMC" time: 0.03 seconds
## tidying output into a tibble. time: 0.03 secondsAnd look at the results for the mixing proportions:
no_kc$mixing_proportions %>%
arrange(mixture_collection, desc(pi))## # A tibble: 18 × 4
## mixture_collection repunit collection pi
## <chr> <chr> <chr> <dbl>
## 1 rec1 CentralValleyfa Feather_H_fa 0.849
## 2 rec1 CentralValleysp Deer_Cr_sp 0.0508
## 3 rec1 CaliforniaCoast Eel_R 0.0400
## 4 rec1 KlamathR Klamath_IGH_fa 0.0305
## 5 rec1 MidOregonCoast Umpqua_sp 0.0252
## 6 rec1 CentralValleywi Sacramento_H 0.00463
## 7 rec2 CentralValleyfa Feather_H_fa 0.809
## 8 rec2 KlamathR Klamath_IGH_fa 0.103
## 9 rec2 MidOregonCoast Umpqua_sp 0.0733
## 10 rec2 CentralValleysp Deer_Cr_sp 0.00550
## 11 rec2 CaliforniaCoast Eel_R 0.00479
## 12 rec2 CentralValleywi Sacramento_H 0.00474
## 13 rec3 CentralValleyfa Feather_H_fa 0.839
## 14 rec3 CaliforniaCoast Eel_R 0.0714
## 15 rec3 MidOregonCoast Umpqua_sp 0.0496
## 16 rec3 KlamathR Klamath_IGH_fa 0.0262
## 17 rec3 CentralValleysp Deer_Cr_sp 0.00875
## 18 rec3 CentralValleywi Sacramento_H 0.00542Now, we will do the same analysis, but pretend that we know that the first 8 of the 36 fish in fishery rec1 are from the Deer_Cr_sp collection.
First we have to add the known_collection column to the reference.
# make reference file that includes the known_collection column
kc_ref <- small_chinook_ref %>%
mutate(known_collection = collection) %>%
select(known_collection, everything())
# see what that looks like
kc_ref[1:10, 1:8]## # A tibble: 10 × 8
## known_collection sample_type repunit collection indiv Ots_94857.232
## <chr> <chr> <chr> <chr> <chr> <int>
## 1 Deer_Cr_sp reference CentralValleysp Deer_Cr_sp Deer_Cr_sp:0001 2
## 2 Deer_Cr_sp reference CentralValleysp Deer_Cr_sp Deer_Cr_sp:0002 2
## 3 Deer_Cr_sp reference CentralValleysp Deer_Cr_sp Deer_Cr_sp:0003 2
## 4 Deer_Cr_sp reference CentralValleysp Deer_Cr_sp Deer_Cr_sp:0004 4
## 5 Deer_Cr_sp reference CentralValleysp Deer_Cr_sp Deer_Cr_sp:0005 2
## 6 Deer_Cr_sp reference CentralValleysp Deer_Cr_sp Deer_Cr_sp:0006 4
## 7 Deer_Cr_sp reference CentralValleysp Deer_Cr_sp Deer_Cr_sp:0007 2
## 8 Deer_Cr_sp reference CentralValleysp Deer_Cr_sp Deer_Cr_sp:0008 4
## 9 Deer_Cr_sp reference CentralValleysp Deer_Cr_sp Deer_Cr_sp:0009 2
## 10 Deer_Cr_sp reference CentralValleysp Deer_Cr_sp Deer_Cr_sp:0010 2
## # … with 2 more variables: Ots_94857.232.1 <int>, Ots_102213.210 <int>Then we add the known collection column to the mixture. We start out making it all NAs, and then we change that to Deer_Cr_sp for 8 of the rec1 fish:
kc_mix <- small_chinook_mix %>%
mutate(known_collection = NA) %>%
select(known_collection, everything())
kc_mix$known_collection[kc_mix$collection == "rec1"][1:8] <- "Deer_Cr_sp"
# here is what that looks like now (dropping most of the genetic data columns)
kc_mix[1:20, 1:7]## # A tibble: 20 × 7
## known_collection sample_type repunit collection indiv Ots_94857.232
## <chr> <chr> <chr> <chr> <chr> <int>
## 1 <NA> mixture <NA> rec3 T125347 4
## 2 Deer_Cr_sp mixture <NA> rec1 T127759 4
## 3 <NA> mixture <NA> rec2 T124955 4
## 4 <NA> mixture <NA> rec2 T127564 2
## 5 <NA> mixture <NA> rec3 T127392 4
## 6 Deer_Cr_sp mixture <NA> rec1 T127414 4
## 7 <NA> mixture <NA> rec3 T124839 4
## 8 Deer_Cr_sp mixture <NA> rec1 T126414 2
## 9 <NA> mixture <NA> rec3 T125252 4
## 10 Deer_Cr_sp mixture <NA> rec1 T127765 4
## 11 Deer_Cr_sp mixture <NA> rec1 T127293 2
## 12 Deer_Cr_sp mixture <NA> rec1 T127577 2
## 13 <NA> mixture <NA> rec2 T126766 4
## 14 <NA> mixture <NA> rec3 T126494 4
## 15 <NA> mixture <NA> rec2 T125205 4
## 16 Deer_Cr_sp mixture <NA> rec1 T126584 4
## 17 <NA> mixture <NA> rec2 T124821 4
## 18 <NA> mixture <NA> rec3 T124905 2
## 19 Deer_Cr_sp mixture <NA> rec1 T124735 4
## 20 <NA> mixture <NA> rec3 T125624 2
## # … with 1 more variable: Ots_94857.232.1 <int>And now we can do the mixture analysis:
# note that the genetic data start in column 6 now
with_kc <- infer_mixture(kc_ref, kc_mix, 6)## Collating data; compiling reference allele frequencies, etc. time: 0.19 seconds
## Computing reference locus specific means and variances for computing mixture z-scores time: 0.02 seconds
## Working on mixture collection: rec3 with 29 individuals
## calculating log-likelihoods of the mixture individuals. time: 0.01 seconds
## performing 2000 total sweeps, 100 of which are burn-in and will not be used in computing averages in method "MCMC" time: 0.03 seconds
## tidying output into a tibble. time: 0.02 seconds
## Working on mixture collection: rec1 with 36 individuals
## calculating log-likelihoods of the mixture individuals. time: 0.01 seconds
## performing 2000 total sweeps, 100 of which are burn-in and will not be used in computing averages in method "MCMC" time: 0.03 seconds
## tidying output into a tibble. time: 0.03 seconds
## Working on mixture collection: rec2 with 35 individuals
## calculating log-likelihoods of the mixture individuals. time: 0.00 seconds
## performing 2000 total sweeps, 100 of which are burn-in and will not be used in computing averages in method "MCMC" time: 0.03 seconds
## tidying output into a tibble. time: 0.03 secondsAnd, when we look at the estimated proportions, we see that for rec1 they reflect the fact that 8 of those fish were singled out as known fish from Deer_Cr_sp:
with_kc$mixing_proportions %>%
arrange(mixture_collection, desc(pi))## # A tibble: 18 × 4
## mixture_collection repunit collection pi
## <chr> <chr> <chr> <dbl>
## 1 rec1 CentralValleyfa Feather_H_fa 0.546
## 2 rec1 CentralValleysp Deer_Cr_sp 0.355
## 3 rec1 CaliforniaCoast Eel_R 0.0411
## 4 rec1 KlamathR Klamath_IGH_fa 0.0318
## 5 rec1 MidOregonCoast Umpqua_sp 0.0220
## 6 rec1 CentralValleywi Sacramento_H 0.00438
## 7 rec2 CentralValleyfa Feather_H_fa 0.806
## 8 rec2 KlamathR Klamath_IGH_fa 0.104
## 9 rec2 MidOregonCoast Umpqua_sp 0.0732
## 10 rec2 CentralValleysp Deer_Cr_sp 0.00688
## 11 rec2 CaliforniaCoast Eel_R 0.00551
## 12 rec2 CentralValleywi Sacramento_H 0.00456
## 13 rec3 CentralValleyfa Feather_H_fa 0.835
## 14 rec3 CaliforniaCoast Eel_R 0.0732
## 15 rec3 MidOregonCoast Umpqua_sp 0.0494
## 16 rec3 KlamathR Klamath_IGH_fa 0.0281
## 17 rec3 CentralValleysp Deer_Cr_sp 0.00883
## 18 rec3 CentralValleywi Sacramento_H 0.00565The output from infer_mixture() in this case can be used just like it was before without known individuals in the baseline.
Fully Bayesian model (with updating of allele freqencies)
The default model in rubias is a conditional model in which inference is done with the baseline allele counts fixed. In a fully Bayesian version, fish from within the mixture that are allocated (on any particular step of the MCMC) to one of the reference samples have their alleles added to that reference sample, thus (one hopes) refining the estimate of allele frequencies in that sample. This is more computationally intensive, and, is done using parallel computation, by default running one thread for every core on your machine.
The basic way to invoke the fully Bayesian model is to use the infer_mixture function with the method option set to “BR”. For example:
full_model_results <- infer_mixture(
reference = chinook,
mixture = chinook_mix,
gen_start_col = 5,
method = "BR"
)More details about different options for working with the fully Bayesian model are available in the vignette about the fully Bayesian model.
