From 46c8911e29b564a33a6347f7dddc8771a86a72eb Mon Sep 17 00:00:00 2001 From: jiabowang Date: Wed, 17 Apr 2024 11:11:45 +0800 Subject: [PATCH] DeBUG taxa index and GTindex --- R/GAPIT.Bread.R | 4 +- R/GAPIT.EMMAxP3D.R | 2 +- R/GAPIT.IC.R | 12 +- R/GAPIT.Main.R | 164 +++++++++++---------- R/GAPIT.Power.compare.R | 1 + R/GAPIT.SS.R | 16 ++- R/GAPIT.SUPER.GS.R | 312 ++++++++++++++++++++-------------------- 7 files changed, 262 insertions(+), 249 deletions(-) diff --git a/R/GAPIT.Bread.R b/R/GAPIT.Bread.R index 905051b..6ae77b5 100644 --- a/R/GAPIT.Bread.R +++ b/R/GAPIT.Bread.R @@ -55,7 +55,7 @@ if(method=="MLM"){ GM=GM, model="MLM", QC=FALSE, - CV.Extragenetic=CV.Extragenetic, + CV.Extragenetic=CV.Extragenetic, # GTindex=GTindex, file.output=file.output ) @@ -110,7 +110,7 @@ if(method=="FaST" | method=="SUPER"| method=="DC") #GPS=myFaSTREML$GPS } -mySUPERFaST=GAPIT.SUPER.FastMLM(ys=matrix(Y[,-1],nrow(Y),1),X0=as.matrix(cbind(matrix(1,nrow(CV),1),CV[,-1])),snp.pool=as.matrix(GK[-1]), xs=as.matrix(GD[GTindex,-1]),vg=vg,delta=delta,LD=LD,method=method) +mySUPERFaST=GAPIT.SUPER.FastMLM(ys=matrix(Y[,-1],nrow(Y),1),X0=as.matrix(cbind(matrix(1,nrow(CV),1),CV[,-1])),snp.pool=as.matrix(GK[-1]), xs=as.matrix(GD[,-1]),vg=vg,delta=delta,LD=LD,method=method) GWAS=cbind(GM,mySUPERFaST$ps,mySUPERFaST$stats,mySUPERFaST$dfs,mySUPERFaST$effect) }#End of if(method=="FaST" | method=="SUPER") diff --git a/R/GAPIT.EMMAxP3D.R b/R/GAPIT.EMMAxP3D.R index 9eb38b1..03ed630 100644 --- a/R/GAPIT.EMMAxP3D.R +++ b/R/GAPIT.EMMAxP3D.R @@ -387,7 +387,7 @@ { xs=NULL }else{ - xs=myFRG$GD[GTindex,] + xs=myFRG$GD } if(!is.null(myFRG$GI)) diff --git a/R/GAPIT.IC.R b/R/GAPIT.IC.R index ff9ba6a..bad9acc 100644 --- a/R/GAPIT.IC.R +++ b/R/GAPIT.IC.R @@ -12,6 +12,7 @@ print("GAPIT.IC in process...") CV=DP$CV GD=DP$GD noCV=FALSE + # print(dim(CV)) if(is.null(CV)) { noCV=TRUE @@ -38,16 +39,21 @@ print("GAPIT.IC in process...") taxa_GD=as.character(GD[,1]) taxa_KI=as.character(DP$KI[,1]) taxa_CV=as.character(CV[,1]) + # print(Y) + # print(tail(taxa_Y)) + # print(tail(taxa_GD)) + # print(tail(taxa_CV)) + taxa_comall=intersect(intersect(intersect(taxa_KI,taxa_GD),taxa_Y),taxa_CV) taxa_g_cv=intersect(intersect(taxa_KI,taxa_GD),taxa_CV) # print(length(taxa_comall)) + # print(length(taxa_g_cv)) comCV=CV[taxa_CV%in%taxa_comall,] comCV <- comCV[match(taxa_comall,as.character(comCV[,1])),] comY=Y[taxa_Y%in%taxa_comall,] comY <- comY[match(taxa_comall,as.character(comY[,1])),] comGD=GD[taxa_GD%in%taxa_comall,] comGD <- comGD[match(taxa_comall,as.character(comGD[,1])),]# comGD=NULL - }else{ # print("@@@@") taxa_GD=as.character(GD[,1]) @@ -94,6 +100,10 @@ print("GAPIT.IC in process...") }#end of K }# end of GD # print(DP$KI[1:5,1:5]) + # print(tail(comY[,1])) + # print(tail(comGD[,1])) + # print(tail(comCV[,1])) + # print(tail(GD[,1])) GT=as.matrix(as.character(taxa_comall)) print(paste("There are ",length(GT)," common individuals in genotype , phenotype and CV files.",sep="")) if(nrow(comCV)!=length(GT))stop ("GAPIT says: The number of individuals in CV does not match to the number of individuals in genotype files.") diff --git a/R/GAPIT.Main.R b/R/GAPIT.Main.R index 453a0d9..dbc7fbf 100644 --- a/R/GAPIT.Main.R +++ b/R/GAPIT.Main.R @@ -104,6 +104,8 @@ function(Y, G=NULL, GD=NULL, + allGD=NULL, + allCV=NULL, GM=NULL, KI=NULL, Z=NULL, @@ -243,10 +245,12 @@ function(Y, #print(model) SUPER_GS_GAPIT = GAPIT.SUPER.GS(Y=Y, GD=GD, + allGD=allGD, GM=GM, KI=KI, Z=Z, CV=CV, + allCV=allCV, GK=GK, kinship.algorithm=kinship.algorithm, bin.from=bin.from, @@ -271,7 +275,7 @@ function(Y, sangwich.bottom=sangwich.bottom, QC=QC, QTN.gs=QTN.gs, - GTindex=GTindex, + # GTindex=GTindex, LD=LD, file.output=FALSE, GAPIT3.output=GAPIT3.output, @@ -347,7 +351,7 @@ function(Y, } if(is.null(GT)) { - GT=as.character(GD[,1]) + GT=as.character(CV[,1]) } #print("@@@@@@@@") #print(GD) @@ -355,19 +359,19 @@ function(Y, # if(PCA.total>0&!is.null(CV))CV=GAPIT.CVMergePC(CV,PC) # if(PCA.total>0&is.null(CV))CV=PC #for GS merge CV with GD name - if (is.null(CV)){ - my_allCV=CV - }else{ - taxa_GD=rownames(GD) - my_allCV=CV[order(CV[,1]),] - my_allCV=my_allCV[my_allCV[,1]%in%taxa_GD,] + # if (is.null(CV)){ + # my_allCV=CV + # }else{ + # taxa_GD=rownames(GD) + my_allCV=allCV + # my_allCV=my_allCV[my_allCV[,1]%in%taxa_GD,] #print(dim(my_allCV)) - } + # } #Handler of CV.Extragenetic - if(is.null(CV) & !is.null(CV.Extragenetic)){ - stop ("GAPIT says: CV.Extragenetic is more than avaiables.") - } + # if(is.null(CV) & !is.null(CV.Extragenetic)){ + # stop ("GAPIT says: CV.Extragenetic is more than avaiables.") + # } if(!is.null(CV)& !is.null(CV.Extragenetic)){ if(CV.Extragenetic>(ncol(CV)-1)){ @@ -377,8 +381,8 @@ function(Y, #Create Z as identity matrix from Y if it is not provided if(kinship.algorithm!="None" & kinship.algorithm!="SUPER" & is.null(Z)){ - taxa=as.character(Y[,1]) #this part will make GS without CV not present all prediction - Z=as.data.frame(diag(1,nrow(Y))) + taxa=as.character(CV[,1]) #this part will make GS without CV not present all prediction + Z=as.data.frame(diag(1,nrow(CV))) #taxa=as.character(KI[,1]) #Z=as.data.frame(diag(1,nrow(KI))) Z=rbind(taxa,Z) @@ -394,16 +398,16 @@ function(Y, } #Create CV with all 1's if it is not provided - noCV=FALSE - if(is.null(CV)){ - noCV=TRUE - CV=Y[,1:2] - CV[,2]=1 - colnames(CV)=c("taxa","overall") - } + # noCV=FALSE + # if(is.null(CV)){ + # noCV=TRUE + # CV=Y[,1:2] + # CV[,2]=1 + # colnames(CV)=c("taxa","overall") + # } #Remove duplicat and integragation of data - print("QC is in process...") + # print("QC is in process...") CVI <- CV @@ -422,13 +426,13 @@ function(Y, # } # print(length(GK)) - GTindex=match(as.character(KI[,1]),as.character(Y[,1])) - GTindex=GTindex[!is.na(GTindex)] + # GTindex=match(as.character(KI[,1]),as.character(Y[,1])) + # GTindex=GTindex[!is.na(GTindex)] # KI=KI[GTindex,GTindex+1] - my_taxa=as.character(KI[,1]) - CV=CV[as.character(CV[,1])%in%as.character(Y[,1]),] + # my_taxa=as.character(KI[,1]) + # CV=CV[as.character(CV[,1])%in%as.character(Y[,1]),] #Output phenotype - colnames(Y)=c("Taxa",name.of.trait) + # colnames(Y)=c("Taxa",name.of.trait) # if(file.output){ # try(utils::write.table(Y, paste("GAPIT.", name.of.trait,".phenotype.csv" ,sep = ""), # quote = FALSE, sep = ",", @@ -439,45 +443,45 @@ function(Y, # Default kinship.algorithm = "VanRaden". # This if() may be seldom used. #TDP - if( kinship.algorithm =="None" ){ - if(min(CV[,2])==max(CV[,2])) CV=NULL +# if( kinship.algorithm =="None" ){ +# if(min(CV[,2])==max(CV[,2])) CV=NULL - # GAPIT.TDP() does not appear to exist. - # theTDP = GAPIT.TDP(Y=Y, - # CV=CV, - # SNP = as.data.frame(cbind(GT[GTindex],as.matrix(as.data.frame(GD[GTindex,])))), - # QTN=QTN, - # Round=QTN.round, - # QTN.limit=QTN.limit, - # QTN.update=QTN.update, - # Method=QTN.method - # ) -#print(dim(GM)) -#print(length(theTDP$p)) - - #theGWAS=cbind(GM,theTDP$p,NA,NA,NA) - theGWAS=cbind(GM,NA,NA,NA,NA) +# # GAPIT.TDP() does not appear to exist. +# # theTDP = GAPIT.TDP(Y=Y, +# # CV=CV, +# # SNP = as.data.frame(cbind(GT[GTindex],as.matrix(as.data.frame(GD[GTindex,])))), +# # QTN=QTN, +# # Round=QTN.round, +# # QTN.limit=QTN.limit, +# # QTN.update=QTN.update, +# # Method=QTN.method +# # ) +# #print(dim(GM)) +# #print(length(theTDP$p)) + +# #theGWAS=cbind(GM,theTDP$p,NA,NA,NA) +# theGWAS=cbind(GM,NA,NA,NA,NA) - return (list(Compression = NULL, - kinship.optimum = NULL, - kinship = NULL, - PC = NULL, - GWAS = theGWAS, - GPS = NULL, - Pred = NULL, - REMLs = NULL, - #QTN = theTDP$QTN, - QTN = NULL, - Timmer = Timmer, - Memory = Memory, - h2 = NULL)) - } +# return (list(Compression = NULL, +# kinship.optimum = NULL, +# kinship = NULL, +# PC = NULL, +# GWAS = theGWAS, +# GPS = NULL, +# Pred = NULL, +# REMLs = NULL, +# #QTN = theTDP$QTN, +# QTN = NULL, +# Timmer = Timmer, +# Memory = Memory, +# h2 = NULL)) +# } #rm(qc) - gc() + # gc() - Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="QC") - Memory=GAPIT.Memory(Memory=Memory,Infor="QC") + # Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="QC") + # Memory=GAPIT.Memory(Memory=Memory,Infor="QC") #Get indicator of sangwich top and bottom byPass.top=FALSE @@ -509,29 +513,32 @@ function(Y, #print(GD[1:5,1:5]) #Create GK if not provided - if(is.null(GK)){ + if(is.null(GK)) + { # set.seed(1) nY=floor(nrow(Y)*.9) nG=ncol(GD) - if(nG>nY){ + if(nG>nY) + { snpsam=sample(1:nG,nY) }else{ snpsam=1:nG } - GK=GD[GTindex,snpsam] + GK=GD # print(dim(GK)) # print(GK[270:279,1:5]) SNPVar=apply(as.matrix(GK), 2, stats::var) # print(SNPVar) GK=GK[,SNPVar>0] - GK=cbind(as.data.frame(GT[GTindex]),as.data.frame(GK)) #add taxa + GK=cbind(as.data.frame(GT),as.data.frame(GK)) #add taxa } #myGD=cbind(as.data.frame(GT),as.data.frame(GD)) - # file.output.temp=file.output - # file.output=FALSE - #print(sangwich.top)[GTindex,c(1,GTindex+1)] - GP=GAPIT.Bread(Y=Y,CV=CV,Z=Z,KI=KI,GK=GK,GD=cbind(as.data.frame(GT[GTindex]),as.data.frame(GD[GTindex,])),GM=GI,method=sangwich.top,GTindex=GTindex,LD=LD,file.output=FALSE)$GWAS + print(tail(Y[,1])) + print(tail(CV[,1])) + print(tail(GT)) + # print(dim(cbind(as.data.frame(GT[GTindex]),as.data.frame(GD[GTindex,])))) + GP=GAPIT.Bread(Y=Y,CV=CV,Z=Z,KI=KI,GK=GK,GD=cbind(as.data.frame(GT),as.data.frame(GD)),GM=GI,method=sangwich.top,LD=LD,file.output=FALSE)$GWAS # file.output=file.output.temp @@ -695,7 +702,7 @@ function(Y, # print(optOnly) # print(head(CVI)) p3d <- GAPIT.EMMAxP3D(ys=ys, - xs=as.matrix(as.data.frame(GD[GTindex,colInclude])), + xs=as.matrix(as.data.frame(GD[,colInclude])), K = as.matrix(bk$KW), Z=matrix(as.numeric(as.matrix(zc$Z[,-1])),nrow=zrow,ncol=zcol), X0=X0, @@ -721,7 +728,7 @@ function(Y, file.GM=file.GM, file.Ext.GD=file.Ext.GD, file.Ext.GM=file.Ext.GM, - GTindex=GTindex, + # GTindex=GTindex, genoFormat=genoFormat, optOnly=optOnly, SNP.effect=SNP.effect, @@ -828,16 +835,15 @@ function(Y, Memory=GAPIT.Memory(Memory=Memory,Infor="Genotype for burger") print(paste("bin---",bin,"---inc---",inc,sep="")) - GK=GD[GTindex,myGenotype$SNP.QTN] + GK=GD[,myGenotype$SNP.QTN] SUPER_GD=GD[,myGenotype$SNP.QTN] SNPVar=apply(as.matrix(GK), 2, stats::var) - GK=GK[,SNPVar>0] SUPER_GD=SUPER_GD[,SNPVar>0] - GK=cbind(as.data.frame(GT[GTindex]),as.data.frame(GK)) #add taxa + GK=cbind(as.data.frame(GT),as.data.frame(GK)) #add taxa # print(length(GT)) # print(dim(SUPER_GD)) - SUPER_GD=cbind(as.data.frame(GT[GTindex]),as.data.frame(SUPER_GD)) #add taxa + SUPER_GD=cbind(as.data.frame(GT),as.data.frame(SUPER_GD)) #add taxa # print(dim(GK)) #GP=NULL }# end of if(is.null(GK)) @@ -946,7 +952,7 @@ if(Model.selection == TRUE){ p3d <- GAPIT.EMMAxP3D(ys=ys,xs=as.matrix(as.data.frame(GD[,1])),K = as.matrix(bk$KW) ,Z=Z1,X0=X0.test,CVI=CVI,CV.Extragenetic=CV.Extragenetic,GI=GI,SNP.P3D=SNP.P3D,Timmer=Timmer,Memory=Memory,fullGD=fullGD, SNP.permutation=SNP.permutation, GP=GP, file.path=file.path,file.from=file.from,file.to=file.to,file.total=file.total, file.fragment = file.fragment, byFile=byFile, file.G=file.G,file.Ext.G=file.Ext.G,file.GD=file.GD, file.GM=file.GM, file.Ext.GD=file.Ext.GD,file.Ext.GM=file.Ext.GM, - GTindex=GTindex,genoFormat=genoFormat,optOnly=TRUE,SNP.effect=SNP.effect,SNP.impute=SNP.impute,name.of.trait=name.of.trait, Create.indicator = Create.indicator, Major.allele.zero = Major.allele.zero) + genoFormat=genoFormat,optOnly=TRUE,SNP.effect=SNP.effect,SNP.impute=SNP.impute,name.of.trait=name.of.trait, Create.indicator = Create.indicator, Major.allele.zero = Major.allele.zero) @@ -1122,10 +1128,10 @@ print("-------------- Sandwich bottom with raw burger------------------------") print("--------------EMMAxP3D with the optimum setting-----------------------") #print(dim(ys)) #print(dim(as.matrix(as.data.frame(GD[GTindex,colInclude])))) - p3d <- GAPIT.EMMAxP3D(ys=ys,xs=as.matrix(as.data.frame(GD[GTindex,colInclude])) ,K = as.matrix(bk$KW) ,Z=Z1,X0=as.matrix(X0),CVI=CVI, CV.Extragenetic=CV.Extragenetic,GI=GI,SNP.P3D=SNP.P3D,Timmer=Timmer,Memory=Memory,fullGD=fullGD, + p3d <- GAPIT.EMMAxP3D(ys=ys,xs=as.matrix(as.data.frame(GD[,colInclude])) ,K = as.matrix(bk$KW) ,Z=Z1,X0=as.matrix(X0),CVI=CVI, CV.Extragenetic=CV.Extragenetic,GI=GI,SNP.P3D=SNP.P3D,Timmer=Timmer,Memory=Memory,fullGD=fullGD, SNP.permutation=SNP.permutation, GP=GP, file.path=file.path,file.from=file.from,file.to=file.to,file.total=file.total, file.fragment = file.fragment, byFile=byFile, file.G=file.G,file.Ext.G=file.Ext.G,file.GD=file.GD, file.GM=file.GM, file.Ext.GD=file.Ext.GD,file.Ext.GM=file.Ext.GM, - GTindex=GTindex,genoFormat=genoFormat,optOnly=optOnly,SNP.effect=SNP.effect,SNP.impute=SNP.impute,name.of.trait=name.of.trait, Create.indicator = Create.indicator, Major.allele.zero = Major.allele.zero) + genoFormat=genoFormat,optOnly=optOnly,SNP.effect=SNP.effect,SNP.impute=SNP.impute,name.of.trait=name.of.trait, Create.indicator = Create.indicator, Major.allele.zero = Major.allele.zero) Timmer=GAPIT.Timmer(Timmer=Timmer,Infor="GWAS") Memory=GAPIT.Memory(Memory=Memory,Infor="GWAS") @@ -1144,7 +1150,7 @@ print("---------------Sandwich bottom: reload bins ---------------------------") #SUPER: Final screening GK=GK.save # print(GK) - myBread=GAPIT.Bread(Y=Y,CV=CV,Z=Z,GK=GK,GD=cbind(as.data.frame(GT[GTindex]),as.data.frame(GD)),GM=GI,method=sangwich.bottom,GTindex=GTindex,LD=LD,file.output=FALSE) + myBread=GAPIT.Bread(Y=Y,CV=CV,Z=Z,GK=GK,GD=cbind(as.data.frame(GT),as.data.frame(GD)),GM=GI,method=sangwich.bottom,LD=LD,file.output=FALSE) print("SUPER saving results...") diff --git a/R/GAPIT.Power.compare.R b/R/GAPIT.Power.compare.R index fe05968..879b431 100644 --- a/R/GAPIT.Power.compare.R +++ b/R/GAPIT.Power.compare.R @@ -40,6 +40,7 @@ colnames(Y)=c("Taxa","Simu") GD=GD, GM=GM, PCA.total=PCA.total, + cutOff=0.01, CV=CV, file.output=FALSE, model=all.method[j], diff --git a/R/GAPIT.SS.R b/R/GAPIT.SS.R index 87db69e..a02787c 100644 --- a/R/GAPIT.SS.R +++ b/R/GAPIT.SS.R @@ -129,10 +129,10 @@ if(DP$SNP.test) #CV.Extragenetic specified QTN.gs=ncol(GD2) CV.Extragenetic=DP$CV.Extragenetic - XCVI=XCV[,c((2+CV.Extragenetic):(ncol(XCV)-QTN.gs))] + if(CV.Extragenetic!=0)XCVI=XCV[,c((2+CV.Extragenetic):(ncol(XCV)-QTN.gs))] XCVN=XCV[,c(1:(1+CV.Extragenetic))] if(QTN.gs!=0)XCVqtn=XCV[,c((ncol(XCV)-QTN.gs):ncol(XCV))] - beta.I=lm.coeff[c((2+CV.Extragenetic):(ncol(XCV)-QTN.gs))] + if(CV.Extragenetic!=0)beta.I=lm.coeff[c((2+CV.Extragenetic):(ncol(XCV)-QTN.gs))] beta.N=lm.coeff[c(1:(1+CV.Extragenetic))] if(QTN.gs!=0)beta.QTN=lm.coeff[c((ncol(XCV)-QTN.gs):ncol(XCV))] BLUE.N=XCVN%*%beta.N @@ -257,9 +257,11 @@ if(DP$SNP.test) if(DP$PCA.total==0) ic_PCA=NULL gapitMain <- GAPIT.Main(Y=ic_Y, GD=IC$GD[,-1], + allGD=IC$myallGD[,-1], + allCV=IC$myallCV, GM=DP$GM, KI=ic_KI, - CV=IC$myallCV, + CV=IC$PCA, CV.Extragenetic=DP$CV.Extragenetic, GP=DP$GP, GK=DP$GK, @@ -295,7 +297,7 @@ if(DP$SNP.test) SNP.impute=DP$SNP.impute, PCA.total=DP$PCA.total, #GAPIT.Version=GAPIT.Version, - GT=DP$GT, + GT=IC$GT, SNP.fraction = DP$SNP.fraction, seed =DP$seed, BINS = DP$BINS, @@ -360,11 +362,13 @@ if(!is.null(GWAS))myPower=GAPIT.Power(WS=DP$WS, alpha=DP$alpha, maxOut=DP$maxOut # Here is Genomic Prediction function print("GAPIT will be into GS approach...") gapitMain <- GAPIT.Main(Y=IC$Y, - GD=DP$GD[,-1], + GD=IC$GD[,-1], + allGD=IC$allGD[,-1], GM=DP$GM, KI=IC$KI, Z=DP$Z, - CV=IC$myallCV, + CV=IC$PCA, + allCV=IC$myallCV, CV.Extragenetic=DP$CV.Extragenetic, GP=DP$GP, GK=DP$GK, diff --git a/R/GAPIT.SUPER.GS.R b/R/GAPIT.SUPER.GS.R index 3328fe8..967a4dc 100644 --- a/R/GAPIT.SUPER.GS.R +++ b/R/GAPIT.SUPER.GS.R @@ -46,10 +46,12 @@ `GAPIT.SUPER.GS`<- function(Y, GD = NULL, + allGD=NULL, GM = NULL, KI = NULL, Z = NULL, CV = NULL, + allCV=NULL, GK = NULL, kinship.algorithm = NULL, bin.from = 10000, @@ -110,7 +112,7 @@ print(paste("Y is empty. No GWAS/GS performed for ",name.of.trait,sep="")) return (list(compression=NULL,kinship.optimum=NULL, kinship=KI,PC=PC,GWAS=NULL, GPS=NULL,Pred=NULL, REMLs=NULL,Timmer=Timmer,Memory=Memory)) } print("------------Examining data (QC)------------------------------------------") -if(is.null(Y)) stop ("GAPIT says: Phenotypes must exist.") +# if(is.null(Y)) stop ("GAPIT says: Phenotypes must exist.") if(is.null(KI)&missing(GD) & kinship.algorithm!="SUPER") stop ("GAPIT says: Kinship is required. As genotype is not provided, kinship can not be created.") if(is.null(GD) & is.null(GT)) { GT=as.matrix(Y[,1]) @@ -120,8 +122,8 @@ if(is.null(GD) & is.null(GT)) { colnames(GI)=c("SNP","Chromosome","Position") } # print(cbind(CV,PC)) -if(PCA.total>0&!is.null(CV))CV=GAPIT.CVMergePC(CV,PC) -if(PCA.total>0&is.null(CV))CV=PC +# if(PCA.total>0&!is.null(CV))CV=GAPIT.CVMergePC(CV,PC) +# if(PCA.total>0&is.null(CV))CV=PC if(kinship.algorithm!="None" & kinship.algorithm!="SUPER" & is.null(Z)){ taxa=as.character(Y[,1]) @@ -134,52 +136,52 @@ if(kinship.algorithm!="None" & kinship.algorithm!="SUPER" & !is.null(Z)) { if(nrow(Z)-1nY){snpsam=sample(1:nG,nY)}else{snpsam=1:nG} - GK=GD[GTindex,snpsam] + GK=GD[,snpsam] SNPVar=apply(as.matrix(GK), 2, stats::var) #print(snpsam) # if(snpsam==1)stop ("GAPIT says: SUPER_GS must putin GD and GM.") GK=GK[,SNPVar>0] - GK=cbind(as.data.frame(GT[GTindex]),as.data.frame(GK)) #add taxa - + GK=cbind(as.data.frame(GT),as.data.frame(GK)) #add taxa } #print(head(CV)) #myGD=cbind(as.data.frame(GT),as.data.frame(GD)) @@ -211,116 +212,108 @@ print("-------------------start SUPER BREAD-----------------------------------") # file.output.temp=file.output # file.output=FALSE # print(memory.size()) - GP=GAPIT.Bread(Y=Y,CV=CV,Z=Z,KI=KI,GK=GK,GD=cbind(as.data.frame(GT),as.data.frame(GD)),GM=GI,method=sangwich.top,GTindex=GTindex,LD=LD,file.output=FALSE,CV.Extragenetic=CV.Extragenetic)$GWAS + GP=GAPIT.Bread(Y=Y,CV=CV,Z=Z,KI=KI,GK=GK,GD=cbind(as.data.frame(GT),as.data.frame(GD)),GM=GI,method=sangwich.top,LD=LD,file.output=FALSE,CV.Extragenetic=CV.Extragenetic)$GWAS # file.output=file.output.temp # print(memory.size()) GK=NULL -if(inclosure.to>nrow(Y)) ##########removed by Jiabo Wang ,unlimited number of inclosures -{ -inclosure.to=nrow(Y)-1 -print("the number of choosed inclosure is more than number of individuals") -print("Set the number of choosed incolosure max equal to individuals") -} -if(inclosure.from>inclosure.to) ##########removed by Jiabo Wang ,unlimited number of inclosures -{ -inclosure.from=inclosure.to -} -bin.level=seq(bin.from,bin.to,by=bin.by) -inclosure=seq(inclosure.from,inclosure.to,by=inclosure.by) + if(inclosure.to>nrow(Y)) ##########removed by Jiabo Wang ,unlimited number of inclosures + { + inclosure.to=nrow(Y)-1 + print("the number of choosed inclosure is more than number of individuals") + print("Set the number of choosed incolosure max equal to individuals") + } + if(inclosure.from>inclosure.to) ##########removed by Jiabo Wang ,unlimited number of inclosures + { + inclosure.from=inclosure.to + } + bin.level=seq(bin.from,bin.to,by=bin.by) + inclosure=seq(inclosure.from,inclosure.to,by=inclosure.by) #print(inclosure) -e=1 #################################number of bins and inclosure -count=0 -num_selection=length(bin.level)*length(inclosure) -SUPER_selection=matrix(,num_selection,6) -colnames(SUPER_selection)=c("bin","pseudo_QTNs","Max_pQTNs","REML","VA","VE") + e=1 #################################number of bins and inclosure + count=0 + num_selection=length(bin.level)*length(inclosure) + SUPER_selection=matrix(,num_selection,6) + colnames(SUPER_selection)=c("bin","pseudo_QTNs","Max_pQTNs","REML","VA","VE") #for (bin in bin.level){bin=bin.level[e]} #for (inc in inclosure){inc=inclosure[e]} -for (bin in bin.level){ -for (inc in inclosure){ -count=count+1 - mySpecify=GAPIT.Specify(GI=GI,GP=GP,bin.size=bin,inclosure.size=inc) - SNP.QTN=mySpecify$index - num_pseudo_QTN=length(mySpecify$CB) - num_bins=mySpecify$num_bins -#print(paste("bin---",bin,"---inc---",inc,sep="")) - GK=GD[GTindex,SNP.QTN] - SUPER_GD=GD[,SNP.QTN] - SNPVar=apply(as.matrix(GK), 2, stats::var) - GK=GK[,SNPVar>0] - SUPER_GD=SUPER_GD[,SNPVar>0] - GK=cbind(as.data.frame(GT[GTindex]),as.data.frame(GK)) #add taxa - SUPER_GD=cbind(as.data.frame(GT),as.data.frame(SUPER_GD)) #add taxa - myBurger=GAPIT.Burger(Y=Y,CV=CV,GK=GK) #modifed by Jiabo Wang - myREML=myBurger$REMLs - myVG=myBurger$vg - myVE=myBurger$ve -SUPER_selection[count,1]=bin -SUPER_selection[count,2]=num_pseudo_QTN -SUPER_selection[count,3]=num_bins -SUPER_selection[count,4]=myREML -SUPER_selection[count,5]=myVG -SUPER_selection[count,6]=myVE + for (bin in bin.level) + { + for (inc in inclosure) + { + count=count+1 + mySpecify=GAPIT.Specify(GI=GI,GP=GP,bin.size=bin,inclosure.size=inc) + SNP.QTN=mySpecify$index + num_pseudo_QTN=length(mySpecify$CB) + num_bins=mySpecify$num_bins +#print(paste("bin---",bin,"---inc---",inc,sep="")) + GK=GD[,SNP.QTN] + SUPER_GD=GD[,SNP.QTN] + SNPVar=apply(as.matrix(GK), 2, stats::var) + GK=GK[,SNPVar>0] + SUPER_GD=SUPER_GD[,SNPVar>0] + GK=cbind(as.data.frame(GT),as.data.frame(GK)) #add taxa + SUPER_GD=cbind(as.data.frame(GT),as.data.frame(SUPER_GD)) #add taxa + myBurger=GAPIT.Burger(Y=Y,CV=CV,GK=GK) #modifed by Jiabo Wang + myREML=myBurger$REMLs + myVG=myBurger$vg + myVE=myBurger$ve + SUPER_selection[count,1]=bin + SUPER_selection[count,2]=num_pseudo_QTN + SUPER_selection[count,3]=num_bins + SUPER_selection[count,4]=myREML + SUPER_selection[count,5]=myVG + SUPER_selection[count,6]=myVE #print(SUPER_selection[count,]) - if(count==1){ - GK.save=GK - LL.save=myREML - SUPER_optimum_GD=SUPER_GD ########### get SUPER GD -}else{ - if(myREMLnk){ group.from=nk #warning("The lower bound of groups is too high. It was set to the size of kinship!") - message("The lower bound of groups is too high. It was set to the size of kinship!") + print("The lower bound of groups is too high. It was set to the size of kinship!") } } @@ -371,7 +364,7 @@ if(!is.null(CV)){ group.from=ncol(CV)+2 group.to=ncol(CV)+2 #warning("The upper bound of groups (group.to) is not sufficient. both boundries were set to their minimum and GLM is performed!") - message("The upper bound of groups (group.to) is not sufficient. both boundries were set to their minimum and GLM is performed!") + print("The upper bound of groups (group.to) is not sufficient. both boundries were set to their minimum and GLM is performed!") } } @@ -380,7 +373,9 @@ if(missing("kinship.cluster")) kinship.cluster=c("ward", "single", "complete", " if(missing("kinship.group")) kinship.group=c("Mean", "Max", "Min", "Median") numSetting=length(GROUP)*length(kinship.cluster)*length(kinship.group) ys=as.matrix(Y[2]) -X0=as.matrix(CV[,-1]) +# print(dim(CV)) +# print(dim(allCV)) +X0=as.matrix(CV[,-1,drop=FALSE]) if(min(X0[,1])!=max(X0[,1])) X0 <- cbind(1, X0) #do not add overall mean if X0 has it already at first column hold_Z=Z @@ -400,7 +395,7 @@ for (group in GROUP) #if(!optOnly) {print("Compressing and Genome screening..." )} order_count=order_count+1 if(order_count==1)print("-------Mixed model with Kinship-----------------------------") -if(group