R/LD.sign.test.R

In pould: Phased or Unphased Linkage Disequilibrium

## Performing sign tests on Phased and UnPhased LD values -- Steven J. Mack April 6, 2020 v1.0

#' Perform the sign test on paired LD values for phased and unphased haplotypes 
#' 
#' Applies binom.test() to pairs of {D'}, {Wn}, {WLoc1/Loc2}, {WLoc2/Loc1} and number of haplotypes values in the *_LD_results.csv files generated by LDWrap(). 
#' 
#' This function returns p-values for the sign test comparing the phased and unphased values of each LD measure, as well as the number of haplotypes, for each locus pair in a dataset tested using LDWrap(). It also returns the number of locus pairs for which the value in quesiton is higher in unphased haplotypes than phased haplotypes, the number of locus pairs in which the values are the same, and the total number of locus pairs asssessed. This function writes a results file to the working directory for each dataset, and will optionally display those results in the console. 
#' 
#' @param dataName The "base" name of the _LD_result.csv files generated by LDWrap() without the "_Phased_LD_results.csv" or "_Unphased_LD_results.csv" suffixes. See Examples, below. If the corresponding "<dataName>_Phased_LD_results.csv" or "<dataName>_Unphased_LD_results.csv" files are not found, the funciton will halt with a notification. 
#' @param verbose A boolean identifying if messages about function progress and results should be displayed in the console (verbose=TRUE) or not (verbose=FALSE). The default is verbose=TRUE.
#' @param returnFrame A boolean identifying if a data frame of results should be returned (returnFrame=TRUE). If 'returnFrame=FALSE', a CSV file of results named <dataName>_LD-sign-test_results.csv is written in the directory specified by the 'resultDir' parameter. The default is returnFrame=TRUE.
#' @param resultDir The directory into which the CSV file of results should be written when 'returnFrame=FALSE'. The default is the directory specified by tempdir().
#' @keywords sign test linkage disequilibrium
#' @importFrom utils write.table read.table read.csv
#' @importFrom stats binom.test
#' @importFrom graphics plot
#' @return A data frame of five columns ({D'}, {Wn}, {WLoc1/Loc2}, {WLoc2/Loc1} and N_Haplotypes) and four rows (#unphased > phased, #unphased = phased, #locus pairs and p-values).
#' @note When verbose=TRUE, LD.sign.test() writes a table of results to the console with the column headers "Measure", "#U > P", "#U = P" and "p-value". Column "#U > P" identifies the number locus pairs for which that measure for unphased haplotypes (U) was greater than that measure for phased haplotypes (P). Similarly, column "#U = P" identiifes the number of locus pairs for each meausure where the value of that measure was the same in phased and unphased haplotypes.
#' @note Only the significance of the sign test is reported; when a significant trend is indicated, the directionality of the trend is not reported.
#' @export
#' @examples 
#' # Using LDWrap() to analyze the first 10 rows of the drb1.dqb1.demo dataset.
#' # LDWrap() results are saved in the temporary directory as 
#' # "hla-family-data_Phased_LD_results.csv" and 
#' # "hla-family-data_Unphased_LD_results.csv", respectively. 
#' LDWrap(drb1.dqb1.demo[1:10,])
#' LDWrap(drb1.dqb1.demo[1:10,],phased=FALSE)
#' exampleData <- paste(tempdir(),"hla-family-data",sep=.Platform$file.sep)
#' LD.res <- LD.sign.test(exampleData)
#' # Return only a data frame of results.
#' LD.res <- LD.sign.test(exampleData,verbose=FALSE)
#' @references Osoegawa et al. Hum Immunol. 2019;80(9):633 (https://doi.org/10.1016/j.humimm.2019.01.010)

LD.sign.test <- function(dataName,verbose=TRUE,returnFrame=TRUE,resultDir=tempdir()) {
  
if(file.exists(paste(dataName,"_Phased_LD_results.csv",sep="")) && file.exists(paste(dataName,"_Unphased_LD_results.csv",sep=""))) {

phaseDataRaw <- read.csv(file=paste(dataName,"_Phased_LD_results.csv",sep=""),header = TRUE)
unphaseDataRaw <- read.csv(file=paste(dataName,"_Unphased_LD_results.csv",sep=""),header=TRUE)

## need to ignore any locus pairs for which no calculations were performed
phaseData <- phaseDataRaw[!is.na(phaseDataRaw[,6]),]
unphaseData <- unphaseDataRaw[!is.na(unphaseDataRaw[,6]),]

if(nrow(phaseData) < nrow(phaseDataRaw)) { # Need to convert factors to numeric values when all 55 locus pairs were not calculated
  phaseData[,1:6] <- lapply(phaseData[,1:6],as.character)
  unphaseData[,1:6] <- lapply(unphaseData[,1:6],as.character)
  if(verbose){cat(paste("LD.sign.test: Removing ",nrow(phaseDataRaw) - nrow(phaseData), " uncalculated Locus Pairs.",'\n',sep=""))}
  phaseData[,2:6] <- lapply(phaseData[,2:6],as.numeric)
  unphaseData[,2:6] <- lapply(unphaseData[,2:6],as.numeric)
  
    if (length(phaseData[,1]) != length(unphaseData[,1]) || sum((phaseData[,1] != unphaseData[,1]) == TRUE) > 0) { # in this case, there are some locus pairs with data in one dataset, but not the other
      
    locMatch <- intersect(phaseData[,1],unphaseData[,1])
    allPairs <- union(phaseData[,1],unphaseData[,1])
    toIgnore <- allPairs[!allPairs %in% locMatch]
    
    phaseData <- phaseData[!phaseData[,1] %in% toIgnore,]
    unphaseData <- unphaseData[!unphaseData[,1] %in% toIgnore,]

    if(verbose){cat(paste("Mismatch between phased and unphased locus pairs. Excluding ",paste(toIgnore,collapse = " "),".","\n",sep=""))}
    
  }
}
diffs <- unphaseData[,2:6] - phaseData[,2:6]
nPos <- colSums(diffs > 0)
nZero <- colSums(diffs == 0)
nTrials <- nrow(phaseData)
pvals <- rep(NA,5)

for(n in 1:5) {pvals[n] <- binom.test(nPos[n],nTrials)$p.value}

result <- rbind(nPos,nZero,nTrials,pvals)
colnames(result) <- c("D'","Wn","WLoc1/Loc2","WLoc2/Loc1","N_Haplotypes")
  
rownames(result) <- c("#unphased > phased","#unphased = phased","#locus pairs","p-values")

if (!returnFrame) {write.table(result,paste(resultDir,paste(basename(dataName),"_LD-sign-test_results.csv",sep=""),sep=.Platform$file.sep),append = FALSE,quote = FALSE,sep = ",",row.names = TRUE,col.names = NA)}

if(verbose){
cat(paste("Sign Test results for the ",dataName," dataset for ",nTrials," locus pairs.",'\n',sep=""))
cat(paste("Measure     ","#U > P","#U = P","p-value",'\n',sep="\t"))
cat(paste("D'          ",nPos[1],nZero[1],signif(pvals[1],4),'\n',sep="\t"))  
cat(paste("Wn          ",nPos[2],nZero[2],signif(pvals[2],4),'\n',sep="\t"))  
cat(paste("WLoc1/Loc2  ",nPos[3],nZero[3],signif(pvals[3],4),'\n',sep="\t"))  
cat(paste("WLoc2/Loc1  ",nPos[4],nZero[4],signif(pvals[4],4),'\n',sep="\t"))  
cat(paste("# Haplotypes",nPos[5],nZero[5],signif(pvals[5],4),'\n',sep="\t"))  
}

if (returnFrame) {return(result)}

} else {
  
  if(!file.exists(paste(dataName,"_Phased_LD_results.csv",sep="")) && !file.exists(paste(dataName,"_Unphased_LD_results.csv",sep=""))) {
      warning(paste("LD.sign.test: Files ",dataName,"_Phased_LD_results.csv"," and ",dataName,"_Unphased_LD_results.csv"," are not found.",sep=""))
  } else {if(!file.exists(paste(dataName,"_Phased_LD_results.csv",sep=""))) { 
      warning(paste("LD.sign.test: File ",dataName,"_Phased_LD_results.csv"," is not found.",sep=""))
      } else {warning(paste("LD.sign.test: File ",dataName,"_Unphased_LD_results.csv"," is not found.",sep=""))}
    }
    
    }
}