#-------------------------------------------------------------------------------
# Written by: Kory R Johnson
# Email: johnsonko@mail.nih.gov
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
# Covariance-based PCA using Final Selected Head Age Markers
#-------------------------------------------------------------------------------

par(mfrow=c(1,2))
palette(rainbow(20))
Final.working.data <- read.table("Final_Normalized_Expression_For_Head_Age_Markers.txt",header=T,row.names=1)
col.2.use <- c(5,5,5,5,4,4,4,4,4,11,11,11,11,11,14,14,14,14,14,3,3,3,3,3,2,2,2,2,2)
pca.results <- princomp(Final.working.data,cor=F)
pca.loadings <- loadings(pca.results)
totalvar <- (pca.results$sdev^2)
varianceAll <- round(totalvar/sum(totalvar)*100,1)
variancePer <- varianceAll[1:3]
x.range <- range(pca.loadings[,1])
y.range <- range(pca.loadings[,2])
xlab <- paste(c("Principal Component 1 (",variancePer[1],"%)"),collapse="")	
ylab <- paste(c("Principal Component 2 (",variancePer[2],"%)"),collapse="")
plot(pca.loadings[,1],pca.loadings[,2],type='n',xlab=xlab,ylab=ylab)
points(pca.loadings[,1],pca.loadings[,2],lwd=2,bg=col.2.use,cex=3,pch=21)
x.range <- range(pca.loadings[,2])
y.range <- range(pca.loadings[,3])
xlab <- paste(c("Principal Component 2 (",variancePer[2],"%)"),collapse="")	
ylab <- paste(c("Principal Component 3 (",variancePer[3],"%)"),collapse="")
plot(pca.loadings[,2],pca.loadings[,3],type='n',xlab=xlab,ylab=ylab)
points(pca.loadings[,2],pca.loadings[,3],lwd=2,bg=col.2.use,cex=3,pch=21)
write.table(pca.results$scores,"PCA_Marker_Scores_For_Head_Age_Markers.txt",sep="\t")
write.table(varianceAll,"PCA_Variance_Per_Component_Using_Head_Age_Markers.txt",sep="\t")