Small code changes relevant to vegan package

[BrennaF]

A vegan contributor warned me that "summary(rda.object)" would be broken in a future vegan package update; these small code changes will prevent vignette errors after these updates occur.

[hlapp]

So we're still getting this error:

Quitting from lines 119-143 (2016-01-26-SNP-selection.Rmd) 
Error in UseMethod("pcadapt") : 
  no applicable method for 'pcadapt' applied to an object of class "character"
Calls: <Anonymous> ... withCallingHandlers -> withVisible -> eval -> eval -> pcadapt

That's even though the hlapp/rpopgen container builds successfully and uses the latest pcadapt version from Github. Any idea what the issue might be?

[BrennaF]

Hi Hilmar! That is strange since the RDA vignette doesn't call pcadapt. Looks like the error is referring to the 2016-01-26-SNP-selection.Rmd vignette?

[hlapp]

Looks like the error is referring to the 2016-01-26-SNP-selection.Rmd vignette?

Ah yes, indeed. My mistake. @zkamvar (or perhaps @keurcien?) any idea what may need to be changed due to recent changes in pcadapt? The relevant code chunk is this:

K <- 25
x <- pcadapt(genotype_file, K = K)
plot(x, option = "screeplot") # 19 groups seems to be the correct value
plot(x, option = "scores", pop = sel[, 1]) # how populations are shared among the 19 groups

K <- 19
x <- pcadapt(genotype_file, K = K, min.maf = 0)

summary(x) # numerical quantities obtained after performing a PCA
plot(x, option = "manhattan")
plot(x, option = "qqplot", threshold = 0.1)
plot(x, option = "stat.distribution") # Distribution of Mahalanobis distances.

qval <- qvalue(x$pvalues)$qvalues
alpha <- 0.1
outliers_pcadapt <- which(qval < alpha)
print(outliers_pcadapt)
length(outliers_pcadapt) # 14 outliers

alpha <- 0.05 # use of a more stringent threshold to detect outliers
outliers <- which(qval < alpha)
print(outliers)
length(outliers) # 14 outliers

Maybe pcadapt() doesn't expect the file name but the content in some kind of structure?

[zkamvar]

Maybe pcadapt() doesn't expect the file name but the content in some kind of structure?

I just looked into it and this is the case. In 3.0.3, support for passing matrices to the pcadapt() function was removed so users had to supply an input file (as was noted in the preceeding chunk). According to the current documents on the pcadapt vignette, this functionality was restored. The only catch is that it expects the input to be a class of pcadapt_matrix, which, for all intents and purposes is a matrix with a different class assignment.

The new workflow for this now looks to be:

1. read in matrix with read.csv() and clean
2. convert to pcadapt matrix with read.pcadapt()
3. Run pcadapt() on the new pcadapt_matrix object.

Since Brenna's PR is a maintenance PR, I think it would be fine to make the changes directly on this branch instead of creating a new one. Would you be okay with that @hlapp and @BrennaF?

[zkamvar]

For posterity, the new PCAdapt chunks should look like this:

genotype <- sel[, 3:ncol(sel)]
dim(genotype)
# PCAdapt requires a pcadapt_class object. You can convert a matrix to 
# pcadapt_class with the read.pcadapt() function.
pca_genotype <- read.pcadapt(t(genotype))

K <- 25
x <- pcadapt(pca_genotype, K = K)
plot(x, option = "screeplot") # 19 groups seems to be the correct value
plot(x, option = "scores", pop = sel[, 1]) # how populations are shared among the 19 groups

K <- 19
x <- pcadapt(pca_genotype, K = K, min.maf = 0)

summary(x) # numerical quantities obtained after performing a PCA
plot(x, option = "manhattan")
plot(x, option = "qqplot", threshold = 0.1)
plot(x, option = "stat.distribution") # Distribution of Mahalanobis distances.

qval <- qvalue(x$pvalues)$qvalues
alpha <- 0.1
outliers_pcadapt <- which(qval < alpha)
print(outliers_pcadapt)
length(outliers_pcadapt) # 14 outliers

alpha <- 0.05 # use of a more stringent threshold to detect outliers
outliers <- which(qval < alpha)
print(outliers)
length(outliers) # 14 outliers

[hlapp]

Since Brenna's PR is a maintenance PR, I think it would be fine to make the changes directly on this branch instead of creating a new one. Would you be okay with that @hlapp and @BrennaF?

Yes I would. (Note BTW that you can edit this file in the PR directly, without having to clone @BrennaF's fork to get access to her patch-1 branch. I believe that only works on the Github UI, though, so you wouldn't be able to test locally this way.)

[zkamvar]

Note BTW that you can edit this file in the PR directly, without having to clone @BrennaF's fork to get access to her patch-1 branch. I believe that only works on the Github UI, though, so you wouldn't be able to test locally this way.

I've actually used this before 😁. It's quite nice and... you CAN clone it. When @BrennaF gives me the go-ahead, I'll do just that.

[BrennaF]

Fine with me all! Thank you! Brenna

Sent from my typo-enabled device.

On May 17, 2018, at 5:59 PM, Zhian N. Kamvar notifications@github.com wrote:

Note BTW that you can edit this file in the PR directly, without having to clone @BrennaF's fork to get access to her patch-1 branch. I believe that only works on the Github UI, though, so you wouldn't be able to test locally this way.

I've actually used this before 😁. It's quite nice and... you CAN clone it. When @BrennaF gives me the go-ahead, I'll do just that.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.

[zkamvar]

I have gone ahead and made the changes listed above.

@BrennaF, I noticed that the following line was giving a warning:

summary(eigenvals(wolf.rda, constrained=T))
## Argument `constrained` is deprecated; use `model` instead.

I changed it in https://github.com/NESCent/popgenInfo/pull/207/commits/bca6fd512a9c24fcc9a02ebf87ee0cf66c8ec5f5 and it no longer emits the warning. Is that okay with you?

[hlapp]

Thanks @BrennaF and @zkamvar !!

[BrennaF]

Sorry all; I’m traveling and without good internet service. This all sounds good - let me know if you need anything from me!

Brenna

Sent from my typo-enabled device.

On May 18, 2018, at 10:15 AM, Hilmar Lapp notifications@github.com wrote:

Thanks @BrennaF and @zkamvar !!

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or mute the thread.