Nothing
MOJOV.read<-function(x=NULL,genofile=NULL,phenofile=NULL,indelfile=NULL,
header=T,column=1,...)
{
if(is.null(genofile)&class(x)!="MOJOV")
stop("You must input a correct genotype file name and path.")
if(is.null(phenofile))
stop("You must input a correct phenotype file name and path.")
if(class(x)=="MOJOV")
{
genoFile=genoFile(x)
genoIID=genoIID(x)
genoChr=genoChr(x)
genoPosi=genoPosi(x)
genoVariant=genoVariant(x)
genoH=genoH(x)
genoGene=genoGene(x)
}
else
{
x1<-read.csv(genofile,header=T,...)
names(x1)<-c("IID","chromosome","Position","SNP","Mut","Gene")
if(!is.null(indelfile)){
if(file.access(indelfile)>=0){
x2<-read.csv(indelfile,header=T,...)
genofile<-c(genofile,indelfile)
names(x2)<-c("IID","chromosome","Position","SNP","Mut","Gene")
x<-rbind(x1,x2)
}else{
warning("Your indel file has some error, it doesn't work.")
}
}else{
x=x1
}
g<-x[,5]
gsub("het","Het",g)
gsub("hom","Hom",g)
geno<-ifelse(g=="het",1,2)
geno<-as.numeric(geno)
genoFile=genofile
genoIID=as.character(x[,1])
genoChr=as.character(x[,2])
genoPosi=as.numeric(x[,3])
genoVariant=as.character(x[,4])
genoH=geno
genoGene=as.character(x[,6])
}
y<-read.csv(phenofile,header=header,...)
MOJOVObj<-new("MOJOV",genoFile=genoFile,genoIID=genoIID,
genoChr=genoChr,genoPosi=genoPosi,
genoVariant=genoVariant,genoH=genoH,
genoGene=genoGene,
phenoFile=phenofile,phenoIID=as.character(y[,1]),
phenoGender=as.numeric(y[,3]),phenoAge=as.numeric(y[,2]),
phenoHeight=as.numeric(y[,4]),phenoWeight=as.numeric(y[,5]),
phenoLabel=names(y)[5+column],phenoData=as.numeric(y[,5+column]),
phenoColumn=column,phenoIndNum=dim(y)[1])
return(MOJOVObj)
}
MOJOV.analysis<-function(x=NULL,MAF=0.01,ROI="scan",savefile=NULL,
codeMethod=c("Proportion","Indicator","ChuanhuaXin"),
weightMethod=NULL,
testMethod=c("FTest","WaldTest","LRT","Sandwich","all"),...)
{
if(class(x)!="MOJOV")
stop("Must specify a MOJOV object.")
codeMethod<-match.arg(codeMethod)
weightMethod<-match.arg(weightMethod)
testMethod<-match.arg(testMethod)
genoIID=genoIID(x)
genoVariant=genoVariant(x)
genoH=genoH(x)
genoGene=genoGene(x)
genoChr=genoChr(x)
genoPosi=genoPosi(x)
phenoIID=phenoIID(x)
phenoIndNum=phenoIndNum(x)
adjustData=adjustData(x)
varTot=unique(genoVariant)
varMAF=MAF
if(is.null(ROI))
{
tmp<-MOJOV.genoMatrix(genoIID=genoIID,genoVariant=genoVariant,
genoH=genoH,phenoIID=phenoIID,MAF=varMAF,...)
varFreq<-apply(tmp,2,sum)/(phenoIndNum*2)
tmp<-MOJOV.genoVector(x=tmp,y=adjustData,codeMethod=codeMethod,
weightMethod=weightMethod,...)
stat<-MOJOV.linearRegAnalysis(x=tmp,y=adjustData,
testMethod=testMethod,...)
resultPvalue=unlist(stat)
cat("P-value:",stat)
regionType="Whole"
regionChr=unique(genoChr)
regionPStart=numeric()
regionPStop=numeric()
regionGene=genoVariant
regionFile=character()
}
else
{
if(ROI%in%genoGene)
{
regionGene=ROI
index<-(1:length(genoGene))[genoGene==regionGene]
tmp<-MOJOV.genoMatrix(genoIID=genoIID[index],genoVariant=genoVariant[index],
genoH=genoH[index],phenoIID=phenoIID,MAF=varMAF,...)
varFreq<-apply(tmp,2,sum)/(phenoIndNum*2)
tmp<-MOJOV.genoVector(x=tmp,y=adjustData,codeMethod=codeMethod,
weightMethod=weightMethod,...)
return(lm(adjustData~tmp))
#This step has been edited. It will return a lm object directly.
stat<-MOJOV.linearRegAnalysis(x=tmp,y=adjustData,
testMethod=testMethod,...)
resultPvalue=unlist(stat)
cat("P-value:",stat)
regionType="Gene"
regionChr=genoChr[index][1]
regionPStart=min(genoPosi[index])
regionPStop=max(genoPosi[index])
regionFile=character()
}
else
{
if(file.exists(ROI))
{
regionFile=ROI
r<-read.csv(regionFile,header=FALSE,...)
if(dim(r)[2]!=4)
stop("ROI file has an incorrect format.")
pb<-txtProgressBar(min=1,max=length(r[,2]),style=3)
pbi=1
regionType="Position"
regionChr=as.character(r[,2])
regionPStart=r[,3]
regionPStop=r[,4]
regionGene=as.character(r[,1])
if(length((regionPStop>regionPStart)==FALSE)>0)
stop("Stop should greater than start position for ROI file.")
stat<-numeric()
for ( ii in 1:length(regionChr))
{
index<-genoChr==regionChr[ii]
IID_tmp<-genoIID[index]
Variant_tmp<-genoPosi[index]
H_tmp<-genoH[index]
Posi_tmp<-genoVariant[index]
index<-(Posi_tmp>regionPStart[ii]&Posi_tmp<regionPStart[ii])
IID_tmp<-IID_tmp[index]
Variant_tmp<-Variant_tmp[index]
H_tmp<-H_tmp[index]
index<-genoPosi[genoVariant==ii]
tmp<-MOJOV.genoMatrix(genoIID=IID_tmp,genoVariant=Variant_tmp,
genoH=H_tmp,phenoIID=phenoIID,MAF=varMAF,...)
varFreq<-apply(tmp,2,sum)/(phenoIndNum*2)
tmp<-MOJOV.genoVector(x=tmp,y=adjustData,codeMethod=codeMethod,
weightMethod=weightMethod,...)
stat<-rbind(stat,MOJOV.linearRegAnalysis(x=tmp,y=adjustData,
testMethod=testMethod,...))
setTxtProgressBar(pb,pbi)
pbi<-pbi+1
}
resultPvalue=stat
}
else
{
if(ROI=="scan")
{
geneTmp=unique(genoGene)
pb<-txtProgressBar(min=1,max=length(geneTmp),style=3)
pbi=1
stat<-numeric()
for( ROIT in geneTmp)
{
regionGene=ROIT
index<-(1:length(genoGene))[genoGene==regionGene]
tmp<-MOJOV.genoMatrix(genoIID=genoIID[index],genoVariant=genoVariant[index],
genoH=genoH[index],phenoIID=phenoIID,MAF=varMAF,...)
varFreq<-apply(tmp,2,sum)/(phenoIndNum*2)
tmp<-MOJOV.genoVector(x=tmp,y=adjustData,codeMethod=codeMethod,
weightMethod=weightMethod,...)
stat<-rbind(stat,MOJOV.linearRegAnalysis(x=tmp,y=adjustData,
testMethod=testMethod,...))
setTxtProgressBar(pb,pbi)
pbi<-pbi+1
}
resultPvalue=as.vector(stat)
regionType=geneTmp
regionChr=character()
regionPStart=numeric()
regionPStop=numeric()
regionFile=character()
}
else
stop("ROI should be a gene symbol from your file or a ROI file path.")
}
}
}
varRare<-varTot[varFreq<=varMAF]
if(!is.null(savefile))
{
stat<-matrix(stat,,9,byrow=T)
write.csv(x=data.frame(stat,phenotype=rep(phenoLabel(x),dim(as.matrix(stat))[1])),file=savefile,row.names=FALSE)
}
if(is.null(weightMethod))
weightMethod=character()
MOJOVObj<-new("MOJOV",varMAF=varMAF,varFreq=varFreq,varRare=varRare,
regionFile=regionFile,regionType=regionType,
regionChr=regionChr,regionPStart=regionPStart,
regionGene=regionGene,analyCode=codeMethod,
analyWeighted=weightMethod,
resultMethod=testMethod,resultPvalue=resultPvalue,
adjustAuto=adjustAuto(x),adjustData=adjustData(x),
adjustGender=adjustGender(x),adjustPower=adjustPower(x),
adjustPowerPvalue=adjustPowerPvalue(x),
adjustTerms=adjustTerms(x),
adjustPvalue=adjustPvalue(x),
genoFile=genoFile(x),genoIID=genoIID(x),
genoChr=genoChr(x),genoPosi=genoPosi(x),
genoVariant=genoVariant(x),genoH=genoH(x),
genoGene=genoGene(x),phenoFile=phenoFile(x),
phenoIID=phenoIID(x),phenoGender=phenoGender(x),phenoAge=phenoAge(x),
phenoHeight=phenoHeight(x),phenoWeight=phenoWeight(x),
phenoLabel=phenoLabel(x),phenoData=phenoData(x),
phenoColumn=phenoColumn(x),phenoIndNum=phenoIndNum(x))
return(MOJOVObj)
}
MOJOV.genoMatrix<-function(genoIID=NULL,genoVariant=NULL,
genoH=NULL,phenoIID=NULL,MAF=NULL)
{
if(is.null(genoIID))
stop("You should specify your genotype ID infomation.")
if(is.null(genoVariant))
stop("You should specify your genotype varaint lablels.")
if(is.null(genoH))
stop("You should specify your genotype varaint infomation.")
if(is.null(phenoIID))
stop("You should specify your phenotype ID infomation.")
variantLabels<-unique(genoVariant)
variantNum<-length(variantLabels)
individualsNum<-length(phenoIID)
I<-matrix(0,individualsNum,variantNum)
for( i in 1:length(genoH))
{
col=match(genoVariant[i],variantLabels)
row=match(genoIID[i],phenoIID)
I[row,col]<-genoH[i]
}
return(I)
}
MOJOV.genoVector<-function(x=NULL,y=NULL,
codeMethod=c("Proportion","Indicator","ChuanhuaXin"),
weightMethod=c("ChuanhuaXin"))
{
if(is.null(x)|is.null(y))
stop("You need specify your data.")
if(class(x)!="matrix")
stop("The data x should be a matrix.")
if(length(y)!=dim(x)[1])
stop("The genotype matrix x is not conformable.")
codeMethod<-match.arg(codeMethod)
if(codeMethod=="ChuanhuaXin")
{
weight<-MOJOV.weight(x=x,y=y,weightMethod=weightMethod)
X<-apply(x,1,function(x)sum(x*weight))
}
if(codeMethod=="Proportion")
{
X<-apply(x,1,sum)/(dim(x)[2]*2)
}
if(codeMethod=="Indicator")
{
X<-apply(x,1,sum)
X<-unlist(lapply(X,function(x)ifelse(x>0,1,0)))
}
return(X)
}
MOJOV.linearRegAnalysis<-function(x=NULL,y=NULL,
testMethod=c("FTest","WaldTest","LRT","Sandwich","all"),...)
{
if(is.null(x)|is.null(y))
stop("You need input correct data for x and y.")
if(dim(as.matrix(x))[2]!=1|dim(as.matrix(y))[2]!=1)
stop("The vector : X or Y is not one-dimension.")
if(length(x)!=length(y))
stop("X and Y is not same length.")
testMethod<-match.arg(testMethod)
if(testMethod=="FTest")
{
if(length(unique(x))==1)
return(1)
return(summary(lm(y~x))$coefficient[2,4])
}
if(testMethod=="WaldTest")
{
if(length(unique(x))==1)
return(1)
return(MOJOV.wald.test(x,y,...))
}
if(testMethod=="Sandwich")
{
if(length(unique(x))==1)
return(1)
return(MOJOV.saws(x,y,...))
}
if(testMethod=="SurveyWald")
{
if(length(unique(x))==1)
return(1)
return(MOJOV.regTermTest(x,y,...))
}
if(testMethod=="LRT")
{
if(length(unique(x))==1)
return(1)
return(MOJOV.regTermTest(x,y,method="LRT",...))
}
if(testMethod=="all")
{
if(length(unique(x))==1)
return(rep(1,9))
statistic.1<-MOJOV.wald.test(x,y,...)
statistic.2<-MOJOV.saws(x,y,...)
statistic.3<-MOJOV.regTermTest(x,y)
statistic.4<-MOJOV.regTermTest(x,y,method="LRT")
statistic<-unlist(c(wald.test.chi2=statistic.1$statistic[6]$result$chi2[1],
wald.test.p=statistic.1$statistic[6]$result$chi2[3],
saws.p1=statistic.2$statistic[9][1],
saws.p2=statistic.2$statistic[9][2],
regTT.wald.F=statistic.3$statistic[3],
regTT.wald.P=statistic.3$statistic[7],
regTT.LRT=statistic.4$statistic[3],
regTT.LRT=statistic.4$statistic[6],
pr=summary(statistic.1$fit)[12]$coefficient[2,4]
))
return(statistic)
}
}
MOJOV.phenotype<-function(x=NULL,auto=FALSE,gender=TRUE,power=1,Terms=1:3)
{
if(class(x)!="MOJOV")
stop("Must specify a MOJOV object.")
if(!is.logical(auto))
stop("auto should be logical.")
if(!is.logical(gender))
stop("gender should be logical.")
adjustAuto=auto
x_which<-phenoData(x)
phenoAge=phenoAge(x)
phenoHeight=phenoHeight(x)
phenoWeight=phenoWeight(x)
phenoGender=phenoGender(x)
x_x<-matrix(c(phenoAge,phenoHeight,phenoWeight),,3)
labels=c("Age","Height","Weight")
if(adjustAuto==T)
{
if(is.numeric(power))
{
if(length(power)<=1)
power<-c(seq(-5,5,0.01))
power<-power[power!=0]
pvalue_tmp<-numeric()
power_tmp<-numeric()
for ( i in power )
{
pvalue_tmp<-c(pvalue_tmp,unlist(shapiro.test(x_which^i)[2]))
power_tmp<-c(power_tmp,i)
}
pvalue_max<-pvalue_tmp[pvalue_tmp==max(pvalue_tmp)]
power_max<-power_tmp[pvalue_tmp==max(pvalue_tmp)]
x_which<-x_which^power_max
}
else
{
if(power!=FALSE)
power<-c(seq(-5,5,0.01))
if(power==F)
{
power_max=1
pvalue_max=unlist(shapiro.test(x_which)[2])
}
}
fit0<-glm(x_which~x_x)
p_vector<-summary(fit0)[12]$coefficients[(2:4),4]
Terms<-p_vector<=0.05
if(length(Terms[Terms==T]>0))
{
fit<-glm(x_which~x_x[,(1:3)[Terms]])
adjustTerms=names(x_x)[1:3][Terms]
adjustPvalue=unlist(summary(fit)[12]$coefficients[2:(length(adjustTerms)+1),4])
x_res<-fit$residuals
MOJOVObj<-new("MOJOV",adjustAuto=adjustAuto,adjustData=x_res,
adjustGender=FALSE,adjustPower=power_max,
adjustPowerPvalue=pvalue_max,
adjustTerms=labels[Terms],
adjustPvalue=adjustPvalue,
genoFile=genoFile(x),genoIID=genoIID(x),
genoChr=genoChr(x),genoPosi=genoPosi(x),
genoVariant=genoVariant(x),genoH=genoH(x),
genoGene=genoGene(x),phenoFile=phenoFile(x),
phenoIID=phenoIID(x),phenoGender=phenoGender(x),phenoAge=phenoAge(x),
phenoHeight=phenoHeight(x),phenoWeight=phenoWeight(x),
phenoLabel=phenoLabel(x),phenoData=phenoData(x),
phenoColumn=phenoColumn(x),phenoIndNum=phenoIndNum(x))
return(MOJOVObj)
}
else
{
MOJOVObj<-new("MOJOV",adjustAuto=adjustAuto,adjustData=x_which,
adjustGender=FALSE,adjustPower=power_max,
adjustPowerPvalue=pvalue_max,
adjustTerms="NA",adjustPvalue=Inf,
genoFile=genoFile(x),genoIID=genoIID(x),
genoChr=genoChr(x),genoPosi=genoPosi(x),
genoVariant=genoVariant(x),genoH=genoH(x),
genoGene=genoGene(x),phenoFile=phenoFile(x),
phenoIID=phenoIID(x),phenoGender=phenoGender(x),phenoAge=phenoAge(x),
phenoHeight=phenoHeight(x),phenoWeight=phenoWeight(x),
phenoLabel=phenoLabel(x),phenoData=phenoData(x),
phenoColumn=phenoColumn(x),phenoIndNum=phenoIndNum(x))
return(MOJOVObj)
}
}
else
{
if(mode(power)=="numeric"&length(power)==1)
{
x_which=x_which^power
}else{stop("power should be a number.")}
if(gender==TRUE)
{
#For primary sort
x_ID<-1:length(x_which)
Gender_1<-phenoGender==1
Gender_2<-phenoGender==2
x_1<-x_which[Gender_1]
x_ID_1<-x_ID[Gender_1]
x_2<-x_which[Gender_2]
x_ID_2<-x_ID[Gender_2]
if(mode(Terms)=="numeric")
{
x_fit_1<-glm(x_1~x_x[Gender_1,Terms])
x_fit_2<-glm(x_2~x_x[Gender_2,Terms])
x_res_1<-x_fit_1$residuals
x_res_2<-x_fit_2$residuals
adjustTerms=labels[Terms]
}
x_res<-c(x_res_1,x_res_2)[order(c(x_ID_1,x_ID_2))]
fit<-list(x_fit_1,x_fit_2)
adjustPvalue=c(unlist(summary(x_fit_1)[12]$coefficients[2:(length(adjustTerms)+1),4]),
unlist(summary(x_fit_2)[12]$coefficients[2:(length(adjustTerms)+1),4]))
names(adjustPvalue)<-c(paste(adjustTerms,"Gender1",sep=""),paste(adjustTerms,"Gender2",sep=""))
}
if(gender==FALSE)
{
fit<-glm(x_which~x_x[,Terms])
adjustPvalue=unlist(summary(fit)[12][2:(length(adjustTerms)+1),4])
x_res<-fit$residuals
}
if(power==F)power=1
adjustPowerPvalue=unlist(shapiro.test(x_which)[2])
MOJOVObj<-new("MOJOV",adjustAuto=adjustAuto,adjustData=x_res,
adjustGender=gender,adjustPower=power,
adjustPowerPvalue=adjustPowerPvalue,
adjustTerms=labels[Terms],
adjustPvalue=adjustPvalue,
genoFile=genoFile(x),genoIID=genoIID(x),
genoChr=genoChr(x),genoPosi=genoPosi(x),
genoVariant=genoVariant(x),genoH=genoH(x),
genoGene=genoGene(x),phenoFile=phenoFile(x),
phenoIID=phenoIID(x),phenoGender=phenoGender(x),phenoAge=phenoAge(x),
phenoHeight=phenoHeight(x),phenoWeight=phenoWeight(x),
phenoLabel=phenoLabel(x),phenoData=phenoData(x),
phenoColumn=phenoColumn(x),phenoIndNum=phenoIndNum(x))
return(MOJOVObj)
}
x_x<-as.matrix(x[,c(2,4,5)])
if(auto==TRUE)
{
if(mode(power)=="numeric")
{
if(length(power)<=1)
power<-c(seq(-5,5,0.01))
}
else
{
if(mode(power)=="logical"&length(power)==1)
{
if(power!=FALSE)
{
power<-c(seq(-5,5,0.01))
}
}
stop("power should be TRUE or FALSE or numeric vector.")
}
if(mode(power)=="numeric")
{
power<-power[power!=0]
pvalue_tmp<-numeric()
power_tmp<-numeric()
for ( i in power )
{
pvalue_tmp<-c(pvalue_tmp,unlist(shapiro.test(x_which^i)[2]))
power_tmp<-c(power_tmp,i)
}
pvalue_max<-pvalue_tmp[pvalue_tmp==max(pvalue_tmp)]
power_max<-power_tmp[pvalue_tmp==max(pvalue_tmp)]
x_which<-x_which^power_max
}
fit0<-glm(x[,6]~x_x)
p_vector<-summary(fit0)[12]$coefficients[(2:4),4]
Terms<-p_vector<=0.05
fit<-glm(x_which~x_x[,(1:3)[Terms]])
x_res<-fit$residuals
phenotype<-list(phenotype=x_res,fit=fit,fit0=fit0,
power=data.frame(power=power_max,pvalue=pvalue_max))
class(phenotype)<-"MOJOV.ad.Phenotype"
return(phenotype)
}else{
phenotype<-list(phenotype=x_res,fit=fit)
class(phenotype)<-"MOJOV.ad.Phenotype"
return(phenotype)
}
}
MOJOV.regTermTest<-function(x,y,...)
{
fit<-glm(y~x)
statistic<-regTermTest(fit,"x",...)
return(list(fit=fit,statistic=statistic))
}
MOJOV.saws<-function(x,y,...)
{
out<-lmfitSaws(model.matrix(~x),y)
statistic<-saws(out,...)
fit<-glm(y~x)
return(list(fit=fit,statistic=statistic))
}
#### MOJOV class definition ####
setClass(Class="MOJOV",
representation = representation(
genoFile="character",
genoIID="character",
genoChr="character",
genoPosi="numeric",
genoVariant="character",
genoH="numeric",
genoGene="character",
phenoFile="character",
phenoIID="character",
phenoGender="numeric",
phenoAge="numeric",
phenoHeight="numeric",
phenoWeight="numeric",
phenoLabel="character",
phenoData="numeric",
phenoColumn="numeric",
phenoIndNum="numeric",
adjustAuto="logical",
adjustData="numeric",
adjustGender="logical",
adjustPower="numeric",
adjustPowerPvalue="numeric",
adjustTerms="character",
adjustPvalue="numeric",
varMAF="numeric",
varFreq="numeric",
varTot="character",
varRare="character",
regionFile="character",
regionType="character",
regionChr="character",
regionPStart="numeric",
regionPStop="numeric",
regionGene="character",
analyCode="character",
analyWeighted="character",
resultMethod="character",
resultPvalue="numeric"
),
prototype=prototype()
)
#### Methods for extracting slots ###
# genotype file
setGeneric('genoFile', function(object) standardGeneric('genoFile'))
setMethod('genoFile','MOJOV', function(object) object@genoFile)
# genotype individuals ID
setGeneric('genoIID', function(object) standardGeneric('genoIID'))
setMethod('genoIID','MOJOV', function(object) object@genoIID)
# genotype chromosome
setGeneric('genoChr', function(object) standardGeneric('genoChr'))
setMethod('genoChr','MOJOV', function(object) object@genoChr)
# genotype position
setGeneric('genoPosi', function(object) standardGeneric('genoPosi'))
setMethod('genoPosi','MOJOV', function(object) object@genoPosi)
# genotype variants labels
setGeneric('genoVariant', function(object) standardGeneric('genoVariant'))
setMethod('genoVariant','MOJOV', function(object) object@genoVariant)
# genotype status
setGeneric('genoH', function(object) standardGeneric('genoH'))
setMethod('genoH','MOJOV', function(object) object@genoH)
# genotype gene labels
setGeneric('genoGene', function(object) standardGeneric('genoGene'))
setMethod('genoGene','MOJOV', function(object) object@genoGene)
# phenotype file
setGeneric('phenoFile', function(object) standardGeneric('phenoFile'))
setMethod('phenoFile','MOJOV', function(object) object@phenoFile)
# phenotype individuals ID
setGeneric('phenoIID', function(object) standardGeneric('phenoIID'))
setMethod('phenoIID','MOJOV', function(object) object@phenoIID)
# phenotype gender
setGeneric('phenoGender', function(object) standardGeneric('phenoGender'))
setMethod('phenoGender','MOJOV', function(object) object@phenoGender)
# phenotype age
setGeneric('phenoAge', function(object) standardGeneric('phenoAge'))
setMethod('phenoAge','MOJOV', function(object) object@phenoAge)
# phenotype height
setGeneric('phenoHeight', function(object) standardGeneric('phenoHeight'))
setMethod('phenoHeight','MOJOV', function(object) object@phenoHeight)
# phenotype weight
setGeneric('phenoWeight', function(object) standardGeneric('phenoWeight'))
setMethod('phenoWeight','MOJOV', function(object) object@phenoWeight)
# phenotype type label
setGeneric('phenoLabel', function(object) standardGeneric('phenoLabel'))
setMethod('phenoLabel','MOJOV', function(object) object@phenoLabel)
# phenotype data
setGeneric('phenoData', function(object) standardGeneric('phenoData'))
setMethod('phenoData','MOJOV', function(object) object@phenoData)
# which column of source phenotype file
setGeneric('phenoColumn', function(object) standardGeneric('phenoColumn'))
setMethod('phenoColumn','MOJOV', function(object) object@phenoColumn)
# the number of individuals from phenotype file
setGeneric('phenoIndNum', function(object) standardGeneric('phenoIndNum'))
setMethod('phenoIndNum','MOJOV', function(object) object@phenoIndNum)
# if auto choose the best parameters for adjusting phenotype
setGeneric('adjustAuto', function(object) standardGeneric('adjustAuto'))
setMethod('adjustAuto','MOJOV', function(object) object@adjustAuto)
# the phenotype data after adjusted
setGeneric('adjustData', function(object) standardGeneric('adjustData'))
setMethod('adjustData','MOJOV', function(object) object@adjustData)
# if adjusting the phenotype based different gender
setGeneric('adjustGender', function(object) standardGeneric('adjustGender'))
setMethod('adjustGender','MOJOV', function(object) object@adjustGender)
# the power for adjusting
setGeneric('adjustPower', function(object) standardGeneric('adjustPower'))
setMethod('adjustPower','MOJOV', function(object) object@adjustPower)
# the pvalue from shapiro test for the power
setGeneric('adjustPowerPvalue', function(object) standardGeneric('adjustPowerPvalue'))
setMethod('adjustPowerPvalue','MOJOV', function(object) object@adjustPowerPvalue)
# the terms for adjusting
setGeneric('adjustTerms', function(object) standardGeneric('adjustTerms'))
setMethod('adjustTerms','MOJOV', function(object) object@adjustTerms)
# the pvalue from t test for linear regression
setGeneric('adjustPvalue', function(object) standardGeneric('adjustPvalue'))
setMethod('adjustPvalue','MOJOV', function(object) object@adjustPvalue)
# the file name of region of interest
setGeneric('regionFile', function(object) standardGeneric('regionFile'))
setMethod('regionFile','MOJOV', function(object) object@regionFile)
# the type of region position or gene
setGeneric('regionType', function(object) standardGeneric('regionType'))
setMethod('regionType','MOJOV', function(object) object@regionType)
# the chromosome of ROI
setGeneric('regionChr', function(object) standardGeneric('regionChr'))
setMethod('regionChr','MOJOV', function(object) object@regionChr)
# the position start of ROI
setGeneric('regionPStart', function(object) standardGeneric('regionPStart'))
setMethod('regionPStart','MOJOV', function(object) object@regionPStart)
# the position stop of ROI
setGeneric('regionPStop', function(object) standardGeneric('regionPStop'))
setMethod('regionPStop','MOJOV', function(object) object@regionPStop)
# the labels of genes for ROI
setGeneric('regionGene', function(object) standardGeneric('regionGene'))
setMethod('regionGene','MOJOV', function(object) object@regionGene)
# the model for analysis
setGeneric('analyCode', function(object) standardGeneric('analyCode'))
setMethod('analyCode','MOJOV', function(object) object@analyCode)
# the method of weighted for model
setGeneric('analyWeighted', function(object) standardGeneric('analyWeighted'))
setMethod('analyWeighted','MOJOV', function(object) object@analyWeighted)
# the significance test for linear regression
setGeneric('resultMethod', function(object) standardGeneric('resultMethod'))
setMethod('resultMethod','MOJOV', function(object) object@resultMethod)
# the significance P value for linear regression
setGeneric('resultPvalue', function(object) standardGeneric('resultPvalue'))
setMethod('resultPvalue','MOJOV', function(object) object@resultPvalue)
### show ###
setMethod('show', signature='MOJOV', definition=function(object)
{
symbols1=paste(c(rep("#",40),"\n"),collapse="")
symbols2=paste(c(rep("=",40),"\n"),collapse="")
region<-1:(min(3,length(object@phenoIID)))
cat(symbols1)
cat(" MOJOV is a R package for rare variant\n")
cat(" Author:Conda E-mail:wkh983@sina.com\n")
cat(symbols2)
cat("Genotype Infomation\n")
cat(symbols2)
cat("File:")
cat(object@genoFile)
cat("\nIndividuals:")
cat(object@genoIID[region])
cat("... ...\nChromosome(s):")
region<-1:(min(3,length(unique(object@genoChr))))
cat(unique(object@genoChr)[region])
cat("\nGene Labels:")
region<-1:(min(3,length(unique(object@genoGene))))
cat(unique(object@genoGene)[region])
cat("\nVariant Labels:")
region<-1:(min(2,length(unique(object@genoVariant))))
cat(unique(object@genoVariant)[region])
cat("... ...\nPositions:")
region<-1:(min(3,length(object@genoPosi)))
cat(unique(object@genoPosi)[region])
cat(".\nGenotype Status:")
region<-1:(min(10,length(object@genoH)))
cat(object@genoH[region])
cat("... \n\n")
cat(symbols2)
cat("Phenotype Infomation\n")
cat(symbols2)
cat("File:")
cat(object@phenoFile)
cat("\nColumn:")
cat(object@phenoColumn)
cat("\nPhenotype Label:")
cat(object@phenoLabel)
cat("\n\n")
cat(symbols2)
cat("Correction Infomation\n")
if(length(object@adjustAuto)!=0)
{
if(object@adjustAuto==T)
cat("Auto correction\n")
else
cat("Non-Auto correction\n")
cat("Gender:")
cat(object@adjustGender)
cat(paste("\nExponent:",object@adjustPower," ",object@adjustPowerPvalue,sep=""))
cat("\nTerms:\n")
cat(object@adjustTerms)
cat("\n")
cat(object@adjustPvalue)
}
cat("\n\n")
cat(symbols2)
cat("Region of interest Infomation\n")
cat(symbols2)
if(length(object@regionType)>1)
{
cat("Type:")
cat(object@regionType[1:min(3,length(object@regionType))])
cat("...")
}
if(length(object@regionType)==1)
{
cat("Type:")
cat(object@regionType)
if(object@regionType=="Single")
{
cat("\nGene:")
cat(object@regionGene)
}
if(object@regionType=="Gene")
{
cat("File:")
cat(object@regionFile)
cat("\nGene:")
region<-1:(min(3,length(object@regionGene)))
cat(object@regionGene[region])
}
if(object@regionType=="Position")
{
cat("File:")
cat(object@regionFile)
cat("\nChromosome:")
region<-1:(min(3,length(object@regionChr)))
cat(object@regionChr[region])
cat("\nPosition Start:")
cat(object@regionPStart[region])
cat("\nPosition Stop:")
cat(object@regionPStopt[region])
cat("\nGene:")
region<-1:(min(3,length(object@regionGene)))
cat(object@regionGene[region])
}
}
cat("\n\n")
cat(symbols2)
cat("Result after analysis\n")
cat(symbols2)
cat("Model:")
cat(object@analyCode)
cat("\nWeight Method:")
cat(object@analyWeighted)
cat("\nSignificance Test Methods:")
cat(object@resultMethod)
cat("\nSignificance Test Result:\n")
region<-1:(min(3,length(object@resultPvalue)))
cat(object@resultPvalue[region])
cat("...\n")
cat(symbols1)
}
)
MOJOV.simulation<-function(cohortSize=500,nReps=2,theta=10,sites=NULL,affectNum=NULL,
MAF=0.01,totalMAF=0.05,lambda=1,sd=NULL,type = c("alpha","belta"),
sampleNum=100,outFile=NULL,plot=FALSE,
codeMethod=c("Proportion","Indicator","ChuanhuaXin"),
weightMethod=c("ChuanhuaXin"),
testMethod=c("FTest","WaldTest","LRT","Sandwich","all"),
save=NULL)
{
type<-match.arg(type)
if(!is.null(sites))
theta<-paste(theta,"-s",sites)
if(is.null(outFile))
outFile<-"MOJOV.simulation.out"
if(sampleNum>cohortSize)
stop("The sample number should not greater than population.")
cohortSize<-cohortSize*2
commandCharacter<-paste("ms",cohortSize,nReps,"-t",theta,">",outFile)
cat("ms is running...\nIt may take lots of time, please wait for it patiently...\n")
system(commandCharacter)
inFile<-file(outFile,"r")
systemInfo<-readLines(inFile,n=2)
pb<-txtProgressBar(min=0,max=nReps,style=3)
#data collection
result<-numeric()
resultSummary<-numeric()
for ( i in 1:nReps)
{
garbage<-readLines(inFile,n=2)
if(!is.null(sites))
garbage<-readLines(inFile,n=1)
segsites<-as.numeric(unlist(strsplit(readLines(inFile,n=1),split=" "))[2])
positions<-as.numeric(unlist(strsplit(readLines(inFile,n=1),split=" "))[-1])
tmp<-readLines(inFile,n=cohortSize)
tmp<-as.numeric(unlist(strsplit(tmp,split="")))
x<-matrix(tmp,cohortSize,segsites,byrow=T)
sitesLabel<-paste("SITE",1:segsites,sep="")
#pairing together to form an individual at random
sampleIndex<-sample(1:cohortSize,cohortSize)
sampleIndex<-matrix(sampleIndex,,2)
sample<-unlist(lapply(1:dim(sampleIndex)[1],function(i)apply(x[sampleIndex[i,],],2,sum)))
sample<-matrix(sample,,segsites,byrow=T)
dimnames(sample)<-list(paste("ID",(1:(cohortSize/2)),sep=""),sitesLabel)
#Generating phenotype vector
freq<-apply(sample,2,sum)/cohortSize
resultSummary<-rbind(resultSummary,summary.freq(freq))
if(is.null(sd))
stdd<-sd(freq)*lambda
else
stdd<-sd
if(type=="alpha")
{
index_tmp<-(1:segsites)[freq<=MAF]
index<-numeric()
if(is.null(affectNum))
{
total<-0
for(ii in sample(index_tmp,length(index_tmp)))
{
total=total+freq[ii]
index=c(index,ii)
if(total>totalMAF)break
}
}
else
{
if(affectNum<=1)
affect<-floor(affectNum*length(index_tmp))
else
affect<-affectNum
index<-sample(index_tmp,affect)
}
if(length(index)<2)
{
warning("The affect sites is less than 2.")
next
}
rareSample<-sample[,index]
meanPhenotype<-apply(rareSample,1,sum)
phenotype<-unlist(lapply(meanPhenotype,function(x)rnorm(1,x,stdd)))
}
else
phenotype<-rnorm((cohortSize/2),0,stdd)
#sample at random
index<-sample((1:(cohortSize/2)),sampleNum)
x<-sample[index,]
y<-phenotype[index]
result<-rbind(result,unlist(MOJOV.linearRegAnalysis(x=MOJOV.genoVector(x,y,
weightMethod=weightMethod,
codeMethod=codeMethod),y=y,
testMethod=testMethod)))
# return(list(sample=sample,rareSample=rareSample,phenotype=phenotype,
# x=x,y=y,index=index,result=result))
setTxtProgressBar(pb,i)
}
apply(as.matrix(result),2,mean)
cat("\n")
if(!is.null(save))
write.csv(x=result,file=save,row.names=FALSE)
close(inFile)
if(is.null(outFile))
file.remove("MOJOV.simulation.out")
Summary=apply(as.matrix(resultSummary),2,mean)
if(plot==TRUE)
barchart(Summary,col="pink",horizontal=T,xlab="Mean",ylab="Frequency",main="Distribution of Frequency")
list(TestResult=result,
SummaryResult=resultSummary,
Summary=Summary)
}
MOJOV.wald.test<-function(x,y,...)
{
fit<-glm(y~x)
statistic<-wald.test(b=coef(fit),Sigma=vcov(fit),Terms=2,...)
return(list(fit=fit,statistic=statistic))
}
MOJOV.weight<-function(x=NULL,y=NULL,weightMethod="ChuanhuaXin")
{
weightMethodVector<-c("ChuanhuaXin")
if(length(weightMethod)!=1&is.character(weightMethod))
stop("weightMethod should be one character.")
if(!weightMethod%in%weightMethodVector)
stop("weightMethod is not proper.")
sitesNum<-dim(x)[2]
individualsNum<-dim(x)[1]
if(weightMethod == "ChuanhuaXin")
{
delta<-ifelse(y<mean(y)+sd(y)&y<mean(y)-sd(y),1,0)
weight<-unlist(lapply(1:sitesNum,function(col.index)
{
p<-(sum(x[,col.index]*delta)+1)/(2*(sum(delta)+1))
1/(2*individualsNum*p*(1-p))^(1/2)
}))
return(weight)
}
if(weightMethod == "")
{
weight<-rep(1,individualsNum)
return(weight)
}
}
MOJOV.export<-function(x=NULL,file="MOJOV.result.csv",p=0.001)
{
if(class(x)!="MOJOV")
stop("Must specify a MOJOV object.")
regionType<-regionType(x)
resultPvalue<-resultPvalue(x)
resultMethod<-resultMethod(x)
analyCode<-analyCode(x)
genoFile<-genoFile(x)
phenoFile<-phenoFile(x)
phenoLabel<-phenoLabel(x)
if(length(regionType)>1)
{
note1<-paste("##Genotype File :",genoFile)
note2<-paste("##Phenotype File :",phenoFile)
note3<-paste("##Phenotype Name :",phenoLabel)
note4<-paste("##Code Method :",analyCode)
note5<-paste("##Test Method :",resultMethod)
note6<-paste("##Time Stamp :", Sys.time())
note<-c("##This file is exported by MOJOV package.",
note1,note2,note3,note4,note5,note6)
write.table(x=data.frame(note),row.names=F,col.names=F,
file=file,quote=FALSE)
dataSet<-data.frame(Gene=regionType,pvalue=resultPvalue)
write.table(x=dataSet,row.names=F,append=T,file=file,sep=",")
}
return(x@regionType[x@resultPvalue<=p])
}
summary.freq<-function(x, bin=c(0,0.01,.05,.1,.2,.3,.4,.5))
{
len<-length(bin)-1
total<-length(x)
labels<-character()
result<-numeric()
x[x>0.5]<-1-x[x>0.5]
for(i in 1:len)
{
result<-c(result,length(x[x>bin[i]&x<=bin[i+1]]))
labels<-c(labels,paste("[",bin[i],",",bin[i+1],"]",sep=""))
}
result<-result/total
names(result)<-labels
result
}
MOJOV.summary<-function(x=NULL,ROI="scan",bin=c(0,0.01,.05,.1,.2,.3,.4,.5),
MAF=0.05,...)
{
if(class(x)!="MOJOV")
stop("Must specify a MOJOV object.")
genoIID=genoIID(x)
genoChr=genoChr(x)
genoPosi=genoPosi(x)
genoVariant=genoVariant(x)
genoH=genoH(x)
genoGene=genoGene(x)
phenoFile=phenoFile(x)
phenoIID=phenoIID(x)
phenoGender=phenoGender(x)
phenoAge=phenoAge(x)
phenoHeight=phenoHeight(x)
phenoWeight=phenoWeight(x)
phenoLabel=phenoLabel(x)
phenoData=phenoData(x)
phenoColumn=phenoColumn(x)
phenoIndNum=phenoIndNum(x)
if(ROI=="scan")
{
geneList=unique(genoGene)
pb<-txtProgressBar(min=1,max=length(geneList),style=3)
pbi=1
sumFreq<-numeric()
geneVarNum<-numeric()
geneRareNum<-numeric()
for( i in geneList)
{
regionGene=i
index<-(1:length(genoGene))[genoGene==regionGene]
tmp<-MOJOV.genoMatrix(genoIID=genoIID[index],genoVariant=genoVariant[index],
genoH=genoH[index],phenoIID=phenoIID,MAF=MAF,...)
varFreq<-apply(tmp,2,sum)/(phenoIndNum*2)
sumFreq<-rbind(sumFreq,summary.freq(varFreq,bin=bin))
geneVarNum<-c(geneVarNum,length(varFreq))
geneRareNum<-c(geneRareNum,length(varFreq[varFreq<=MAF]))
setTxtProgressBar(pb,pbi)
pbi<-pbi+1
}
bin=names(summary.freq(varFreq,bin=bin))
#total information
totalVariant<-sum(geneVarNum)
totalRareVar<-sum(geneRareNum)
totalGene<-length(geneList)
#gene total
geneRatio<-geneRareNum/geneVarNum
geneDist<-apply(sumFreq,2,mean)
#gene variants
geneVarTmpValue<-geneVarNum[order(geneVarNum,decreasing=T)]
geneVarTmpGene<-geneList[order(geneVarNum,decreasing=T)]
if(geneVarTmpGene[1]=="NULL")
{
geneVarTmpValue<-geneVarTmpValue[-1]
geneVarTmpGene<-geneVarTmpGene[-1]
}
geneVarMaxValue<-geneVarTmpValue[1]
geneVarMaxGene<-geneVarTmpGene[1]
geneVarMinValue<-geneVarTmpValue[totalGene]
geneVarMinNum<-length(geneVarTmpValue[geneVarTmpValue==geneVarMinValue])
geneVarMeanValue<-mean(geneVarTmpValue)
#gene rare variants
geneRareTmpValue<-geneRareNum[order(geneRareNum,decreasing=T)]
geneRareTmpGene<-geneList[order(geneRareNum,decreasing=T)]
if(geneRareTmpGene[1]=="NULL")
{
geneRareTmpValue<-geneRareTmpValue[-1]
geneRareTmpGene<-geneRareTmpGene[-1]
}
geneRareMaxValue<-geneRareTmpValue[1]
geneRareMaxGene<-geneRareTmpGene[1]
geneRareMinValue<-geneRareTmpValue[length(geneRareTmpValue)]
geneRareMinNum<-length(geneRareTmpValue[geneRareTmpValue==geneRareMinValue])
geneRareMeanValue<-mean(geneRareTmpValue)
#gene ratio
geneRatioMaxValue<-max(geneRatio)
geneRatioMinValue<-min(geneRatio)
geneRatioMaxNum<-length(geneRatio[geneRatio==geneRatioMaxValue])
geneRatioMinNum<-length(geneRatio[geneRatio==geneRatioMinValue])
geneRatioMeanValue<-mean(geneRatio)
MOJOV.SummaryObject<-new("MOJOV.Summary",
ROI=ROI,
MAF=MAF,
bin=bin,
#total information
totalGene=totalGene,
totalVariant=totalVariant,
totalRareVar=totalRareVar,
#gene total
geneVarNum=geneVarNum,
geneRareNum=geneRareNum,
geneList=geneList,
geneDist=geneDist,
geneRatio=geneRatio,
#gene variants
geneVarMaxValue=geneVarMaxValue,
geneVarMaxGene=geneVarMaxGene,
geneVarMinValue=geneVarMinValue,
geneVarMinNum=geneVarMinNum,
geneVarMeanValue=geneVarMeanValue,
#gene rare variants
geneRareMaxValue=geneRareMaxValue,
geneRareMaxGene=geneRareMaxGene,
geneRareMinValue=geneRareMinValue,
geneRareMinNum=geneRareMinNum,
geneRareMeanValue=geneRareMeanValue,
#gene ratio
geneRatioMaxValue=geneRatioMaxValue,
geneRatioMaxNum=geneRatioMaxNum,
geneRatioMinValue=geneRatioMinValue,
geneRatioMinNum=geneRatioMinNum,
geneRatioMeanValue=geneRatioMeanValue
)
return(MOJOV.SummaryObject)
}
else
{
if(ROI%in%genoGene)
{
index<-(1:length(genoGene))[genoGene==regionGene]
tmp<-MOJOV.genoMatrix(genoIID=genoIID[index],genoVariant=genoVariant[index],
genoH=genoH[index],phenoIID=phenoIID,MAF=MAF,...)
varFreq<-apply(tmp,2,sum)/(phenoIndNum*2)
geneDist<-summary.freq(varFreq,bin=bin)
geneVarNum<-length(varFreq)
geneRareNum<-length(varFreq[varFreq<=MAF])
geneRatio<-geneRareNum/geneVarNum
bin<-names(geneDist(varFreq,bin=bin))
MOJOV.SummaryObject<-new("MOJOV.Summary",
ROI=ROI,
MAF=MAF,
bin=bin,
#total information
totalGene=numeric(),
totalVariant=numeric(),
totalRareVar=numeric(),
#gene total
geneVarNum=geneVarNum,
geneRareNum=geneRareNum,
geneList="character",
geneDist=geneDist,
geneRatio=geneRatio,
#gene variants
geneVarMaxValue=numeric(),
geneVarMaxGene=numeric(),
geneVarMinValue=numeric(),
geneVarMinNum=numeric(),
geneVarMeanValue=numeric(),
#gene rare variants
geneRareMaxValue=numeric(),
geneRareMaxGene=numeric(),
geneRareMinValue=numeric(),
geneRareMinNum=numeric(),
geneRareMeanValue=numeric(),
#gene ratio
geneRatioMaxValue=numeric(),
geneRatioMaxNum=numeric(),
geneRatioMinValue=numeric(),
geneRatioMinNum=numeric(),
geneRatioMeanValue=numeric()
)
return(MOJOV.SummaryObject)
}
else
stop("ROI should be a gene symbol from your file or \"scan\".")
}
}
setClass(Class="MOJOV.Summary",
representation=representation(
ROI="character",
MAF="numeric",
bin="character",
#total information
totalGene="numeric",
totalVariant="numeric",
totalRareVar="numeric",
#gene total
geneVarNum="numeric",
geneRareNum="numeric",
geneList="character",
geneDist="numeric",
geneRatio="numeric",
#gene variants
geneVarMaxValue="numeric",
geneVarMaxGene="character",
geneVarMinValue="numeric",
geneVarMinNum="numeric",
geneVarMeanValue="numeric",
#gene rare variants
geneRareMaxValue="numeric",
geneRareMaxGene="character",
geneRareMinValue="numeric",
geneRareMinNum="numeric",
geneRareMeanValue="numeric",
#gene ratio
geneRatioMaxValue="numeric",
geneRatioMaxNum="numeric",
geneRatioMinValue="numeric",
geneRatioMinNum="numeric",
geneRatioMeanValue="numeric"
),
prototype=prototype()
)
setMethod('show', signature="MOJOV.Summary",
definition=function(object)
{
if(object@ROI=="scan")
{
symbols1=paste(c(rep("#",40),"\n"),collapse="")
cat(symbols1)
cat("Summary statistic infomations")
cat("\nGene :", object@ROI)
cat("\nThreshold value :",object@MAF)
cat("\nThe total number of variants :",object@totalVariant)
cat("\nThe total number of rare variants :",object@totalRareVar)
cat("\nThe total number of gene :",object@totalGene)
cat("\nThe average distribution for all genes :\n")
for(i in object@bin)
{
cat(i,":",object@geneDist[object@bin==i],"\n")
}
cat("\nInformation for one gene.")
cat("\nAverage number of variants :",object@geneVarMeanValue)
cat("\nMax number of variants :",object@geneVarMaxGene,
object@geneVarMinValue)
cat("\nMin number of variants :",object@geneVarMinValue,
" COUNTS :",object@geneVarMinNum)
cat("\nAverage number of rare variants :",object@geneRareMeanValue)
cat("\nMax number of rare variants :",object@geneRareMaxGene,
object@geneRareMinValue)
cat("\nMin number of rare variants :",object@geneRareMinValue,
" COUNTS :",object@geneRareMinNum)
cat("\nAverage ratio for all genes :",object@geneRatioMeanValue)
cat("\nMax ratio for all genes :",object@geneRatioMaxValue,
" COUNTS :",object@geneRatioMaxNum)
cat("\nMin ratio for all genes :",object@geneRatioMinValue,
" COUNTS :",object@geneRatioMinNum)
}
else
{
symbols1=paste(c(rep("#",40),"\n"),collapse="")
cat(symbols1)
cat("Summary statistic infomations")
cat("\nGene :", object@ROI)
cat("\nThreshold value :",object@MAF)
cat("\nThe number of variants :",object@geneVarNum)
cat("\nThe number of rare variants :",object@geneRareNum)
cat("\nThe ratio :",object@geneRatio)
cat("\nThe average distribution for all genes :\n")
for(i in object@bin)
{
cat(i,":",object@geneDist[object@bin==i],"\n")
}
}
})
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.