Add files via upload

DESCRIPTION

@@ -1,12 +1,26 @@
 Package: LinRace
 Type: Package
-Title: What the Package Does (Title Case)
+Title: cell division history reconstruction of single cells using paired lineage barcode and gene expression data
 Version: 0.1.0
-Author: Who wrote it
+Author: Xinhai Pan
 Maintainer: The package maintainer <[email protected]>
-Description: More about what it does (maybe more than one line)
-    Use four spaces when indenting paragraphs within the Description.
-License: What license is it under?
+Description: Understanding how single cells divide and differentiate into different cell types in devel-7
+oped organs is one of the major tasks of developmental and stem cell biology. Recently, lineage tracing8
+technology using CRISPR/Cas9 genome editing has enabled simultaneous readouts of gene expressions9
+and lineage barcodes in single cells, which allows for the reconstruction of the cell division tree, and10
+even the detection of cell types and differentiation trajectories at the whole organism level. While11
+most state-of-the-art methods for lineage reconstruction utilize only the lineage barcode data, methods12
+that incorporate gene expression data are emerging, aiming to improve the accuracy of lineage recon-13
+struction. However, effectively incorporating the gene expression data requires a reasonable model on14
+how gene expression data changes along generations of divisions. Here, we present LinRace (Lineage15
+Reconstruction with asymmetric cell division model), a method that integrates the lineage barcode16
+and gene expression data using the asymmetric cell division model and infers cell lineage under a17
+framework combining Neighbor Joining and maximum-likelihood heuristics. On both simulated and18
+real data, LinRace outputs more accurate cell division trees than existing methods. Moreover, Lin-19
+Race can output the cell states (cell types) of ancestral cells, which is rarely performed with existing20
+lineage reconstruction methods. The information on ancestral cells can be used to analyze how a pro-21
+genitor cell generates a large population of cells with various functionalities. 
+License: MIT license
 Encoding: UTF-8
 LazyData: true
 Suggests: 

LinRace.Rproj

@@ -1,20 +1,20 @@
-Version: 1.0
-
-RestoreWorkspace: Default
-SaveWorkspace: Default
-AlwaysSaveHistory: Default
-
-EnableCodeIndexing: Yes
-UseSpacesForTab: Yes
-NumSpacesForTab: 2
-Encoding: UTF-8
-
-RnwWeave: Sweave
-LaTeX: pdfLaTeX
-
-AutoAppendNewline: Yes
-StripTrailingWhitespace: Yes
-
-BuildType: Package
-PackageUseDevtools: Yes
-PackageInstallArgs: --no-multiarch --with-keep.source
+Version: 1.0
+
+RestoreWorkspace: Default
+SaveWorkspace: Default
+AlwaysSaveHistory: Default
+
+EnableCodeIndexing: Yes
+UseSpacesForTab: Yes
+NumSpacesForTab: 2
+Encoding: UTF-8
+
+RnwWeave: Sweave
+LaTeX: pdfLaTeX
+
+AutoAppendNewline: Yes
+StripTrailingWhitespace: Yes
+
+BuildType: Package
+PackageUseDevtools: Yes
+PackageInstallArgs: --no-multiarch --with-keep.source

NAMESPACE

@@ -1,12 +1,16 @@
 # Generated by roxygen2: do not edit by hand
 
+export(DivideMut_DCLEAR)
+export(Generate_Mut_Unique)
 export(LikelihoodCal)
 export(LikelihoodCal2)
 export(LikelihoodCal_Dif)
 export(NJ_asym)
 export(NJ_asym_Cas_Dif)
+export(NJ_asym_DCLEAR)
 export(NJ_asym_Dif)
 export(NJ_asym_alter)
+import(DCLEAR)
 import(TreeTools)
 import(ape)
 import(castor)

R/TreeBuild_Cas.R

@@ -35,6 +35,22 @@ DivideMut_Cas <- function(X){
   return(list(tree_backbone,dt))
 }
 
+#' Generate unique barcode data for cas hybrid runs
+#'
+#' @param mut_dir directory to lineage barcode matrix
+
+#' @export
+Generate_Mut_Unique <- function(mut_dir) {
+  X <- read.csv(file = mut_dir,row.names = 1,stringsAsFactors = FALSE)
+
+  X_unique <- X[!duplicated(X),]
+  l <- dim(X_unique)[1]
+  rownames(X_unique) <- 1:l
+
+  mut_unique_dir <- paste0(gsub(".csv","",mut_dir),"_unique.csv")
+  write.csv(X_unique, mut_unique_dir)
+}
+
 Cas_greedy <- function(X,prob_table){
   if (nrow(X) == 1){return(rownames(X))}
   #if (nrow(X) == 1){return(paste(X, collapse = '|'))}
@@ -154,3 +170,46 @@ NJ_asym_Cas_Dif <- function(muts,states,counts,state_lineages,max_Iter = 200){
 
   tree_final <- ConstructTree(tree_backbone,subtree_list)
 }
+
+
+
+#' LinRace main function with Cassiopeia backbone: asymmetric division based Cassiopeia-hybrid
+#'
+#' @param muts lineage barcode matrix
+#' @param states cell states of single cells
+#' @param counts gene expression data of single cells
+#' @param state_lineages the lineages that makes the state network from the root state to leaf states
+#' @param max_Iter the maximum iterations for local search
+#' @import data.table
+#' @import phangorn
+#' @export
+NJ_asym_Cas_hybrid <- function(muts,states,counts,state_lineages,max_Iter = 200,tree_backbone){
+  labels <- rownames(muts)
+  returnList <- DivideMut_Cas(muts)
+
+  dt <- returnList[[2]]
+
+  subtree_list <- list()
+
+  dm <- DiffusionMap(data=log2(counts+1),n_pcs = 30)
+  prob_t <- dm@transitions
+
+  for (i in 1:nrow(dt)){
+    cellids <- unlist(dt$cellids[i])
+    muts_sub <- muts[cellids,]
+    counts_sub <- counts[cellids,]
+    states_sub <- states[cellids]
+    labels_sub <- labels[cellids]
+    if (length(labels_sub)>1){
+      #res <- FindExpTree_Dif(states,labels = labels_sub,state_lineages,muts,prob_t,maxIter = max_Iter)
+      res <- FindExpTree_Dif_Newick(states, labels = labels_sub, state_lineages, muts, prob_t, maxIter = max_Iter, lambda1 = 1,lambda2 = 1)
+      subtree_opt <- res[[1]]
+      subtree_opt$name <- as.character(i)
+      subtree_list[[length(subtree_list)+1]] <- subtree_opt
+    }else {
+      tree_backbone$tip.label[tree_backbone$tip.label==as.character(i)] <- labels_sub
+    }
+  }
+
+  tree_final <- ConstructTree(tree_backbone,subtree_list)
+}

R/TreeBuild_DCLEAR.R

@@ -0,0 +1,111 @@
+#' LinRace main function with DCLEAR: asymmetric division based DCLEAR-kmer
+#'
+#' @param muts lineage barcode matrix
+#' @param states cell states of single cells
+#' @param counts gene expression data of single cells
+#' @param state_lineages the lineages that makes the state network from the root state to leaf states
+#' @param max_Iter the maximum iterations for local search
+#' @param lambda1 hyperparameter for asymmetric division likelihood
+#' @param lambda2 hyperparameter for neighbor distance likelihood
+#' @import data.table
+#' @import phangorn
+#' @export
+NJ_asym_DCLEAR <- function(muts, states, counts, state_lineages, max_Iter = 200,lambda1 = 10,lambda2 = 1) {
+  labels <- rownames(muts)
+  X <- muts
+  #X <- BarToMut(muts)
+  returnList <- DivideMut_DCLEAR(X)
+
+  tree_backbone <- returnList[[1]]
+  dt <- returnList[[2]]
+  dt$merged <- FALSE
+
+  subtree_list <- list()
+
+  dm <- DiffusionMap(data = log2(counts + 1), n_pcs = min(nrow(counts),30))
+  prob_t <- dm@transitions
+
+  #names_subtree <- c()
+
+  #df_pairs <- FindPairs(tree_backbone, dt)
+  #res <- MergeDT(tree_backbone, df_pairs, dt)
+  #tree_backbone <- res[[1]]
+  #dt <- res[[2]]
+
+  for (i in 1:nrow(dt)) {
+    if (!dt$merged[i]) {
+      cellids <- unlist(dt$cellids[i])
+      muts_sub <- muts[cellids,]
+      counts_sub <- counts[cellids,]
+      states_sub <- states[cellids]
+      labels_sub <- labels[cellids]
+
+      if (length(labels_sub) > 1) {
+        #res <- FindExpTree_Dif(states,labels = labels_sub,state_lineages,muts,prob_t,maxIter = max_Iter)
+        res <- FindExpTree_Dif_Newick(states, labels = labels_sub, state_lineages, muts, prob_t, maxIter = max_Iter, lambda1 = lambda1,lambda2 = lambda2)
+
+        subtree_opt <- res[[1]]
+        subtree_opt$name <- as.character(i)
+        #names_subtree <- c(names_subtree,as.character(i))
+        subtree_list[[length(subtree_list) + 1]] <- subtree_opt
+      }else {
+        tree_backbone$tip.label[tree_backbone$tip.label == as.character(i)] <- labels_sub
+      }
+    }
+  }
+  #names(subtree_list) <- names_subtree
+
+  tree_final <- ConstructTree(tree_backbone, subtree_list)
+}
+
+#' Divide muts using DCLEAR-kmer
+#'
+#' @param muts lineage barcode matrix
+#' @import DCLEAR
+#' @import phangorn
+#' @export
+DivideMut_DCLEAR <- function(X) {
+  mut_table_str <- c()
+  for (i in 1:nrow(X)) {
+    barcode <- X[i,]
+    cs <- paste(X[i,], collapse = '|')
+    mut_table_str <- c(mut_table_str, cs)
+  }
+  dt <- as.data.table(mut_table_str)[, list(list(.I)), by = mut_table_str]
+  colnames(dt) <- c("barcode", "cellids")
+
+
+  X_unique <- X[!duplicated(X),]
+  l <- dim(X_unique)[1]
+  rownames(X_unique) <- 1:l
+
+  sim_tree <- list()
+  states_barcode <- unique(unlist(X_unique))
+  states_mutated <- states_barcode[states_barcode != "0"]
+  nstates <- length(states_mutated)
+  #overflow_states <- c()
+  #if (nstates > 26){
+  #  overflow_states <- states_mutated[27:nstates]
+  #}
+  for (i in 1:dim(X_unique)[1]){
+    vec <- as.character(X_unique[i,])
+    #vec[vec %in% overflow_states] <- "0"
+    sim_tree[[i]] <- vec
+    #names(sim_tree[[i]]) <- rownames(muts_leaves)[i]
+  }
+  #print(length(sim_tree))
+  names(sim_tree) <- rownames(X_unique)
+  #nmstrings <- c('0','-',c(1:100))
+  nmstrings <- unique(unlist(sim_tree))
+
+  sim_data<- phyDat(sim_tree,type = 'USER',levels = nmstrings)
+  sim_data <- subset(sim_data,1:length(sim_tree))
+
+  dist_kmer <- dist_replacement(sim_data,reps = 20L, k = 2L)
+  Treenj <- nj(dist_kmer)
+  tree_backbone <- multi2di(Treenj)
+  tree_backbone$edge.length <- rep(1, length(tree_backbone$edge.length))
+  #tree_backbone$tip.label <- 1
+
+  return(list(tree_backbone, dt))
+}