R/Qassignate.R

In gammadnamix: Package for intepretation of quantitative DNA mixtures

#' @title Qassignate
#' @author Oyvind Bleka <Oyvind.Bleka.at.fhi.no>
#' @description Q-assignation. It also takes care of include new alleles into popFreq by assigning it lowest observed frequnce.
#' @details Assignes non-shown alleles as one single allele.
#' @param samples A List with samples which for each samples has locus-list elements with list elements adata and hdata. 'adata' is a qualitative (allele) data vector and 'hdata' is a quantitative (peak heights) data vector.
#' @param popFreq A list of allele frequencies for a given population.
#' @param refData Reference objects with list element [[s]]$adata[[i]]. The list element has reference-list with list-element 's' having a loci-list adata with list-element 'i storing qualitative data.
#' @return ret A list(popFreq,refData) with Q-assignated alleles. 
#' @export

Qassignate <- function(samples,popFreq,refData=NULL,doQ=TRUE) {
 popFreq2 <- popFreq
 refData2 <- refData
 locs <- names(popFreq)
 minF <- min(unlist(popFreq)) #lowest observed frequency
 for(loc in locs) { #make Q-assignation for each loci
  evid <- evid0 <- unique(unlist( lapply(samples,function(x) x[[loc]]$adata) )) #vectorize alleles for all replicates
  if(!is.null(refData)) evid  <- unique(c(evid,unlist(refData[[loc]])))

  #if new alleles in evidence or references not in popFreq
  newa <- evid[!evid%in%names(popFreq[[loc]])]   
  if(length(newa)>0) {
   tmp <- names(popFreq[[loc]])
   popFreq[[loc]] <- c(popFreq[[loc]],rep(minF,length(newa)))
   names(popFreq[[loc]]) <-  c(tmp,newa)
   print(paste0("Locus ",loc,": Allele(s) ",paste0(newa,collapse=",")," was inserted with frequency ",minF))
  }
  popFreq[[loc]] <- popFreq[[loc]]/sum(popFreq[[loc]]) #normalize

  if(doQ) {#if Q-assignate
   tmp <- popFreq[[loc]][names(popFreq[[loc]])%in%evid0] #find observed alleles
   if(length(tmp)<length(popFreq[[loc]])) { 
    tmp <- c(tmp,1-sum(tmp))
    names(tmp)[length(tmp)] <- "99"
   }
   popFreq2[[loc]] <- tmp
   if(!is.null(refData)) { #insert 99 as default allele of missing refs
    newP <- names(popFreq2[[loc]]) 
    if(!all(unlist(refData[[loc]])%in%newP)) { #there was some missing alleles
     for(k in 1:length(refData[[loc]])) {
      refData2[[loc]][[k]][!refData[[loc]][[k]]%in%newP] <- "99" #insert missing
     }
    }
   }
  } else { #if don't Q-assignate
   popFreq2 <- popFreq
   refData2 <- refData
  }
 }
 return(list(popFreq=popFreq2,refData=refData2))
}