做芯片PCA主成分分析可以选择使用affycoretools包的plotPCA方法,以样品"GSM363445_LNTT.CEL"、"GSM362948_LTT.CEL"、"GSM363447_LNTT.CEL"、"GSM362949_LTT.CEL"、"GSM363449_LNTT.CEL"、"GSM362947_LTT.CEL"为例:
library(affy)
library(affycoretools)
rawData<-ReadAffy("GSM363445_LNTT.CEL","GSM362948_LTT.CEL","GSM363447_LNTT.CEL", "GSM362949_LTT.CEL","GSM363449_LNTT.CEL","GSM362947_LTT.CEL")
plotPCA(exprs(rawData))
不过affycoretools包比较大,也可以自行编程而不通过affycoretools包:
library(affy)
rawData<-ReadAffy("GSM363445_LNTT.CEL","GSM362948_LTT.CEL","GSM363447_LNTT.CEL","GSM362949_LTT.CEL","GSM363449_LNTT.CEL","GSM362947_LTT.CEL")
x <- t(exprs(rawData)[apply(exprs(rawData),1,function(r) {sum(is.na(r))==0}),])
x <- as.matrix(x)
x <- scale(x, center = T, scale = FALSE) ## arrayanalysis的PCA会默认scale = TRUE
s <- svd(x, nu = 0) ## svd奇异值分解,不需要计算u
pca1 <- list(sdev = s$d, rotation = s$v)
pca1$x <- x %*% s$v
plot(pca1$x[,1],pca1$x[,2],col=c(1:6))
legend("topleft",legend=sampleNames(rawData), lwd=1, col=c(1:6))
可看到两段代码的结果是一样的