FOR PCA ANALYSIS
vcp<-matrix(scan("g_VCP.txt",n=184287*20, sep=","), nrow=20,ncol=184287,byrow=T)
sla<-matrix(scan("g_SLA.txt",n=184287*18, sep=","), nrow=18,ncol=184287,byrow=T)
combined<-rbind(vcp,sla)
vars<-apply(combined, 2, var)
which(vars==0)
comb<-combined[,-which(vars==0)]
>dim(combined)
> dim(comb)
> pca1<-prcomp(comb, scale=TRUE)
> plot(pca1$x[,1],pca1$x[,2])
FOR LDA ANALYSIS
library(MASS)
> comb.lda<-lda(comb)
Error in lda.default(x, grouping, ...) :
argument "grouping" is missing, with no default
> comb.lda<-lda(rep(c("A","B"),19)comb)
Error: unexpected symbol in "comb.lda<-lda(rep(c("A","B"),19)comb"
> grp<-c(rep("SVC",20),rep("VCP",18)
>lda.out<-lda(grp ~ pca1$x[,1] + pca1$x[,2] + pca1$x[,3],CV=TRUE)
> summary(lda.out)
> lda.out$posterior
>barplot(t(pca1_lda$posterior),beside=FALSE)
To save the plot:
jpeg('pca2_lda_plot.jpg')
> barplot(t(pca2_lda$posterior),beside=FALSE)
> dev.off()
pca1: vcp(20 individuals) + sla(18 individuals) (File in R: comb)
pca2: wal(20 individuals) + gvl(18 individuals) (File in R: combined2)
pca3: wal(20 individuals) + vcp (20 individuals) (File in R: combined3)
pca4: svc (20 individuals) + vcp (20 individuals) (File in R: combined4)
pca5: svc (20 individuals) + sla (18 individuals) (File in R: combined5)
I saved the pca and pca_lda plots in dpca folder.
Summaries of pca and lda are saved in file summaries_pca_lda in folder dpca.
Means for pca_lda:
pca1_lda:
SLA VCP
0.4739034 0.5260966
pca2_lda
GVL WAL
0.5419368 0.4580632
colMeans(pca3_lda$posterior)
VCP WAL
0.5 0.5
colMeans(pca4_lda$posterior)
SVC VCP
0.5 0.5
colMeans(pca5_lda$posterior)
SLA SVC
0.4736842 0.5263158
WORK done on June 2, 2015
> vcp<-matrix(scan("g_VCP.txt",n=184287*20, sep=","), nrow=20,ncol=184287,byrow=T)
Read 3685740 items
> sla<-matrix(scan("g_SLA.txt",n=184287*18, sep=","), nrow=18,ncol=184287,byrow=T)
Read 3317166 items
> svc<-matrix(scan("g_SVC.txt",n=184287*20, sep=","), nrow=20,ncol=184287,byrow=T)
Read 3685740 items
> gvl<-matrix(scan("g_GVL.txt",n=184287*18, sep=","), nrow=18,ncol=184287,byrow=T)
Read 3317166 items
> wal<-matrix(scan("g_WAL.txt",n=184287*20, sep=","), nrow=20,ncol=184287,byrow=T)
Read 3685740 items
> vcp_sla<-rbind(vcp,sla)
> wal_gvl<-rbind(wal,gvl)
> wal_vcp<-rbind(wal,vcp)
> svc_vcp<-rbind(svc,vcp)
> svc_sla<-rbind(svc,sla)
> vars_vs<-apply(vcp_s
vcp_sla vcp_svc
> vars_vs<-apply(vcp_sla,2,var)
> pop_vs<-vcp_sla[,-which(vars_vs==0)]
> pca_vs<-prcomb(pop_vs,center=TRUE,scale=FALSE)
Error: could not find function "prcomb"
> pca_vs<-prcomp(pop_vs,center=TRUE,scale=FALSE)
> plot(pca_vs$x[,1],pca_vs$x[,2])
> vars_wg<-apply(wal_gvl,2,var)
> pop_wg<-wal_gvl[,-which(vars_wg==0)]
> pca_wg<-prcomp(pop_wg,center=TRUE,scale=FALSE)
Error in svd(x, nu = 0) : a dimension is zero
> pca_wg<-prcomp(pop_wg,scale=FALSE)
Error in svd(x, nu = 0) : a dimension is zero
> vars_wv<-apply(wal_vcp,2,var)
> pop_wv<-wal_vcp[,-which(vars_wg==0)]
> pca_wv<-prcomp(pop_wv,scale=FALSE)
Error in svd(x, nu = 0) : a dimension is zero
> cov_wg<-matrix(NA,nrow=38,ncol=38)
> for(i in 1:38){for(j in 1:38){
+ cov_wg[i,j]<-cov(wal_gvl[i,],wal_gvl[j,])
+ }}
> pca_wg<-prcomp(cov_wg,center=TRUE,scale=FALSE)
> plot(cov_wg$x[,1],cov_wg$x[,2])
Error in cov_wg$x : $ operator is invalid for atomic vectors
> plot(cov_wg$x[,1] + cov_wg$x[,2])
Error in cov_wg$x : $ operator is invalid for atomic vectors
> plot(cov_wg$x[,1] + cov_wg$x[,2],col=c(rep("orange",20),rep("blue",18))
+ )
Error in cov_wg$x : $ operator is invalid for atomic vectors
>
Standard deviation 1.314e-13
Proportion of Variance 0.000e+00
Cumulative Proportion 1.000e+00
> library(MASS)
> pca_wg<-prcomp(wal_gvl,center=TRUE,scale=FALSE)
> pca_wv<-prcomp(wal_vcp,center=TRUE,scale=FALSE)
> pca_sv<-prcomp(svc_vcp,center=TRUE,scale=FALSE)
> pca_ss<-prcomp(svc_sla,center=TRUE,scale=FALSE)
> plot(pca_wg$x[,1],pca_wg$x[,2])
> plot(pca_wv$x[,1],pca_wv$x[,2])
> plot(pca_sv$x[,1],pca_sv$x[,2])
> plot(pca_ss$x[,1],pca_ss$x[,2])
> barplot(t(pca2_lda$posterior),beside=FALSE)
> popid_wg<-c(rep("WAL",20),rep("GVL",18))
> lda_wg<-lda(popid_wg ~ pca_wg$[,1] + pca_w
pca_wg pca_wv
> lda_wg<-lda(popid_wg ~ pca_wg$[,1] + pca_wg$[,2] + pca_wg$[,3], CV=TRUE)
Error: unexpected '[' in "lda_wg<-lda(popid_wg ~ pca_wg$["
> ldout<-lda(popid_wg ~ pca_wg$x[,1] + pca_wg$x[,2] + pca_wg$x[,3],CV=TRUE)
> lda_wg<-lda(popid_wg ~ pca_wg$x[,1] + pca_wg$x[,2] + pca_wg$x[,3],CV=TRUE)
>
> popid_sv<-c(rep("SVC",20),rep("VCP",20))
> lda_sv<-lda(popid_sv ~ pca_sv$x[,1] + pca_sv$x[,2] + pca_sv$x[,3],CV=TRUE)
> popid_ss<-c(rep("SVC",20),rep("SLA",20))
> lda_ss<-lda(popid_ss ~ pca_ss$x[,1] + pca_ss$x[,2] + pca_ss$x[,3],CV=TRUE)
Error in model.frame.default(formula = popid_ss ~ pca_ss$x[, 1] + pca_ss$x[, :
variable lengths differ (found for 'pca_ss$x[, 1]')
> popid_ss<-c(rep("SVC",20),rep("SLA",18))
> lda_ss<-lda(popid_ss ~ pca_ss$x[,1] + pca_ss$x[,2] + pca_ss$x[,3],CV=TRUE)
> cov_wv<-matrix(NA,nrow=40,ncol=40)
> for(i in 1:40){for(j in 1:40){
+ cov_wv[i,j]<-cov(wal_vcp[i,],wal_vcp[j,])
+ }}
> pca_cov__wv<-prcomp(cov_wv,center=TRUE,scale=FALSE)
> pca_cov__wv<-prcomp(cov_wv,center=TRUE,scale=FALSE)
> pca2wv<-prcomp(cov_wv,center=TRUE,scale=FALSE)
> lda_ss<-lda(popid_ss ~ pca_ss$x[,1] + pca_ss$x[,2] + pca2$x[,3],CV=TRUE)
pca2 pca2_lda pca2wv
> lda2wv<-lda(popid_wv ~ pca2wv$x[,1] + pca2wv$x[,2] + pca2wv$x[,3],CV=TRUE)
> plot(pca2wv$x[,1],pca2wv$x[,2])
> cov_sv<-matrix(NA,nrow=40,ncol=40)
> for(i in 1:40){for(j in 1:40){
+ cov_sv[i,j]<-cov(svc_vcp[i,],svc_vcp[j,])
+ }}
> pca2sv<-prcomp(cov_sv,center=TRUE,scale=FALSE)
> lda2sv<-lda(popid_sv ~ pca2sv$x[,1] + pca2sv$x[,2] + pca2sv$x[,3],CV=TRUE)
> plot(pca2sv$x[,1],pca2sv$x[,2])
> cov_ss<-matrix(NA,nrow=38,ncol=38)
> for(i in 1:38){for(j in 1:38){
+ cov_ss[i,j]<-cov(svc_sla[i,],svc_sla[j,])
+ }}
> plot(pca2ss$x[,1],pca2ss$x[,2])
Error in plot(pca2ss$x[, 1], pca2ss$x[, 2]) : object 'pca2ss' not found
> pca2ss<-prcomp(cov_ss,center=TRUE,scale=FALSE)
> lda2ss<-lda(popid_ss ~ pca2ss$x[,1] + pca2ss$x[,2] + pca2ss$x[,3],CV=TRUE)
> plot(pca2ss$x[,1],pca2ss$x[,2])
> svc_wal<-rbind(svc,wal)
> svc_gvl<-rbind(svc,gvl)
> sla_gvl<-rbind(sla,gvl)
> vcp_gvl<-rbind(vcp,gvl)
> pc_sw<-prcomp(svc_wal,center=TRUE,sc
scale scan screen scale.=
scale.default scatter.smooth screeplot
> pc_sw<-prcomp(svc_wal,center=TRUE,scale=F
F Filter Find Formaldehyde FALSE
> pc_sw<-prcomp(svc_wal,center=TRUE,scale=FALSE)
> pca_sw<-prcomp(svc_wal,center=TRUE,scale=FALSE)
> plot(pca_sw$x[,1],pca_sw$x[,2])
> pca_sg<-prcomp(svc_gvl,center=TRUE,scale=FALSE)
> plot(pca_sg$x[,1],pca_sg$x[,2])
> pca_slg<-prcomp(sla_gvl,center=TRUE,scale=FALSE)
> plot(pca_slg$x[,1],pca_slg$x[,2])
> pca_vcg<-prcomp(vcp_gvl,center=TRUE,scale=FALSE)
> plot(pca_vcg$x[,1],pca_vcg$x[,2])