-
Notifications
You must be signed in to change notification settings - Fork 64
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
SNP weight values #357
Comments
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
I am working with adegenet and I am wondering about if there is a way to calculate SNP weight values. I am working with covariance matrices produced by ANGSD and PCangsd. ANGSD outputs a covariance matrix without any SNP weights. PCangsd has the ability to output SNP weights. It would be really cool if I could output the SNPs with the highest weights for each discriminant function. The code I am currently running is posted at the bottom.
Thanks so much,
Tanner
library(adegenet)
library(tibble)
set.seed(123)
Load the covariance matrix. Don't forget to change the file path if you have downloaded the data to your own computer.
cov = as.matrix(read.table("/Users/tannervanorden/Dropbox/Cutthroat_Research/PCA_Graphs/Original_updated_errorProof.covMat"), header = F)
We will also add a column with population assignments
Hendrys <- rep("Hendry's Creek", 42)
PineR <- rep('Pine & Ridge Creeks', 29)
Mill <- rep("Mill Creek", 18)
Trout = rep("Trout Creek", 4)
Mountain_Dell = rep("Mountain Dell", 2)
Birch = rep("Birch Creek", 3)
Willard = rep("Willard Creek", 30)
population_labels <- c(Hendrys, PineR, Mill, Trout, Mountain_Dell, Birch, Willard)
Perform PCA
mme.pca <- eigen(cov)
Extract eigenvectors and combine with population assignments
eigenvectors = mme.pca$vectors
pca_data <- as_tibble(cbind(population_labels, data.frame(eigenvectors)))
dapc_pops <- dapc(pca_data[, -1], grp = population_labels, perc.pca = 100, n.da = 4)
a_score_species <- optim.a.score(dapc_pops, k = 7) ##Optimum number of principle components for 4 discriminant functions is 75 for this data
Perform DAPC on the PCA results
dapc_result <- dapc(pca_data[, -1], grp = population_labels, n.pca = 75, n.da = 4)
The text was updated successfully, but these errors were encountered: