From 12bfc1f61ae5cc7ee8b8b6867275d4672056afbe Mon Sep 17 00:00:00 2001 From: jiabowang Date: Thu, 22 Sep 2022 08:17:24 +0800 Subject: [PATCH] add matching taxa of KI in IC --- R/GAPIT.Blink.R | 1 + R/GAPIT.Genotype.R | 4 +- R/GAPIT.Genotype.View.R | 20 ++++++++-- R/GAPIT.IC.R | 84 +++++++++++++++++++++++------------------ R/GAPIT.Main.R | 8 ++-- R/GAPIT.RandomModel.R | 11 +++--- 6 files changed, 78 insertions(+), 50 deletions(-) diff --git a/R/GAPIT.Blink.R b/R/GAPIT.Blink.R index ba389f2..f15663b 100644 --- a/R/GAPIT.Blink.R +++ b/R/GAPIT.Blink.R @@ -335,6 +335,7 @@ P[P==0] <- min(P[P!=0],na.rm=TRUE)*0.01 P[is.na(P)] =1 # print(str(myGLM)) + gc() nf=ncol(myGLM$P)/4 tvalue=myGLM$P[,nf*2-shift] diff --git a/R/GAPIT.Genotype.R b/R/GAPIT.Genotype.R index 4846a3c..9b9b115 100644 --- a/R/GAPIT.Genotype.R +++ b/R/GAPIT.Genotype.R @@ -115,8 +115,8 @@ if(!needKinPC &SNP.fraction<1) stop("GAPIT says: You did not require calculate if(!SNP.test & is.null(KI) & !byData & !byFile) stop("GAPIT says: For SNP.test optioin, please input either use KI or use genotype") #if(is.null(file.path) & !byData & byFile) stop("GAPIT Ssays: A path for genotype data should be provided!") -if(is.null(file.total) & !byData & byFile) stop("GAPIT Ssays: Number of file should be provided: >=1") -if(!is.null(G) & !is.null(GD)) stop("GAPIT Ssays: Both hapmap and EMMA format exist, choose one only.") +if(is.null(file.total) & !byData & byFile) stop("GAPIT Says: Number of file should be provided: >=1") +if(!is.null(G) & !is.null(GD)) stop("GAPIT Says: Both hapmap and EMMA format exist, choose one only.") if(!is.null(file.GD) & is.null(file.GM) & (!is.null(GP)|!is.null(GK)) ) stop("GAPIT Ssays: Genotype data and map files should be in pair") if(is.null(file.GD) & !is.null(file.GM) & (!is.null(GP)|!is.null(GK)) ) stop("GAPIT Ssays: Genotype data and map files should be in pair") diff --git a/R/GAPIT.Genotype.View.R b/R/GAPIT.Genotype.View.R index bb3c9d7..602ea34 100644 --- a/R/GAPIT.Genotype.View.R +++ b/R/GAPIT.Genotype.View.R @@ -1,5 +1,5 @@ `GAPIT.Genotype.View` <-function(GI=NULL,X=NULL,chr=NULL, - WS0=NULL,ws=200,Aver.Dis=1000,...){ + WS0=NULL,ws=20,Aver.Dis=1000,...){ # Object: Analysis for Genotype data:Distribution of SNP density,Accumulation,Moving Average of density,result:a pdf of the scree plot # myG:Genotype data # chr: chromosome value @@ -132,9 +132,22 @@ for(i in 1:length(chr)) } odd=seq(1,length(chr),2) r=mapply(GAPIT.Cor.matrix,as.data.frame(x1),as.data.frame(x2)) +d.V=dist/Aver.Dis + +# fig.d=cbind(d.V[rs.index],(r^2)[rs.index]) +# dv=d.V[rs.index] +# r2=(r^2)[rs.index] +# dv2=ceiling(dv/ws)*ws +# dv2.un=as.character(unique(dv2)) +# fig.d=NULL +# for(i in 1:length(dv2.un)) +# { +# index=dv2==dv2.un[i] +# fig.d=rbind(fig.d,cbind(dv2.un[i],mean(r2[index],na.rm=TRUE))) +# } + grDevices::pdf("GAPIT.Genotype.Density_R_sqaure.pdf", width =10, height = 6) -d.V=dist/Aver.Dis # print(summary(d.V)) par(mfcol=c(2,3),mar = c(5,5,2,2)) plot(r[rs.index], xlab="Marker",las=1,xlim=c(1,mm), @@ -171,6 +184,7 @@ d.V.hist$counts=d.V0/d.V0.demo plot(d.V[rs.index],r[rs.index], las=1,xlab="Distance (Kb)", ylab="R", main="c",cex=.5,col="gray60",xlim=c(0,WS0/Aver.Dis)) abline(h=0,col="darkred") plot(d.V[rs.index],(r^2)[rs.index], las=1,xlab="Distance (Kb)", ylab="R sqaure", main="d",cex=.5,col="gray60",xlim=c(0,WS0/Aver.Dis)) +# plot(as.numeric(fig.d[,1]),as.numeric(fig.d[,2]), las=1,xlab="Distance (Kb)", ylab="R sqaure", main="d",cex=.5,col="gray60",xlim=c(0,WS0/Aver.Dis)) #Moving average dist2[dist2>WS0]=NA @@ -181,7 +195,7 @@ maPure=ma[!indRM,] maPure=maPure[!is.na(maPure[,1]),] ns=nrow(maPure) # ws=ws -slide=10 +slide=ws loc=matrix(NA,floor(ns/slide),2) for (i in 1:floor(ns/slide)){ diff --git a/R/GAPIT.IC.R b/R/GAPIT.IC.R index 4cde3fa..f4ae5b9 100644 --- a/R/GAPIT.IC.R +++ b/R/GAPIT.IC.R @@ -12,58 +12,68 @@ print("GAPIT.IC in process...") CV=DP$CV GD=DP$GD noCV=FALSE - if(is.null(CV)){ - noCV=TRUE - if(is.null(GD)) + if(is.null(CV)) { - CV=Y[,1:2] - }else{ - CV=GD[,1:2] - } - CV[,2]=1 - colnames(CV)=c("taxa","overall") - print(paste("There is 0 Covarinces.",sep="")) + noCV=TRUE + if(is.null(GD)) + { + CV=Y[,1:2] + }else{ + CV=GD[,1:2] + } + CV[,2]=1 + colnames(CV)=c("taxa","overall") + print(paste("There is 0 Covarinces.",sep="")) } Y=Y[!is.na(Y[,2]),] taxa_Y=as.character(Y[,1]) # print(head(Y)) if(DP$PCA.total>0&!is.null(DP$CV))CV=GAPIT.CVMergePC(DP$CV,PC) if(DP$PCA.total>0&is.null(DP$CV))CV=PC - if(is.null(GD)&!is.null(DP$KI)) + if(!is.null(GD)) { - taxa_KI=as.character(DP$KI[,1]) - taxa_CV=as.character(CV[,1]) - taxa_comall=intersect(intersect(taxa_KI,taxa_Y),taxa_CV) + if(!is.null(DP$KI)) + { + taxa_GD=as.character(GD[,1]) + taxa_KI=as.character(DP$KI[,1]) + taxa_CV=as.character(CV[,1]) + taxa_comall=intersect(intersect(intersect(taxa_KI,taxa_GD),taxa_Y),taxa_CV) # print(length(taxa_comall)) - 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=NULL + 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{ + }else{ # print("@@@@") - taxa_GD=as.character(GD[,1]) - taxa_comGD=as.character(GD[,1]) - taxa_CV=as.character(CV[,1]) - taxa_comall=intersect(intersect(taxa_GD,taxa_Y),taxa_CV) - comCV=CV[taxa_CV%in%taxa_comall,] - comCV <- comCV[match(taxa_comall,as.character(comCV[,1])),] - - comGD=GD[taxa_GD%in%taxa_comall,] - comGD <- comGD[match(taxa_comall,as.character(comGD[,1])),] - - comY=Y[taxa_Y%in%taxa_comall,] - comY <- comY[match(taxa_comall,as.character(comY[,1])),] + taxa_GD=as.character(GD[,1]) + taxa_comGD=as.character(GD[,1]) + taxa_CV=as.character(CV[,1]) + taxa_comall=intersect(intersect(taxa_GD,taxa_Y),taxa_CV) + comCV=CV[taxa_CV%in%taxa_comall,] + comCV <- comCV[match(taxa_comall,as.character(comCV[,1])),] + comGD=GD[taxa_GD%in%taxa_comall,] + comGD <- comGD[match(taxa_comall,as.character(comGD[,1])),] + comY=Y[taxa_Y%in%taxa_comall,] + comY <- comY[match(taxa_comall,as.character(comY[,1])),] + } + }else{ + # taxa_GD=as.character(GD[,1]) + taxa_KI=as.character(DP$KI[,1]) + taxa_CV=as.character(CV[,1]) + taxa_comall=intersect(intersect(taxa_KI,taxa_Y),taxa_CV) + # print(length(taxa_comall)) + 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=NULL } - - 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.") - print("The dimension of total CV is ") print(dim(comCV)) diff --git a/R/GAPIT.Main.R b/R/GAPIT.Main.R index c6baa5f..c172ad2 100644 --- a/R/GAPIT.Main.R +++ b/R/GAPIT.Main.R @@ -425,11 +425,13 @@ function(Y, # print(length(GK)) GTindex=match(as.character(KI[,1]),as.character(Y[,1])) GTindex=GTindex[!is.na(GTindex)] -# KI=KI[GTindex,GTindex+1] + # KI=KI[GTindex,GTindex+1] my_taxa=as.character(KI[,1]) CV=CV[as.character(CV[,1])%in%as.character(Y[,1]),] -# print(dim(KI)) -# print(dim(CV)) +print(dim(KI)) +print(dim(CV)) +print(dim(Y)) +print(length(GTindex)) #Output phenotype colnames(Y)=c("Taxa",name.of.trait) diff --git a/R/GAPIT.RandomModel.R b/R/GAPIT.RandomModel.R index d8d1b7c..f80a04e 100644 --- a/R/GAPIT.RandomModel.R +++ b/R/GAPIT.RandomModel.R @@ -29,16 +29,17 @@ function(GWAS,Y,CV=NULL,X,cutOff=0.01,GT=NULL,name.of.trait=NULL,N.sig=NULL,n_ra cutoff=max(sort.p[1:N.sig]) index=P.value<=cutoff } - geneGD=X[,index,drop=FALSE] - geneGWAS=GWAS[index,,drop=FALSE] - gene.licols=GAPIT.Licols(X=geneGD) - geneGD=gene.licols$Xsub - geneGWAS=geneGWAS[gene.licols$idx,] + if(length(unique(index))==1) { print("There is no significant marker for VE !!") return(list(GVs=NULL)) } + geneGD=X[,index,drop=FALSE] + geneGWAS=GWAS[index,,drop=FALSE] + gene.licols=GAPIT.Licols(X=geneGD) + geneGD=gene.licols$Xsub + geneGWAS=geneGWAS[gene.licols$idx,] index_T=as.matrix(table(index)) # print(index_T) in_True=ncol(geneGD)