knitr::opts_chunk$set(fig.width=7, fig.asp=0.618)
library("breastCancerVDX")
library(Biobase)
library(gpowerr)
data(vdx)
patientdata <- pData(vdx)
expressions <- t(exprs(vdx))

PCA results

pca_results <- prcomp(expressions)
pca_pev <- pca_results$sdev[1]^2/sum(pca_results$sdev^2)

Index of sparseness

results <- rep(0, length(1:50))

for (i in 1:50) {
  pow <- gpower(data = expressions, k=1, rho = i/50, reg = "l1", center = TRUE)
  results[i] <-  pow$exp_var * pca_pev * pow$prop_sparse
}

plot

png("figures/ISplot.png", width = 600, height = 370, res=80)
plot(1:50/50, results, xlab="rho", ylab="Index of Sparseness")
dev.off()

max(results)
which.max(results) / 50

final weights

power <- gpower(data = expressions, k = 1, rho = 0.16)
scores <- power$scores
scores <- scores[!is.na(patientdata$grade)]
graded <- patientdata$grade[!is.na(patientdata$grade)]
cor(scores, graded)
cor.test(scores, graded)
genes <- row.names(power$weights)[power$weights != 0]
genes

write.table(genes, "genes.txt", quote = FALSE, col.names = FALSE, row.names = FALSE)


plofknaapje/gpowerpca documentation built on July 27, 2021, 4:17 a.m.