/usr/lib/R/site-library/snpStats/doc/pca-vignette.R is in r-bioc-snpstats 1.24.0+dfsg-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 | ### R code from vignette source 'pca-vignette.Rnw'
###################################################
### code chunk number 1: get-data
###################################################
require(snpStats)
data(for.exercise)
controls <- rownames(subject.support)[subject.support$cc==0]
use <- seq(1, ncol(snps.10), 10)
ctl.10 <- snps.10[controls,use]
###################################################
### code chunk number 2: xxt-matrix
###################################################
xxmat <- xxt(ctl.10, correct.for.missing=FALSE)
###################################################
### code chunk number 3: eigen
###################################################
evv <- eigen(xxmat, symmetric=TRUE)
pcs <- evv$vectors[,1:5]
evals <- evv$values[1:5]
evals
###################################################
### code chunk number 4: pc-one
###################################################
pop <- subject.support[controls,"stratum"]
par(mfrow=c(1,2))
boxplot(pcs[,1]~pop)
boxplot(pcs[,2]~pop)
###################################################
### code chunk number 5: pre-multiply
###################################################
btr <- snp.pre.multiply(ctl.10, diag(1/sqrt(evals)) %*% t(pcs))
###################################################
### code chunk number 6: post-multiply
###################################################
pcs <- snp.post.multiply(snps.10[,use], t(btr))
###################################################
### code chunk number 7: testing
###################################################
cc <- subject.support$cc
uncorrected <- single.snp.tests(cc, snp.data=snps.10)
corrected <- snp.rhs.tests(cc~pcs[,1], snp.data=snps.10)
par(mfrow=c(1,2),cex.sub=0.85)
qq.chisq(chi.squared(uncorrected,1), df=1)
qq.chisq(chi.squared(corrected), df=1)
|