### R code from vignette source 'LEA.Rnw' ################################################### ### code chunk number 1: LEA.Rnw:48-53 ################################################### # creation of a directory for the LEA analyses dir.create("LEA_analyses") # set this new directory as the working directory setwd("LEA_analyses") ################################################### ### code chunk number 2: LEA.Rnw:77-88 ################################################### library(LEA) # Creation of the genotypic file "genotypes.lfmm" # with 400 SNPs for 50 individuals. data("tutorial") # in the lfmm format write.lfmm(tutorial.R, "genotypes.lfmm") # in the geno format write.geno(tutorial.R, "genotypes.geno") # creation of the environment file, gradient.env. # It contains 1 environmental variable for 50 individuals. write.env(tutorial.C, "gradients.env") ################################################### ### code chunk number 3: LEA.Rnw:99-110 ################################################### # run of pca # Available options, K (the number of PCs calculated), # center and scale. # Creation of genotypes.pcaProject - the pcaProject object. # a directory genotypes.pca containing: # Create files: genotypes.eigenvalues - eigenvalues, # genotypes.eigenvectors - eigenvectors, # genotypes.sdev - standard deviations, # genotypes.projections - projections, # Create a pcaProject object: pc. pc = pca("genotypes.lfmm", scale = TRUE) ################################################### ### code chunk number 4: LEA.Rnw:115-118 ################################################### # Perfom Tracy-Widom tests on all eigenvalues. # create file: tuto.tracyWidom - tracy-widom test information. tw = tracy.widom(pc) ################################################### ### code chunk number 5: LEA.Rnw:121-123 ################################################### # display the p-values for the Tracy-Widom tests. tw$pvalues[1:10] ################################################### ### code chunk number 6: LEA.Rnw:128-130 ################################################### # plot the percentage of variance explained by each component plot(tw$percentage) ################################################### ### code chunk number 7: LEA.Rnw:142-150 ################################################### # main options, K: (the number of ancestral populations), # entropy: calculate the cross-entropy criterion, # CPU: the number of CPUs. # Runs with K between 1 and 10 with cross-entropy and 10 repetitions. project = NULL project = snmf("genotypes.geno", K=1:10, entropy = TRUE, repetitions = 10, project = "new") ################################################### ### code chunk number 8: LEA.Rnw:157-159 ################################################### # plot cross-entropy criterion of all runs of the project plot(project, lwd = 5, col = "red", pch=1) ################################################### ### code chunk number 9: LEA.Rnw:165-170 ################################################### # get the cross-entropy of each run for K = 4 ce = cross.entropy(project, K = 4) # select the run with the lowest cross-entropy best = which.min(ce) ################################################### ### code chunk number 10: LEA.Rnw:189-199 ################################################### # main options, K: (the number of latent factors), # CPU: the number of CPUs. # Runs with K = 6 and 5 repetitions. # The runs are composed of 6000 iterations including 3000 iterations # for burnin. # around 20 seconds per run. project = NULL project = lfmm("genotypes.lfmm", "gradients.env", K = 6, repetitions = 5, project = "new") ################################################### ### code chunk number 11: LEA.Rnw:215-221 ################################################### # calculate adjusted p-values res = adjusted.pvalues(project, K = 6) cp.values = res$p.values gif = res$genomic.inflation.factor ################################################### ### code chunk number 12: LEA.Rnw:229-242 ################################################### for (alpha in c(.05,.1,.15,.2)) { # expected FDR print(paste("expected FDR:", alpha)) L = length(cp.values) # return a list of candidates with an expected FDR of alpha. w = which(sort(cp.values) < alpha * (1:L) / L) candidates = order(cp.values)[w] # estimated FDR and True Positif estimated.FDR = length(which(candidates <= 350))/length(candidates) estimated.TP = length(which(candidates > 350))/50 print(paste("FDR:", estimated.FDR, "True Positive:", estimated.TP)) } ################################################### ### code chunk number 13: LEA.Rnw:246-249 ################################################### # Copy of the pdf figures in the previous directory # for the creation of the vignette. file.copy(list.files(".", pattern = ".pdf"), "..")