smt-0.3.9-BETA: Single Molecule Tracking

## selComponentTracks-methods
##
##
###############################################################################
##' @name selComponentTracks
##' @aliases selComponentTracks
##' @title selComponentTracks
##' @rdname selComponentTracks-methods
##' @docType methods
##'
##' @description Select trajectory based on component fitting on diffusion
##'   coefficient.
##' @usage selComponentTracks(fit,likelihood=0.9,dcoef,log.transformed=F,output=F)
##' @param fit Component fitting result form fitNormDistr() function.
##' @param likelihood The likelihood of a trajecotry to be in fitted group. This parameter specifies the strigency of selecting trajectories to be in the fitted group and therefore influence the number of trajectories been selected.
##' @param dcoef Diffusion coefficent calcualted by Dcoef, which provide the link between trajecotry index and diffusion coefficent.
##' @param log.transformed A flag indicating if the fitting is been log transformed, select TRUE if fitting was done in fitNormDistr(log.transform = T,..). This parameter will be removed in later version by directly read the info from the output of fitNormalDistr() function.
##' @param output A logical indicating if output of selected trajectory index, which can be used for plot individual trajectory using plotTrack.


##' @return
##' \itemize{
##' \item{combined list of trackll} The result is a combined list of selected trajectories. The list is one level higher than trackll, use subsetting to output trackll.e.g. trackll[[1]], or trackll[["SWR1"]]. }

##' @examples
##'
##' ## selComponentTracks() usage
##' # 1. select componentTracks per folder (cross movie) by using compareFolders
##' # 2. select componentTracks per movie base, use plotComponentTracks to plot component tracks back to initial Nuclei image.
##'
##' ## 1. select componentTracks per folder (cross movie) by using compareFolders
##' folder1=system.file("extdata","SWR1",package="smt")
##' folder2=system.file("extdata","HTZ1",package="smt")
##' trackll=compareFolder(c(folder1,folder2))
##' MSD=msd(trackll=trackll)
##' dcoef=Dcoef(MSD,dt=6,plot=T,output=F)
##' # fit dcoef
##' # for replication purpose set seed to fix number
##' fit=fitNormDistr(dcoef,components=2,log.transform=T,combine.plot=F,output=F,seed=484)
##'
##' # select component tracks from fitting
##' trackll.sel=selComponentTracks(fit=fit,likelihood = 0.9,dcoef = dcoef,log.transformed = T,output = F)
##' # subset component tracks to further analyze msd, dcoef
##' trackll.swr1=trackll.sel[["SWR1"]]
##' msd(trackll.swr1,plot=T)
##' msd(trackll.swr1,summarize=T,plot=T)
##' Dcoef(trackll=trackll.swr1,plot=T)
##' plotTrackOverlay(trackll.swr1)
##'
##' # plotNucTrackOverlay(folder=folder1, trackll=trackll.swr1)
##' dwellTime(trackll.swr1)
##'
##' # Output trajectory index to plot individually
##' # trackll.sel=selComponentTracks(fit=fit,likelihood = 0.9,dcoef = dcoef,log.transformed = T,output = T)
##' # specify index file path.
##' index.file=system.file("extdata","INDEX","componentTrackID-SWR1.comp.1.csv",package="smt")
##' index.file2=system.file("extdata","INDEX","componentTrackID-SWR1.comp.2.csv",package="smt")
##' movie.folder=system.file("extdata","SWR1_2",package="smt")
##' plotTrackFromIndex(index.file=index.file,movie.folder = movie.folder)
##' plotTrackFromIndex(index.file=index.file2,movie.folder = movie.folder)
##'
##'
##' ## 2. select componentTracks per movie base, use plotComponentTracks to plot component tracks back to initial Nuclei image.
##' ## plotComponentTrackOverlay
##' folder3=system.file("extdata","SWR1_2",package="smt")
##' trackll=readDiatrack(folder3,merge=F)
##'
##' ## use merge=T for per folder comparison, the analsyis result can't be plot back to original image
##' ## To see component tracks on original nuclei image, set merge=F (for per movie analysis)
##' ## may not make much sense to msd on individual movie, however for plot component track back to original nuclei image.
##'
##' ## compute MSD
##' MSD=msd(trackll=trackll,plot=T)
##' msd(trackll=trackll,summarize=T,plot=T)
##'
##' ## calculate Dcoef
##' dcoef=Dcoef(MSD=MSD,method="static",plot=TRUE)
##'
##' ## fit normal distribution to define component
##' ## set seed to reproduce results (see fitNormalDistr() for details on seed)
##' fit=fitNormDistr(dcoef,components=2,log.transform=T,combine.plot=F,output=F,seed=481)
##'
##' ## select component tracks based on fitting
##' trackll.sel=selComponentTracks(fit=fit,likelihood = 0.9,dcoef = dcoef,log.transformed = T,output = F)
##' ## plot component tracks
##' plotComponentTrackOverlay(folder=folder3,trackll.sel=trackll.sel)
##'
##'



##' @export selComponentTracks
##'

###############################################################################


selComponentTracks=function(
    fit,likelihood=0.9,dcoef,log.transformed=F,output=F){

    comp=list()
    length(comp)=length(fit)
    names(comp)=names(fit)


    comp.track.index.lst=list()
    length(comp.track.index.lst)=length(fit)
    names(comp.track.index.lst)=names(fit)

    # secondary /internal list structure of comp.track.index.lst is created in a
    # loop below

    # lapply(comp.track.index.lst,function(){})

    #---------------------------------------------------------------------------
    # select tracks from fitNormDistr() outputs
    for (i in 1:length(fit)){

        comp[[i]]=data.frame(fit[[i]]$posterior)

        # specify secondary list structure (length and name)
        # length(comp.track.index.lst[i])=dim(comp[[i]])[2]

        # creates internal list structure of comp.track.index.lst
        comp.track.index.lst[[i]]=list()
        length(comp.track.index.lst[[i]])=dim(comp[[i]])[2]
        names(comp.track.index.lst[[i]])=colnames(comp[[i]])

        for (j in 1:dim(comp[[i]])[2]){

            # select likelihood > 0.9
            comp.track.index.lst[[i]][[j]]=which(comp[[i]][,j]>=likelihood)
        }
    }


    # print selected tracks on the console
    l.print=list()
    length(l.print)=length(comp.track.index.lst)
    names(l.print)=names(comp.track.index.lst)
    for (i in 1:length(comp.track.index.lst)){
        l.print[[i]]=lapply(comp.track.index.lst[[i]],length)
    }
    l.print.result=do.call(rbind,l.print)

    cat("\n","at likelihood of",likelihood,", the number of trajectories selected are:\n")
    print(l.print.result)


    # use comp.track.index.lst to extract track index from dcoef

    comp.trackID.lst=list()
    length(comp.trackID.lst)=length(comp.track.index.lst)
    names(comp.trackID.lst)=names(comp.track.index.lst)


    # log.transformed=fit$log.transformed
    # if log.transformed, remove negative
    # values to make the total index corresponding to..

    if (log.transformed==T){
        dcoef.log.trans=lapply(dcoef,function(x){x=x[(x[,"slope"]>0),]})
        dcoef=dcoef.log.trans
    }

    # dcoef is supposed to have the same structure as comp.track.index.list
    # need a sanity check here TODO

    for (i in 1:length(dcoef)){

        comp.trackID.lst[[i]]=list()
        length(comp.trackID.lst[[i]])=length(comp.track.index.lst[[i]])
        names(comp.trackID.lst[[i]])=names(comp.track.index.lst[[i]])

        for (j in 1:length(comp.track.index.lst[[i]])){

            # subset the trackID using comp.track.index.lst[[i]][[j]]
            comp.trackID.lst[[i]][[j]]= dcoef[[i]][comp.track.index.lst[[i]][[j]],]
        }
    }


    ## (optional) remove zero length list

    # remove empty list out of list
    # Filter(function(k) length(k)>0, mylist)

    # remove na list out of list
    # Filter(function(x) !is.na(x), l)

    comp.trackID.lst=lapply(comp.trackID.lst,function(x){
        Filter(function(k) length(k)>0, x)
    })


    #---------------------------------------------------------------------------
    # export as index file to plot individually

    if (output==T){
        for (i in 1:length(comp.trackID.lst)){
            for (j in 1:length(comp.trackID.lst[[i]])){
                fileName=paste("componentTrackID-",
                               paste(names(comp.trackID.lst)[[i]],".",
                                     names(comp.trackID.lst[[i]][j]),sep=""),
                               .timeStamp("-"),".csv",sep="")
                write.csv(file=fileName,comp.trackID.lst[[i]][[j]])
            }
        }

    }

    #---------------------------------------------------------------------------
    # subset selected component tracks from original trackll for msd, dwell time
    # and further analysis

    # extract trackID

    trackID=list()
    length(trackID)=length(comp.trackID.lst)
    names(trackID)=names(comp.trackID.lst)

    for (i in 1:length(comp.trackID.lst)){
        trackID[[i]]=lapply(comp.trackID.lst[[i]],function(x){rownames(x)})
    }

    componentTracks=list()
    length(componentTracks)=length(trackID)
    names(componentTracks)=names(trackID)


    for (i in 1:length(trackID)){

        ## creates internal list structure of componentTracks
        componentTracks[[i]]=list()
        length(componentTracks[[i]])=length(trackID[[i]])
        names(componentTracks[[i]])=names(trackID[[i]])


        for (j in 1:length(trackID[[i]])){
            componentTracks[[i]][j]=lapply(trackll[i],function(x){x[ trackID[[i]][[j]]  ]})
        }

    }

    # remove empty list out of list
    # Filter(function(k) length(k)>0, mylist)

    # remove na list out of list
    # Filter(function(x) !is.na(x), l)

    # remove zero length list
    componentTracks=lapply(componentTracks,function(x){
        Filter(function(k) length(k)>0, x)
    })

    return(componentTracks)
}