Meta_SKAT.Work<-function(re, n.g, combined.weight=TRUE, n1=NULL, weights.beta=c(1,25),
method="davies", r.corr=0, is.separate=FALSE, Group_Idx=NULL, MAF.cutoff=1, Is.MAF=TRUE, missing_cutoff=0.15){
# optimal =optimal.mod
if(method=="optimal"){
method="optimal.mod"
}
# Combined MAF
p<-length(re[[1]]$MAF)
MAF.Combine=0
MAF.Groups<-list()
Map.Groups<-rep(0,n.g)
#re<<-re
MAF.list<-list()
for(i in 1:n.g){
MAF.list[[i]]<-re[[i]]$MAF
}
#MAF.list1<<-MAF.list
for(i in 1:n.g){
MAF.Combine = MAF.Combine + MAF.list[[i]] * n1[i] / sum(n1)
}
# If MAF.Combined==0 for all SNP, return p-value 1
if(sum(MAF.Combine) == 0){
warning("No polymorphic SNPs!",call.=FALSE)
return(list(p.value=1, p.value.resampling= NULL, pval.zero.msg=NULL))
}
# Get MAF.Groups when Group_Idx != NULL
ID.Groups = unique(Group_Idx)
for(j in 1:length(ID.Groups)){
MAF.Groups[[j]] = 0;
temp<-which(Group_Idx == ID.Groups[j])
Map.Groups[temp]<-j
for(i in temp){
MAF.Groups[[j]] = MAF.Groups[[j]] + MAF.list[[i]] * n1[i] / sum(n1[temp])
}
}
for(i in 1:n.g){
if(combined.weight == TRUE){
weight1<-Beta.Weights(MAF.Combine,weights.beta, MAF.cutoff, Is.MAF=Is.MAF)
} else {
j<-Map.Groups[i]
weight1<-Beta.Weights(MAF.Groups[[j]],weights.beta, MAF.cutoff, Is.MAF=Is.MAF)
}
# if missing < missing_cutoff
idx_missing_exclude<-which(re[[i]]$MissingRate > missing_cutoff)
if(length(idx_missing_exclude) > 0){
weight1[idx_missing_exclude]<-0
}
re[[i]]$Score = re[[i]]$Score * weight1
re[[i]]$SMat.Summary = t(t(re[[i]]$SMat.Summary * weight1) * weight1)
if(!is.null(re[[i]]$Score.Resampling)){
re[[i]]$Score.Resampling = re[[i]]$Score.Resampling * weight1
}
}
re.method<-SKAT_Check_Method(method,r.corr)
if(!is.separate){
re.score<-Meta_SKAT.Work.OneUnit(re, n.g)
} else {
re.score<-Meta_SKAT.Work.Groups(re, n.g, ID.Groups, Group_Idx)
}
re<-Met_SKAT_Get_Pvalue(re.score$Score, re.score$SMat.Summary, re.method$r.corr, re.method$method, re.score$Score.Resampling)
return(re)
}
####################################################
#
# Assume SMat and Setinfo are already aligned
MetaSKAT_withlist<-function(SMat.list, Info.list, n.cohort, n.each, combined.weight=TRUE, weights.beta=c(1,25),
method="davies", r.corr=0, is.separate = FALSE, Group_Idx=NULL, MAF.cutoff=1, missing_cutoff=0.15){
#Info.list1<<-Info.list
re<-list()
p<-length(Info.list[[1]]$MAF)
for(i in 1:n.cohort){
idx_miss<-which(is.na(Info.list[[i]]$MAF))
if(length(idx_miss) > 0){
Info.list[[i]]$Score[idx_miss] = 0
Info.list[[i]]$MAF[idx_miss] = 0
Info.list[[i]]$MissingRate[idx_miss] = 1
}
re1<-list( Score=Info.list[[i]]$Score, SMat.Summary = SMat.list[[i]], MAF=Info.list[[i]]$MAF, MissingRate=Info.list[[i]]$MissingRate)
re[[i]]<-re1
}
if(is.null(Group_Idx)){
Group_Idx<-1:n.cohort
}
# Use AlleleFreq2 instead of MAF, hence Is.MAF=FALSE
re = Meta_SKAT.Work(re, n.cohort, combined.weight, n1=n.each, weights.beta=weights.beta, method=method, r.corr=r.corr, is.separate=is.separate, Group_Idx=Group_Idx,
MAF.cutoff=MAF.cutoff, Is.MAF=TRUE, missing_cutoff=missing_cutoff )
return(re)
}
MetaSKAT_MSSD_OneSet<-function(Cohort.Info, SetID, combined.weight=TRUE, weights.beta=c(1,25),
method="davies", r.corr=0, is.separate = FALSE, Group_Idx=NULL, MAF.cutoff=1, missing_cutoff=0.15){
n.cohort = Cohort.Info$n.cohort
n.each=rep(0,n.cohort)
for(i in 1:n.cohort){
n.each[i]=Cohort.Info$EachInfo[[i]]$header$N.Sample
}
temp<-Get_META_Data_OneSet(Cohort.Info, SetID)
temp1<-Get_META_Data_OneSet_Align(temp$SMat.list, temp$Info.list, temp$IsExistSNV, n.cohort)
re<-MetaSKAT_withlist(temp1$SMat.list, temp1$Info.list, n.cohort, n.each, combined.weight=combined.weight,
weights.beta=weights.beta, method=method, r.corr= r.corr, is.separate = is.separate, Group_Idx=Group_Idx,
MAF.cutoff=MAF.cutoff, missing_cutoff=missing_cutoff)
return(re)
}
MetaSKAT_MSSD_ALL<-function(Cohort.Info, ...){
n.cohort = Cohort.Info$n.cohort
n.set<-length(Cohort.Info$Set_unique)
pval<-rep(NA,n.set)
for(i in 1:n.set){
SetID=Cohort.Info$Set_unique[i]
out<-try(MetaSKAT_MSSD_OneSet(Cohort.Info,SetID, ... ), silent = TRUE)
if(!Is_TryError(out)){
pval[i]<-MetaSKAT_MSSD_OneSet(Cohort.Info,SetID, ... )$p.value
}
}
re<-data.frame(SetID=Cohort.Info$Set_unique, p.value=pval)
return(re)
}
##################################################################
#
# Genotype matrix Z should be matched with y and X
#
# change from here
MetaSKAT_wZ<-function(Z, obj, combined.weight=TRUE, weights.beta=c(1,25),
method="davies", r.corr=0, is.separate = FALSE, Group_Idx=NULL, impute.method="fixed", impute.estimate.maf=1, missing_cutoff=0.15){
if(is.matrix(Z)!= TRUE){
stop("ERROR: Z is not a matrix!")
}
if(!Check_Class(obj, "META_NULL_Model") && !Check_Class(obj, "META_NULL_Model_EmmaX")){
stop("ERROR: obj class is not either META_NULL_Model or META_NULL_Model_EmmaX!")
}
IDX_MISS<-union(which(is.na(Z)),which(Z == 9))
if(length(IDX_MISS) > 0){
Z[IDX_MISS]<-NA
}
nSNP<-ncol(Z)
n<-nrow(Z)
n.g<-obj$n.g
# Calculate missing rate : Since Z is imputed before calling Meta_SKAT_SaveData when impute.estimate.maf==1,
# missing rate should be calculated before calling this function.
missing.ratio.Z = matrix(rep(0, nSNP*n.g), nrow=n.g)
for(i in 1:n.g){
ID<-obj$ID[[i]]
Z1<-as.matrix(Z[ID,])
n1<-nrow(Z1)
for(j in 1:nSNP){
missing.ratio.Z[i,j]<-length(which(is.na(Z1[,j])))/n1
}
}
Is.impute.cohortwise= FALSE
if(length(IDX_MISS) > 0){
msg<-sprintf("The missing genotype rate is %f. Imputation is applied.", (length(IDX_MISS))/length(Z) )
warning(msg,call.=FALSE)
if(impute.estimate.maf==1){
Is.impute.cohortwise = TRUE
} else if(impute.estimate.maf==2){
Z<-Impute(Z,impute.method=impute.method)
} else if(impute.estimate.maf==3){
if(is.null(Group_Idx)){
stop("ERROR: Group_Idx should not be NULL when impute.estimate.maf==3")
}
ID.Groups = unique(Group_Idx)
ID.all<-NULL
for(j in 1:length(ID.Groups)){
id.cohort.a<-which(Group_Idx == ID.Groups[j])
for(k in 1:length(id.cohort.a)){
id.cohort<-id.cohort.a[k]
ID<-obj$ID[[id.cohort]]
ID.all<-c(ID.all, ID)
}
Z[ID.all,]<-Impute(Z[ID.all,],impute.method=impute.method)
}
} else {
stop("ERROR: impute.estimate.mat is wrong! it should be either 1 or 2")
}
}
re1<-list()
for(i in 1:n.g){
ID<-obj$ID[[i]]
Z1<-as.matrix(Z[ID,])
re1[[i]]<-Meta_SKAT_SaveData(Z1, obj$out[[i]], SetID=NULL, impute.method = impute.method)
re1[[i]]$MissingRate = missing.ratio.Z[i,]
}
if(is.null(Group_Idx)){
Group_Idx<-1:n.g
}
re = Meta_SKAT.Work(re1, n.g, combined.weight, n1=obj$n.each, weights.beta=weights.beta, method=method, r.corr=r.corr,is.separate=is.separate
, Group_Idx=Group_Idx, missing_cutoff=missing_cutoff)
return(re)
}
MetaSKAT_wZ_OLD<-function(Z, obj, combined.weight=TRUE, weights.beta=c(1,25),
method="davies", r.corr=0, is.separate = FALSE, Group_Idx=NULL, impute.method="fixed", impute.estimate.maf=1, missing_cutoff=0.15){
if(is.matrix(Z)!= TRUE){
stop("ERROR: Z is not a matrix!")
}
if(!Check_Class(obj, "META_NULL_Model") && !Check_Class(obj, "META_NULL_Model_EmmaX")){
stop("ERROR: obj class is not either META_NULL_Model or META_NULL_Model_EmmaX!")
}
IDX_MISS<-union(which(is.na(Z)),which(Z == 9))
if(length(IDX_MISS) > 0){
Z[IDX_MISS]<-NA
}
nSNP<-ncol(Z)
n<-nrow(Z)
Is.impute.cohortwise= FALSE
missing.ratio.Z = rep(0, nSNP)
for(i in 1:nSNP){
missing.ratio.Z[i]<-length(which(is.na(Z[,i])))/n
}
if(length(IDX_MISS) > 0){
msg<-sprintf("The missing genotype rate is %f. Imputation is applied.", (length(IDX_MISS))/length(Z) )
warning(msg,call.=FALSE)
if(impute.estimate.maf==1){
Is.impute.cohortwise = TRUE
} else if(impute.estimate.maf==1){
Z<-Impute(Z,impute.method=impute.method)
} else {
stop("ERROR: impute.estimate.mat is wrong! it should be either 1 or 2")
}
}
n.g<-obj$n.g
re1<-list()
for(i in 1:n.g){
ID<-obj$ID[[i]]
if(Is.impute.cohortwise ){
Z1<-as.matrix(Z[ID,])
} else {
Z1<-Impute(as.matrix(Z[ID,]),impute.method=impute.method)
}
res<-obj$out[[i]]$res
res.out<-obj$out[[i]]$res.out
X1<-obj$out[[i]]$X1
if(obj$out_type=="C"){
s2<-obj$out[[i]]$s2
re1[[i]]<-Meta_SKAT_SaveData_Linear(res,Z1 , X1, s2, res.out)
} else if (obj$out_type=="D"){
pi_1<-obj$out[[i]]$pi_1
re1[[i]]<-Meta_SKAT_SaveData_Logistic(res,Z1 , X1, pi_1, res.out)
} else if (obj$out_type=="K"){
re1[[i]]<-Meta_SKAT_SaveData_Kinship(res,Z1 , obj$out[[i]]$P)
} else {
stop("ERROR: out_type is wrong!")
}
}
if(is.null(Group_Idx)){
Group_Idx<-1:n.g
}
re = Meta_SKAT.Work(re1, n.g, combined.weight, n1=obj$n.each, weights.beta=weights.beta, method=method, r.corr=r.corr,is.separate=is.separate
, Group_Idx=Group_Idx, missing_cutoff=missing_cutoff)
return(re)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.