#' Get infection times
#' @param record MCMC output produced by inferTTree
#' @param k Case whose posterior infection times are to be extracted. Either an integer or a string matching one of the case names in the data
#' @return A vector of posterior infection times for case k
#' @export
#' @examples
#' getInfectionTimes(record)
getInfectionTimes <- function(record,k) {
if (is.numeric(k)) {
mytimes= vapply(1:length(record), function(x) { tt=extractTTree(record[[x]]$ctree); return(tt$ttree[k,1])},FUN.VALUE=1);
return(mytimes)}
else {
mytimes= vapply(1:length(record), function(x) { tt=extractTTree(record[[x]]$ctree); ii=which(tt$nam==k); return(tt$ttree[ii,1])},FUN.VALUE=1);
return(mytimes)}
}
#' Get generation times
#' @param ctree Combined tree object 'ctree' from the record produced by inferTTree
#' @return Vector of times between becoming infected and infecting others (generation times) in the posterior
#' @examples
#' getGenerationTimes(record[[1]]$ctree)
getGenerationTimes <- function(ctree) { tt=extractTTree(ctree)$ttree;
# 3rd column of ttree lists the infectors; exclude source
infectors=tt[,3]; infectors=infectors[infectors!=0];
# times at which each infector infected others:
infothers=tt[tt[,3]!=0,1];
# times at which each infector was herself infected:
gotinfd=tt[infectors,1];
return(infothers-gotinfd);}
#' Get the number of unsampled cases
#' @param ctree Combined tree object 'ctree' from the record produced by inferTTree
#' @return The number of unsampled cases in ctree
#' @examples
#' getNumberUnsampled(ctree)
getNumberUnsampled <- function(ctree) {tt=extractTTree(ctree); return(sum(is.na(tt$ttree[,2])))}
# getTTreeDistInfo not currently used. If so, add burnin parameter.
#' Create list of wiw information in order to compute transmission tree distances
#' @param record MCMC output produced by inferTTree
#' @param skipnum Number of record entries to skip, ie skipnum=10 uses the 1st, 11th, 21st ... elements of the record
#' @return list of MRCI information required by wiwTreeDist, one entry for each transmission tree that is included
getTTreeDistInfo <- function(record,skipnum=1) {
ind=seq(from=1, by=skipnum,to=length(record))
record=record[ind]
matList <- lapply(1:length(record), function(x) {
info <- extractTTree(record[[x]]$ctree)$ttree
wiw <- cbind(info[,3],1:length(info[,1]))
findMRCIs(wiw)$mrciDepths
})
return(matList)
}
#' Get infection likelihoods for a given case
#' @param record MCMC output produced by inferTTree
#' @param case Index or label can be supplied
#' @param burnin Proportion of record entries to use
#' @export
getLikelyInfectors <- function(record, case, burnin=0.5){
# Who-infected-who matrix from TransPhylo
myMat <- computeMatWIW(record, burnin)
return(myMat[,case])
}
#' Get infection-to-sampling times
#' @param ctree TransPhylo tree
#' @return A vector of posterior times to sampling
#' @examples
#' getTimesToSampling(ctree)
getTimesToSampling <- function(ctree) { tt=extractTTree(ctree)$ttree;
ns=sum(!is.na(tt[,2]))
return(tt[1:ns,2]-tt[1:ns,1])
}
#' Get infection dates
#' @param ctree TransPhylo tree
#' @return A vector of posterior infection dates
#' @examples
#' getInfectionDates(ctree)
getInfectionDates <- function(ctree) { tt=extractTTree(ctree)$ttree;
ns=sum(!is.na(tt[,2]))
return(tt[1:ns,1])
}
#' add alignments to list
#' @param alist (possibly empty) list of alignments in DNAbin format
#' @param fileName fasta file name
#' @return list of alignments
#' @export
addAlignment <- function(alist, fileName){
alignment <- read.dna(system.file("extdata", fileName, package = "transnp", mustWork = TRUE), format="fasta")
alist[[length(alist)+1]] <- alignment
return (alist)
}
#' create list of dates to match alignments
#' @param alist (possibly empty) list of alignments in DNAbin format
#' @param allDates data frame containg all IDs and dates as read from a csv file
#' @return list of dates
#' @export
setUpDates <- function(alist, allDates){
dateList <- list()
for (align in alist){
theseDates <- numeric(length(labels(align)))
for (i in seq(length(labels(align)))){
theseDates[[i]] <- allDates[match(labels(align)[[i]], allDates[,1]), 2]
}
dateList[[length(dateList)+1]] <- theseDates
}
return(dateList)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.