R/X_Chrom.R

In SKAT: SNP-Set (Sequence) Kernel Association Test

Impute_XChr<-function(Z, impute.method, SexVar){
	
	p<-dim(Z)[2]
	MAF<-SKAT_Get_MAF(Z, id_include=NULL, Is.chrX=TRUE, SexVar=SexVar)
	
	id.male<-which(SexVar==1)
	id.female<-which(SexVar==2)
	
	for(i in 1:p){
		IDX<-which(is.na(Z[,i]))
		IDX.m<-intersect(IDX, id.male)
		IDX.f<-intersect(IDX, id.female)
		
		maf1<-MAF[i]
		
		if(impute.method =="fixed"){
			if(length(IDX.m) > 0){
				Z[IDX.m,i]<- maf1
			}
			if(length(IDX.f) > 0){
				Z[IDX.f,i]<-2 * maf1
			}
		} else if(impute.method =="bestguess"){
			if(length(IDX.m) > 0){
				Z[IDX.m,i]<- round(maf1)
			}
			if(length(IDX.f) > 0){
				Z[IDX.f,i]<-round(2 * maf1)
			}
		
		} else {
			stop("Error: Imputation method shoud be \"fixed\" or \"bestguess\" ")
		}
		
	}

	return(Z)
}

#
#	X.inact : x inactivation
#
SKAT_Code_XChr<-function(Z, is_X.inact, SexVar){
	
	id.male<-which(SexVar==1)
	if(is_X.inact && length(id.male) > 0){
		Z[id.male,] = 2* Z[id.male,]
	}
	return(cbind(Z))

}



#
# Return Sex variable
#

Check_Sex_Var = function(formula, SexVar, data=NULL){

	# Check SexVar (1=male, 2=female)
	temp.dat<-get_all_vars(formula, data)
	idx.var.sex = which(colnames(temp.dat) == SexVar)
	if(length(idx.var.sex) == 0){
		msg<-sprintf("No variable %s exists in the formula!\n", SexVar);
		stop(msg)
	}
	
	# check missing
	Sex.Var<-temp.dat[,idx.var.sex]
	Sex.Var.Org<-Sex.Var
	id.missing<-which(is.na(Sex.Var))
	if(length(id.missing) == length(Sex.Var)){
		msg<-sprintf("All sex values are NA\n");
		stop(msg);
	} else if(length(id.missing) > 0){
		msg<-sprintf("There are %d missing values in the Sex variable\n", length(id.missing));
		cat(msg);
		Sex.Var<-Sex.Var[id.missing]
	}
	
	# check whether SexVar is 1, 2 variable
	
	id<-intersect(which(Sex.Var!=1), which(Sex.Var!=2))
	if(length(id) >0){
		msg<-sprintf("%d subjects have sex values neither 1 nor 2. They should be either 1 (=male) or 2(=female)!", length(id))
		stop(msg);
	}
	
	return(Sex.Var.Org)
}



SKAT_Null_Model_ChrX = function(formula, SexVar, data=NULL, out_type="C", n.Resampling=0, type.Resampling="bootstrap", Adjustment=TRUE, Model.Y=FALSE){
	
	SKAT_MAIN_Check_OutType(out_type)
	SexVar=Check_Sex_Var(formula, SexVar, data=data)

	
	# check missing 
	obj1<-model.frame(formula,na.action = na.omit,data)
	obj2<-model.frame(formula,na.action = na.pass,data)

	n<-dim(obj2)[1]
	n1<-dim(obj1)[1]
	id_include<-SKAT_Null_Model_Get_Includes(obj1,obj2)
	n1<-length(id_include)
	
	# Check whether n < 2000 and out_type="D", apply the adjustment 
	# if No_Adjustment = FALSE
	if(n1< 2000 && out_type=="D" && Adjustment){
		MSG<-sprintf("Sample size (non-missing y and X) = %d, which is < 2000. The small sample adjustment is applied!\n",n )
		cat(MSG)
		n.Resampling.kurtosis=10000
		re<-SKAT_Null_Model_MomentAdjust(formula, data, n.Resampling, type.Resampling=type.Resampling, is_kurtosis_adj=TRUE, n.Resampling.kurtosis=n.Resampling.kurtosis)
		
		class(re)<-"SKAT_NULL_Model_Adj_ChrX"
		re$SexVar=SexVar
		re$n.all<-n
		return(re)
	}


	if(n - n1 > 0){
		MSG<-sprintf("%d  samples have either missing phenotype or missing covariates. They are excluded from the analysis!",n - n1)
		warning(MSG,call.=FALSE)
	}

	if(out_type=="C"){
		re<-Get_SKAT_Residuals.linear(formula, data, n.Resampling, type.Resampling, id_include )
	} else {
		re<-Get_SKAT_Residuals.logistic (formula, data, n.Resampling, type.Resampling, id_include )
	}
	
	obj.Y=NULL
	if(Model.Y==TRUE){
	  # Run with Male only
	  id.male<-which(SexVar==1)
	  
	  if(length(id.male)==0){
	    stop("No male subject!!") 
	  }
	  
	  datY<-model.frame(formula, na.action = na.pass,data)[id.male,]
	  obj.Y = SKAT_Null_Model(formula, data=datY, 
	                              out_type=out_type, n.Resampling=n.Resampling, type.Resampling=type.Resampling, Adjustment=Adjustment )
	  
	  
	  #if(Check_Class(obj.Y, "SKAT_NULL_Model_ADJ")){
	  #  id_include1 = obj.Y$re1$id_include
	  #} else if(Check_Class(obj.Y, "SKAT_NULL_Model")){
	  #  id_include1 = obj.Y$id_include
	  #} else {
	  #  msg<-sprintf("Class error %s", class(obj.Y))
	  #  stop(msg)
	  #}
	  
	  #check number of individuals
	  #if(length(id_include1) != length(subset_male) ){
	  #  msg<-sprintf("ChrY phenotype processing, numbers do not match[%d], [%d]", length(id_include1),length(subset_male) )
	  #  stop(msg) 
	  #}
	 
	  obj.Y$id.male = id.male
	}


	class(re)<-"SKAT_NULL_Model_ChrX"
	re$SexVar=SexVar
	re$n.all<-n
	re$obj.Y = obj.Y
	return(re)
	
}



SKAT_ChrX<-function(Z, obj, is_X.inact =TRUE, kernel = "linear.weighted", method="davies", weights.beta=c(1,25)
, weights = NULL, impute.method = "fixed", r.corr=0, is_check_genotype=TRUE
, is_dosage = FALSE, missing_cutoff=0.15, max_maf=1, estimate_MAF=1, SetID=NULL){

	
	if(kernel != "linear" && kernel != "linear.weighted"){

		if(Check_Class(obj,  "SKAT_NULL_Model_Adj_ChrX")){
			msg<-sprintf("The small sample adjustment only can be applied for linear and linear.weighted kernel in the current version of SKAT! No adjustment is applied")
			warning(msg,call.=FALSE)
			obj<-obj$re1
		}

	}
	
	if(!Check_Class(obj, "SKAT_NULL_Model_Adj_ChrX") && !Check_Class(obj, "SKAT_NULL_Model_ChrX")){
		msg<-sprintf("The obj parameter should be a returned object from SKAT_Null_Model_ChrX!")
		stop(msg)
	}
	

	SexVar = obj$SexVar
	n.all=obj$n.all

 	if(Check_Class(obj, "SKAT_NULL_Model_Adj_ChrX")){
		
		id_include = obj$re1$id_include
	

	} else if(Check_Class(obj, "SKAT_NULL_Model_ChrX")){

		id_include = obj$id_include

	} else {
		#re<-SKAT_MAIN(Z,obj, ...)
		stop("The old interface is defunct! Please run SKAT_NULL_Model first!")
	}	

	out.z<-SKAT_MAIN_Check_Z(Z, n.all, id_include, SetID
		, weights, weights.beta, impute.method, is_check_genotype
		, is_dosage, missing_cutoff,max_maf=max_maf, estimate_MAF=estimate_MAF, Is.chrX=TRUE, SexVar=SexVar)

	out.z$Z.test<-SKAT_Code_XChr(out.z$Z.test, is_X.inact=is_X.inact, SexVar=SexVar[out.z$id_include.test])
	
	
 	if(Check_Class(obj, "SKAT_NULL_Model_Adj_ChrX")){
		
	
		re<-SKAT_With_NullModel_ADJ(Z, obj, kernel = kernel, method=method, weights.beta=weights.beta
		, weights = weights, impute.method = impute.method,  r.corr=r.corr, is_check_genotype=is_check_genotype
		, is_dosage = is_dosage, missing_cutoff=missing_cutoff, max_maf=max_maf, estimate_MAF=estimate_MAF, out.z=out.z)

	} else if(Check_Class(obj, "SKAT_NULL_Model_ChrX")){

		re<-SKAT_With_NullModel(Z,obj, kernel = kernel, method=method, weights.beta=weights.beta
		, weights = weights, impute.method = impute.method, r.corr=r.corr, is_check_genotype=is_check_genotype
		, is_dosage = is_dosage, missing_cutoff=missing_cutoff, max_maf=max_maf, estimate_MAF=estimate_MAF, out.z=out.z)

	} 
		
	class(re)<-"SKAT_OUT"
	return(re)

}



SKAT_ChrY<-function(Z, obj, kernel = "linear.weighted", method="davies", weights.beta=c(1,25)
                    , weights = NULL, impute.method = "fixed", r.corr=0, is_check_genotype=TRUE
                    , is_dosage = FALSE, missing_cutoff=0.15, max_maf=1, estimate_MAF=1, SetID=NULL){
  
  
  if(kernel != "linear" && kernel != "linear.weighted"){
    
    if(Check_Class(obj,  "SKAT_NULL_Model_Adj_ChrX")){
      msg<-sprintf("The small sample adjustment only can be applied for linear and linear.weighted kernel in the current version of SKAT! No adjustment is applied")
      warning(msg,call.=FALSE)
      obj<-obj$re1
    }
    
  }
  
  if(!Check_Class(obj, "SKAT_NULL_Model_Adj_ChrX") && !Check_Class(obj, "SKAT_NULL_Model_ChrX")){
    msg<-sprintf("The obj parameter should be a returned object from SKAT_Null_Model_ChrX!")
    stop(msg)
  }

  if(is.null(obj$obj.Y)){
    msg<-sprintf("SKAT_Null_Model_ChrX is not fitted with Model.Y==TRUE")
    stop(msg)
  }  
  
  obj.Y = obj$obj.Y
  id.male<-obj.Y$id.male
  SexVar = obj$SexVar[id.male]

  if(Check_Class(obj.Y, "SKAT_NULL_Model_ADJ")){
    id_include = obj.Y$re1$id_include
  } else if(Check_Class(obj.Y, "SKAT_NULL_Model")){
    id_include = obj.Y$id_include
  } else {
      msg<-sprintf("Class error in SKAT_ChrY %s", class(obj.Y))
      stop(msg)
  }
    


  if(length(id.male)==0){
    stop("No male subject!!") 
  }
  
  Z_male<-as.matrix(Z[id.male,])
  n.all<-length(id.male)
  
  out.z<-SKAT_MAIN_Check_Z(Z_male, n.all, id_include, SetID
                           , weights, weights.beta, impute.method, is_check_genotype
                           , is_dosage, missing_cutoff,max_maf=max_maf, estimate_MAF=estimate_MAF, Is.chrX=TRUE, SexVar=SexVar)
  

  if(Check_Class(obj.Y, "SKAT_NULL_Model_ADJ")){
    
    
    re<-SKAT_With_NullModel_ADJ(Z_male, obj.Y, kernel = kernel, method=method, weights.beta=weights.beta
                                , weights = weights, impute.method = impute.method,  r.corr=r.corr, is_check_genotype=is_check_genotype
                                , is_dosage = is_dosage, missing_cutoff=missing_cutoff, max_maf=max_maf, estimate_MAF=estimate_MAF, out.z=out.z)
    
  } else if(Check_Class(obj.Y, "SKAT_NULL_Model")){
    
    re<-SKAT_With_NullModel(Z_male,obj.Y, kernel = kernel, method=method, weights.beta=weights.beta
                            , weights = weights, impute.method = impute.method, r.corr=r.corr, is_check_genotype=is_check_genotype
                            , is_dosage = is_dosage, missing_cutoff=missing_cutoff, max_maf=max_maf, estimate_MAF=estimate_MAF, out.z=out.z)
    
  }  else {
    msg<-sprintf("Class error in SKAT_ChrY %s", class(obj.Y))
    stop(msg)
  }
  class(re)<-"SKAT_OUT"
  return(re)
  
}