Skip to content

Commit

Permalink
add matching taxa of KI in IC
Browse files Browse the repository at this point in the history
  • Loading branch information
jiabowang committed Sep 22, 2022
1 parent 89ff9a1 commit 12bfc1f
Show file tree
Hide file tree
Showing 6 changed files with 78 additions and 50 deletions.
1 change: 1 addition & 0 deletions R/GAPIT.Blink.R
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand Down
4 changes: 2 additions & 2 deletions R/GAPIT.Genotype.R
Original file line number Diff line number Diff line change
Expand Up @@ -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")
Expand Down
20 changes: 17 additions & 3 deletions R/GAPIT.Genotype.View.R
Original file line number Diff line number Diff line change
@@ -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
Expand Down Expand Up @@ -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),
Expand Down Expand Up @@ -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
Expand All @@ -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)){
Expand Down
84 changes: 47 additions & 37 deletions R/GAPIT.IC.R
Original file line number Diff line number Diff line change
Expand Up @@ -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))

Expand Down
8 changes: 5 additions & 3 deletions R/GAPIT.Main.R
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
11 changes: 6 additions & 5 deletions R/GAPIT.RandomModel.R
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down

0 comments on commit 12bfc1f

Please sign in to comment.