4_IBD_Utilitty_Functions.R

#
# dCor_Profile() is the key function to calculate the dCor association profile across the region.
# It takes the following arguments as input:
# 
# 1. cc_sample: A dataframe with two columns: SeqID and Status
# First column, SeqID, shows the labels of the sequences.
# Second column, Status, represents the case/control status of each sequence. 0 = control, 1 = case
# 
dCor_Profile           = function(cc_sample, distance_matrix_list){
  
  
  # Calculate the phenotypic distance matrix by calling pheno_dist_matrix_calc() function.
  phenotye_dist_matrix  = pheno_dist_matrix_calc(cc_sample, prob = 0.05)
  
  
  # Calculate the dCor profile along the genome.
  n         = length(distance_matrix_list)
  dCor_Stat = numeric( length = n )
  for(i in 1:n){
    partition_dist  = distance_matrix_list[[i]][rownames(phenotye_dist_matrix),colnames(phenotye_dist_matrix)]
    dCor_Stat[i]    = perfectphyloR::dCorTest(Dx = partition_dist, Dy = phenotye_dist_matrix, nperm = 0)$Stat
  }
  
  
  # Returns dCor_Stat, a numeric vector, recording the dCor value across the region.
  return(dCor_Stat)
  
}




#
# pheno_dist_matrix_calc() calculates the phenotypic distance matrix. It takes 2 inputs:
#
# 1. A numeric vector y of length n, representing the phenotype of each sequence.
#    In the case of dichomotous trait it is coded as 1 = Case and 0 = Control
# 
# 2. The probability of disease in the population as prob. This must be a number 
#    between 0 and 1.
# 
# It returns a nxn phenotypic distance matrix according to calculation explained in:
# Burkett et al. 2014.(Link: https://www.karger.com/Article/Abstract/363443 )
# 
pheno_dist_matrix_calc     = function(cc_sample, prob){
  
  seq_names = cc_sample[,1]
  y         = cc_sample[,2]
  
  d = matrix(NA, nrow = length(y), ncol = length(y) )
  for( i in 1:length(y) ){
    for( j in 1:length(y) ){
      s_ij = (y[i] - prob) * (y[j] - prob)
      d[i,j] = 1 - s_ij
    }
  }
  
  diag(d) = 0
  rownames(d) = seq_names
  colnames(d) = seq_names
  
  return(d)
}






#
# dCorN_permute() runs the permutation test for Naive dCor (dCorN).
#
# The arguments of this function are as follows:
#
# 1. nperm = The desired number of permutation.
#
# 2. sample_data = the sample object generated by Sample_from_population() function.
#
# 3. distance_matrix_list = A list of distance matrices.
#
dCorN_permute          = function(nperm, sample_data, distance_matrix_list){

  
  # dCorN_permutation is a matrix. Each rows of this matrix represents the dCor profile for each permutation.
  # We add one extra row to fill it in with dCorN observed statistic later.
  dCorN_permutation  = matrix(NA, nrow = nperm + 1, ncol = length(sample_data$Posn$SNV_Position) )
  
  
  # ccLabel is the original case/control labeling of the individuals.
  ccLabel = sample_data$Genos$ccStatus

  
  
  # Load the permutation indices
  load("permute_indx.RData")

  
  # Running the permutation in parallel
  library(foreach)
  library(doParallel)
  
  # To run parallel on compute canada:
  # (Ref: https://docs.computecanada.ca/wiki/R)
  # 
  # Create an array from the NODESLIST environnement variable
  nodeslist = unlist(strsplit(Sys.getenv("NODESLIST"), split=" "))
  
  
  # Create the cluster with the nodes name. One process per count of node name.
  # nodeslist = node1 node1 node2 node2, means we are starting 2 processes on node1, likewise on node2.
  cl = makeCluster(nodeslist, type = "PSOCK") 
  registerDoParallel(cl)
  
  
  # To run parallel on your own PC/Laptop:
  # cl = makeCluster(detectCores() - 1)
  # registerDoParallel(cl)

  
  res = foreach(i = 1:nperm, .export = c("dCor_Profile","pheno_dist_matrix_calc") ) %dopar%{
           SeqID         = c(sample_data$CaseHapID, sample_data$ControlHapID) 
           Status        = unlist( lapply( ccLabel[ permute_indx[i,] ], FUN = function(x){rep(x,2)} ) )
           pCCLabel      = data.frame(SeqID = SeqID, Status = Status)
           dCor_Profile(cc_sample = pCCLabel, distance_matrix_list)
  }
  
  
  stopCluster(cl)
  
  for( i in 1:length(res) ){
    dCorN_permutation[i,] = res[[i]]
  }

  colnames(dCorN_permutation) = sample_data$Posn$SNV_Names
  return( dCorN_permutation )

}




# Case_reclassify() function takes the sample_data as the input. 
# We can use this function to reclassify case sequences based
# on their true carrier status.
# The output of this function is a vector indicating the 
# case control status of sequences after reclassification.
Case_reclassify        = function(sample_data){

  cSNV           = sample_data$poly_cSNV
  case_sequences = sample_data$CaseHapID  
  
  case_seq_carrier_status = colSums( sample_data$Haps$sample_haps[cSNV,case_sequences] )
  case_seq_carrier_status[case_seq_carrier_status >= 1] = 1
  names(case_seq_carrier_status) = NULL
  
  result = c(case_seq_carrier_status, rep(0, length(sample_data$ControlHapID)) )
  names(result) = c(sample_data$CaseHapID ,sample_data$ControlHapID)
  
  return(CarrierStatus = result)
  
}






#
# Mantel_Profile() is the key function to calculate the Mantel association profile across the region.
# It takes the following arguments as input:
# 
# 1. cc_sample: A dataframe with two columns: SeqID and Status
# First column, SeqID, shows the labels of the sequences.
# Second column, Status, represents the case/control status of each sequence. 0 = control, 1 = case
# 
# 2. distance_matrix_list: a list of distance matrices calculated based on partitions
#
Mantel_Profile           = function(cc_sample, distance_matrix_list){
  
  
  # Calculate the phenotypic distance matrix by calling pheno_dist_matrix_calc() function.
  phenotye_dist_matrix  = pheno_dist_matrix_calc(cc_sample, prob = 0.05)
  
  
  # Calculate the dCor profile along the genome.
  n           = length(distance_matrix_list)
  Mantel_Stat = numeric( length = n )
  for(i in 1:n){
    partition_dist  = distance_matrix_list[[i]][rownames(phenotye_dist_matrix),colnames(phenotye_dist_matrix)]
    Mantel_Stat[i]  = perfectphyloR::MantelTest(Dx = partition_dist, Dy = phenotye_dist_matrix, nperm = 0)$Stat
  }
  
  
  # Returns Mantel_Stat, a numeric vector, recording the Mantel value across the region.
  return(Mantel_Stat)
  
}



#
# Mantel_permute() runs the permutation test for Mantel.
#
# The arguments of this function are as follows:
#
# 1. nperm = The desired number of permutation.
#
# 2. sample_data = the sample object generated by Sample_from_population() function.
#
# 3. distance_matrix_list = A list of distance matrices.
#
Mantel_permute          = function(nperm, sample_data, distance_matrix_list){
  
  
  # Mantel_permutation is a matrix. Each rows of this matrix represents the Mantel profile for each permutation.
  # We add one extra row to fill it in with Mantel observed statistic later.
  Mantel_permutation  = matrix(NA, nrow = nperm + 1, ncol = length(sample_data$Posn$SNV_Position) )
  
  
  # ccLabel is the original case/control labeling of the individuals.
  ccLabel = sample_data$Genos$ccStatus
  
  
  
  # Load the permutation indices
  load("permute_indx.RData")
  
  
  # Running the permutation in parallel
  library(foreach)
  library(doParallel)
  
  # To run parallel on compute canada:
  # (Ref: https://docs.computecanada.ca/wiki/R)
  # 
  # Create an array from the NODESLIST environnement variable
  nodeslist = unlist(strsplit(Sys.getenv("NODESLIST"), split=" "))
  
  
  # Create the cluster with the nodes name. One process per count of node name.
  # nodeslist = node1 node1 node2 node2, means we are starting 2 processes on node1, likewise on node2.
  cl = makeCluster(nodeslist, type = "PSOCK") 
  registerDoParallel(cl)
  
  
  # To run parallel on your own PC/Laptop:
  # cl = makeCluster(detectCores() - 1)
  # registerDoParallel(cl)
  
  
  res = foreach(i = 1:nperm, .export = c("Mantel_Profile","pheno_dist_matrix_calc") ) %dopar%{
    SeqID         = c(sample_data$CaseHapID, sample_data$ControlHapID) 
    Status        = unlist( lapply( ccLabel[ permute_indx[i,] ], FUN = function(x){rep(x,2)} ) )
    pCCLabel      = data.frame(SeqID = SeqID, Status = Status)
    Mantel_Profile(cc_sample = pCCLabel, distance_matrix_list)
  }
  
  
  stopCluster(cl)
  
  for( i in 1:length(res) ){
    Mantel_permutation[i,] = res[[i]]
  }
  
  colnames(Mantel_permutation) = sample_data$Posn$SNV_Names
  return( Mantel_permutation )
  
}