SimCopy.R

#
# Copyright (C) 2013 EMBL - European Bioinformatics Institute
#
# This program is free software: you can redistribute it
# and/or modify it under the terms of the GNU General
# Public License as published by the Free Software
# Foundation, either version 3 of the License, or (at your
# option) any later version.
#
# This program is distributed in the hope that it will be
# useful, but WITHOUT ANY WARRANTY; without even the
# implied warranty of MERCHANTABILITY or FITNESS FOR A
# PARTICULAR PURPOSE. See the GNU General Public License
# for more details.
#
# Neither the institution name nor the name simcopy
# can be used to endorse or promote products derived from
# this software without prior written permission. For
# written permission, please contact <sbotond@ebi.ac.uk>.

# Products derived from this software may not be called
# simcopy nor may simcopy appear in their
# names without prior written permission of the developers.
# You should have received a copy of the GNU General Public
# License along with this program. If not, see
# <http://www.gnu.org/licenses/>.

##########################################################################/** 
#
# @RdocClass SimCopy
# \alias{simcopy}
#
# @title "The SimCopy class"
# 
# \description{ 
#
#   The \code{SimCopy} class simulates the evolution of copy number profiles along an input tree. It relies on the \pkg{\link{PhyloSim}}
#   package for performing the simulations by encoding the genomic regions as sites in sequences and using modified processes acting on them. 
#   Please note, that the \code{SimCopy} simulations are restricted to a single chromosome.
#   
#   The genomes are encoded as a sequence of sites containing integers identifying genomic regions. Negative integers represent
#   inverted genomic regions.
#   
#   The following processes are supported:
#
#   \itemize{
#   \item \code{deletion} - deletes genomic regions.
#   \item \code{duplication} - duplicates genomic regions.
#   \item \code{inversion} - changes the orientation of the genomic regions by taking the opposite of the corresponding integer.
#   \item \code{inverted duplication} - duplicates genomic regions and flips their orientation.
#   \item \code{translocation} - translocates a stretch of genomic regions.
#   }
#
#   The processes can be "activated" by specifying their parameters as arguments to the \code{SimCopy} constructor.
#   The processes are parametrised by a list with the following elements:
#
#   \itemize{
#   \item \code{rate} - the rate whereby the process proposes events.
#   \item \code{mean} - the mean of the truncated Geometric+1 distribution determining the number of genomic regions affected by individual 
#   events generated by the process.
#   \item \code{max} - the maximum number of genomic regions affected by a single event generated by the process. The default is \code{mean * 10}.
#   }
#
#   The simulations are run by the \link{Simulate.SimCopy}, which returns a list with the following elements:
#
#   \itemize{
#   \item \code{cnh} - a data frame containing the simulated copy number profiles.
#   \item \code{aln} - the simulated genomic region alignment as reported by \link{Simulate.PhyloSim}. 
#   Note that any element containing only dashes represents a single gap.
#   \item \code{fasta} - a string containing the simulated alignment in a fasta-like format.
#   Note that any element containing only dashes represents a single gap.
#   \item \code{phylo} - the \code{phylo} object used for simulation. It might be in different order than the object specified as argument.
#   \item \code{phylosim} - the \pkg{\link{PhyloSim}} object used for simulation.
#   \item \code{processes} - the \pkg{\link{PhyloSim}} processes used during the simulation.
#   }
#
#   \code{SimCopy} objects can be recycled for simulations along different trees. See the \link{Simulate.SimCopy} method for an example.
#
#	@classhierarchy
# }
#	
# @synopsis
#	
# \arguments{
# 	\item{root.size}{The number of regions in the root genome (100 by default).}
#   \item{deletion}{A list containing the parameters of the deletion process.}
#   \item{duplication}{A list containing the parameters of the duplication process.}
#   \item{inv.duplication}{A list containing the parameters of the inverted duplication process.}
#   \item{inversion}{A list containing the parameters of the inversion process.}
#   \item{translocation}{A list containing the parameters of the translocation process.}
# 	\item{...}{Additional arguments.}
#	}
# 
# \section{Fields and Methods}{ 
# 	@allmethods
# }
# 
# \examples{ 
#   # The following tiny examples illustrate the
#   # effects of individual processes:    
#
#   tree<-rcoal(2)  # We will use this tiny tree in the examples below.
#   rate<-0.08      # Common rate for the small examples.
#
#   # Simulating deletions and dealing with the results:
#   cat("\nSimulating deletions:\n")
#   # Construct a SimCopy object:
#   sc <- SimCopy(
#       root.size=40,
#       deletion=list(rate=rate, mean=2)
#    )
#   # Run simulation:
#   res<-Simulate(sc, tree)
#   # Deal with the simulation results:
#   print(res$aln)  # print out the simulated alignment
#   print(res$cnh)  # print out the simulated copy number history
#   print(res$fasta)# print out the fasta alignment
#   summary(res$phylosim) # get the details of the PhyloSim object used for simulations
#   summary(res$processes[[1]]) # get the details of the deletion process
#   plot(res$processes[[1]])    # plot the distribution of deletion lengths
#
#   # Simulate duplications and print out the resulting alignment:
#   cat("\nSimulating duplications:\n")
#   # Construct a SimCopy object:
#   sc <- SimCopy(
#       root.size=20,
#       duplication=list(rate=rate, mean=2)
#    )
#   print( Simulate(sc, tree)$aln )
#
#   # Simulate inverted duplications and print out the resulting alignment:
#   cat("\nSimulating inverted duplications:\n")
#   # Construct a SimCopy object:
#   sc <- SimCopy(
#       root.size=20,
#       inv.duplication=list(rate=rate, mean=2)
#    )
#   print( Simulate(sc, tree)$aln )
#
#   # Simulate inversions and print out the resulting alignment:
#   cat("\nSimulating inversions:\n")
#   # Construct a SimCopy object:
#   sc <- SimCopy(
#       root.size=20,
#       inversion=list(rate=rate, mean=2)
#    )
#   print( Simulate(sc, tree)$aln )
#
#   # Simulate translocations and print out the resulting alignment:
#   cat("\nSimulating translocations:\n")
#   # Construct a SimCopy object:
#   sc <- SimCopy(
#       root.size=20,
#       translocation=list(rate=rate, mean=2)
#    )
#   print( Simulate(sc, tree)$aln )
#
#   ##
#   ## In the following simulation we will use all the processes above 
#   ## and we will attempt to recover the topology using simple hierarchical
#   ## clustering of the copy number profiles.
#   ##
#
#   tree<-rcoal(6)
#   rate<-0.05
#   sc <- SimCopy(
#       root.size=50,
#       deletion=list(rate=rate, mean=2),
#       duplication=list(rate=rate, mean=2),
#       inv.duplication=list(rate=rate, mean=2),
#       inversion=list(rate=rate, mean=2),
#       translocation=list(rate=rate, mean=2)
#    )
#   res<-Simulate(sc, tree, anc=FALSE) # discard internal nodes
#   # Print out the simulate genomic region alignment through
#   # the underlying PhyloSim object:
#   plot(res$phylosim)
#   # Calculate distances between copy number profiles:
#   d<-dist(res$cnh)
#   # Cluster the copy number profiles:
#   hc<-hclust(d)
#
#   # Relabel the tips of the true tree and plot it out:
#   tree$tip.label<-1:length(tree$tip.label)
#   plot(tree)
#   # Plot out the results of hierarchical clustering:
#   plot(hc)
#
# }
# 
# @author
#
# \seealso{ 
# 	See the \link{Simulate.SimCopy} method for the details of running simulations.
#   See also the \pkg{\link{phylosim}} and \pkg{\link{ape}} packages.
# }
# 
#*/###########################################################################
setConstructorS3(
  "SimCopy",
  function(
    root.size=100,
    deletion=list(rate=NA, mean=NA, max=NA),
    duplication=list(rate=NA, mean=NA, max=NA),
    inv.duplication=list(rate=NA, mean=NA, max=NA),
    inversion=list(rate=NA, mean=NA, max=NA),
    translocation=list(rate=NA, mean=NA, max=NA),
    ... 
    ){

    this <- extend(SCRoot(), "SimCopy",
        root.size=root.size,
        deletion=deletion,
        duplication=duplication,
        inv.duplication=inv.duplication,
        inversion=inversion,
        translocation=translocation
    );

    this$.processes<-.construct.processes(this)
    return(this)
  },
  enforceRCC=TRUE
);

# Function constructing processes:
.construct.processes<-function(this) {

        # Construct processes:
        processes<-list()
        if(!is.na(this$deletion$rate)) {
            del <-.construct.deletor(rate=this$deletion$rate, lenMean=this$deletion$mean, lenMax=this$deletion$max)
            processes<-c(processes, list(del))
        }

        if(!is.na(this$duplication$rate)) {
            dup <-.construct.duplicator(rate=this$duplication$rate, lenMean=this$duplication$mean, lenMax=this$duplication$max)
            processes<-c(processes, list(dup))
        }

        if(!is.na(this$inv.duplication$rate)) {
            inv.dup <-.construct.inv.duplicator(rate=this$inv.duplication$rate, lenMean=this$inv.duplication$mean, lenMax=this$inv.duplication$max)
            processes<-c(processes, list(inv.dup))
        }

        if(!is.na(this$inversion$rate)) {
            inv <-.construct.invertor(rate=this$inversion$rate, lenMean=this$inversion$mean, lenMax=this$inversion$max)
            processes<-c(processes, list(inv))
        }

        if(!is.na(this$translocation$rate)) {
            trl <-.construct.translocator(rate=this$translocation$rate, lenMean=this$translocation$mean, lenMax=this$translocation$max)
            processes<-c(processes, list(trl))
        }

        return(processes)
}

# Set up a truncated geometric+1 distribution according to the
# specified mean and maximum:
.setup.ldist<-function(l.mean, l.max=NULL) {
    if(is.na(l.mean) || l.mean < 1) {
        stop("Illegal mean length: ", l.mean)
    }
    if(is.null(l.max)){
        l.max <- l.mean * 10
    } 
    sizes   <- 0:(l.max-1)
    probs   <- dgeom(sizes, prob=1/l.mean)
    probs   <- probs/sum(probs)
    sizes   <- sizes + 1
    return(
        list(sizes=sizes, probs=probs)
    )
}

###########################################################################/**
#
# @RdocMethod Simulate
# 
# @title "Method for simulating copy number histories" 
# 
# \description{ 
#	@get "title".
#
#   This method takes a \link{SimCopy} object and a \code{phylo} object generated by the \pkg{\link{ape}}
#   package and simulates the evolution of genomic regions using the \pkg{\link{PhyloSim}} package.
# } 
# 
# @synopsis 
# 
# \arguments{ 
# 	\item{this}{A \code{SimCopy} object.} 
#   \item{phylo}{A phylo object constructed by the \pkg{\link{ape}} package.}
#   \item{anc}{Save copy number profiles corresponding to internal nodes (TRUE by default).}
#   \item{quiet}{Suppress \pkg{\link{PhyloSim}} verbose output (FALSE by default).}
# 	\item{...}{Not used.} 
# } 
# 
# \value{ 
#   The return value is a list with the following elements:
#   \itemize{
#   \item \code{cnh} - a data frame containing the simulated copy number profiles.
#   \item \code{aln} - the simulated genomic region alignment as reported by \link{Simulate.PhyloSim}.
#   Note that any element containing only dashes represents a single gap.
#   \item \code{fasta} - a string containing the simulated alignment in a fasta-like format.
#   Note that any element containing only dashes represents a single gap.
#   \item \code{phylo} - the \code{phylo} object used for simulation. It might be in different order than the object specified as argument.
#   \item \code{phylosim} - the \pkg{\link{PhyloSim}} object used for simulation.
#   \item \code{processes} - the \pkg{\link{PhyloSim}} processes used during the simulation.
#   }
# } 
# 
# \examples{
#
#   # The following example illustrates how to simulate copy number 
#   # evolution along many trees by the processes defined in a single
#   # SimCopy object. See the class documentation for more detailed examples.
#
#   # Construct a SimCopy object with a duplication  process:
#   sc <- SimCopy(
#       root.size=5,
#       duplication=list(rate=0.05, mean=1)
#    )
#
#   for(i in 1:4) {
#       cat("\nReplicate: ", i, "\n\n")
#       # Generate a coalescent tree:
#       tree<-rcoal(4)
#       # Run simulation and print out the copy number history:
#       print( Simulate(sc, tree, anc=FALSE, quiet=TRUE)$cnh )
#   }
#
# } 
# 
# @author 
# 
# \seealso{ 
# 	@seeclass 
# } 
# 
#*/###########################################################################
setMethodS3(
  "Simulate",
  class="SimCopy",
  function(
    this,
    phylo,
    anc=TRUE,
    quiet=FALSE,
    ...
  ){

        if(missing(phylo)){
            stop("Tree must be given as an APE phylo object!")
        }

        # Construct the root genome object:
        tmp <- .construct.root(this$root.size)
        root.seq <- tmp$seq
        this$.alphabet <- tmp$alphabet

        # Attach the processes to the root genome:
        setProcesses(root.seq, list(this$.processes))
       
        # Construct the Phylosim simulation object: 
        psim<-PhyloSim(root.seq=root.seq, phylo=phylo)
        # Simulate copy number evolution:
        Simulate(psim, quiet=quiet)
        # Construct copy number history from the PhyloSim alignment:
        cnh <- .getCnh(this, psim, anc)
        phylo <-psim$phylo
        fasta <-.buildFasta(this, psim, anc)
        return(
            list(
                phylo= phylo,
                aln  = psim$alignment,
                cnh  = cnh,
                fasta=fasta,
                phylosim=psim,
                processes=this$.processes
            ) 
        )
  },
  private=FALSE,
  protected=FALSE,
  overwrite=TRUE,
  conflict="warning"
);

# Private method constructing fasta from the alignment matrix:
setMethodS3(
  ".buildFasta",
  class="SimCopy",
  function(
            object,
            this,
    		anc,
    		...
  ){

        fasta <- ""
		if(any(is.na(this$.alignment))){
			warning("Alignment is undefined, nothing to save!\n");
			return();
		}
		else {
			if(anc){
				for(i in 1:dim(this$.alignment)[[1]]){
					fasta <- paste(fasta, ">",rownames(this$.alignment)[[i]],"\n",sep="");
					fasta <- paste(fasta, paste(this$.alignment[i,],collapse="\t"),"\n", sep="");
				}
			} else {
				for(i in 1:dim(this$.alignment)[[1]]){
					name<-rownames(this$.alignment)[[i]];
					if(!any((length(grep("^Node \\d+$",name,perl=TRUE,value=FALSE)) > 0),(length(grep("^Root node \\d+$",name,perl=TRUE,value=FALSE)) > 0))){
						fasta<-paste(fasta, ">",name,"\n", sep="");
						fasta<-paste(fasta, paste(this$.alignment[i,],collapse="\t"),"\n", sep="");
					}

				}
			}
		}

		return(fasta);
  },
  private=TRUE,
  protected=FALSE,
  overwrite=FALSE,
  conflict="warning",
  validators=getOption("R.methodsS3:validators:setMethodS3")
);

# Private method for constructing the copy number history
# from the SimCopy alignment.
setMethodS3(
  ".getCnh",
  class="SimCopy",
  function(
    this,
    sim,
    anc=FALSE,
    ...
  ){
        sym     <- 1:this$root.size
        nodes   <- getNodes(sim)
        node.names   <- c()
        cnh     <- data.frame()
        for (n in nodes) {
            if(!anc && !is.tip(sim, n)){
                next
            }
            s<- abs(as.numeric(getSeqFromNode(sim, n)$states))
            cnh<-rbind(cnh, tabulate(s, this$root.size))
            node.names<-c(node.names, n)
        }
        row.names(cnh) <- node.names
        colnames(cnh) <- 1:this$root.size
        return(cnh)
  },
  private=TRUE,
  protected=FALSE,
  overwrite=TRUE,
  conflict="warning"
);

## Method: summary.SimCopy
##
###########################################################################/**
#
# @RdocMethod summary
#
# @title "Summarize the properties of an object"
#
# \description{
#       @get "title".
# }
#
# @synopsis
#
# \arguments{
#       \item{object}{A SimCopy object}
#       \item{...}{Not used.}
# }
#
# \value{
#  Returns a SCRoot.Summary object.
# }
#
# \examples{
#	# Create a SimCopy object with a deletion process defined:
#	sim<-SimCopy(
#      root.size=10,
#      deletion=list(rate=0.5, mean=20, max=20)
#	);
#       # get a summary
#       summary(sim)
# }
#
# @author
#
# \seealso{
#       @seeclass
# }
#
#*/###########################################################################
setMethodS3(
  "summary",
  class="SimCopy",
  function(
    object,
    ...
  ){
     this<-object;
     this$.summary$"Root genome length"             <-this$root.size;
     this$.summary$"Deletion process"               <-.build.process.summary(this$deletion);
     this$.summary$"Duplication process"            <-.build.process.summary(this$duplication);
     this$.summary$"Inverted duplication process"   <-.build.process.summary(this$inv.duplication);
     this$.summary$"Inversion process"              <-.build.process.summary(this$inversion);
     this$.summary$"Translocation process"          <-.build.process.summary(this$translocation);
     NextMethod();

  },
  private=FALSE,
  protected=FALSE,
  overwrite=FALSE,
  conflict="warning",
  validators=getOption("R.methodsS3:validators:setMethodS3")
);

.build.process.summary<-function(l){
    if(is.na(l$rate)){
        return(NA)
    }
    if(l$rate == 0){
        return(NA)
    }
    t<-paste(
        "\n\tRate:\t\t", l$rate, "\n", 
        "\tMean length:\t", l$mean, "\n", 
        "\tMaximum length:\t", l$max,  
        sep="")
        return(t)
}