14、PCA分析

　　做芯片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))

 

 

可看到两段代码的结果是一样的

 

转载于:https://www.cnblogs.com/xianwen/p/4522383.html