Commit 8709a9bb authored by Jun Zhao's avatar Jun Zhao
Browse files

automate marker detection

parent 9584b374
Loading
Loading
Loading
Loading
+32 −2
Original line number Diff line number Diff line
@@ -132,7 +132,7 @@ getDAregion <- function(
    X.region.plot <- plotCellLabel(
      X = plot.embedding[cell.idx,], label = as.factor(X.tclust), 
      size = size, do.label = F, return.plot = T
    ) + scale_color_manual(values = c("gray",hue_pal()(X.n.da)), breaks = c(1:X.n.da))
    ) + scale_color_manual(values = c(rgb(255,255,255,max = 255,alpha = 0),hue_pal()(X.n.da)), breaks = c(1:X.n.da))
  } else {
    X.region.plot <- NULL
  }
@@ -146,6 +146,36 @@ getDAregion <- function(



## Step 3: detect genes that characterize DA regions from Step 2

#' @param cell.idx result "da.cell.idx" from the output of function getDAcells
#' @param da.region.label result "cluster.res" from the output of function getDAregion
#' @param obj Seurat object that contain ALL cells in the analysis
#' @param ... parameters for Seurat function FindMarkers()
#' 
#' @return a list of matrices with markers for each DA region

findMarkersForDAregion <- function(
  cell.idx, da.region.label, obj, ...
){
  n.da <- length(unique(da.region.label)) - 1
  obj@meta.data$da <- 0
  obj@meta.data$da[cell.idx] <- da.region.label
  obj <- SetAllIdent(obj, id = "da")
  
  da.markers <- list()
  for(ii in 1:n.da){
    da.markers[[ii]] <- FindMarkers(obj, ident.1 = ii, ...)
    da.markers[[ii]]$pct.diff <- da.markers[[ii]]$pct.1 - da.markers[[ii]]$pct.2
  }
  
  return(da.markers)
}





##=======================================================##
## Other functions

@@ -292,7 +322,7 @@ plotCellLabel <- function(X, label, cell.col = NULL, size = 0.5, do.label = T, r
  # Plot cells with labels
  myggplot <- ggplot() + theme_cowplot() +
    geom_point(data = data.frame(Dim1 = X[,1], Dim2 = X[,2], Group = label, stringsAsFactors = F), 
               aes(x = Dim1, y = Dim2, col = Group), size = size, alpha = 0.75) + 
               aes(x = Dim1, y = Dim2, col = Group), size = size) + 
    guides(colour = guide_legend(override.aes = list(size=3), title = NULL))
  
  # Change cell color
+11 −18
Original line number Diff line number Diff line
@@ -120,7 +120,9 @@ da_cells$da.cells.plot

# cluster DA cells to get DA regions
da_regions <- getDAregion(
  X = pca_embedding, cell.idx = da_cells$da.cell.idx, k = 5, alpha = 0.1, 
  X = pca_embedding, 
  cell.idx = da_cells$da.cell.idx, 
  k = 5, alpha = 0.1, 
  cell.labels = data_S@meta.data$lesion,
  labels.1 = labels_res, 
  labels.2 = labels_nonres, 
@@ -130,23 +132,14 @@ da_regions$da.region.plot
da_regions$DA.stat



## Marker detection

# set ident in Seurat
n.da <- length(unique(da_regions$cluster.res)) - 1
data_S@meta.data$da <- 0
data_S@meta.data$da[da_cells$da.cell.idx] <- da_regions$cluster.res
data_S <- SetAllIdent(data_S, id = "da")
TSNEPlot(data_S, pt.size = 0.5)


# identify markers for each DA region
da.markers <- list()
for(i in 1:n.da){
  da.markers[[i]] <- FindMarkers(data_S, ident.1 = i, only.pos = T, min.pct = 0.1, min.diff.pct = 0.09)
  da.markers[[i]]$pct.diff <- da.markers[[i]]$pct.1 - da.markers[[i]]$pct.2
}
# Marker detection for DA regions
da_markers <- findMarkersForDAregion(
  cell.idx = da_cells$da.cell.idx,
  da.region.label = da_regions$cluster.res,
  obj = data_S,
  only.pos = T, min.pct = 0.1, min.diff.pct = 0.09
)
str(da_markers)


# plot some markers
+9 −0
Original line number Diff line number Diff line
@@ -17,6 +17,8 @@ DA-seq can be used as follows:

Let X be a N-by-p matrix of the PCA embeddings of merged scRNA-seq datasets A (A1 and A2) and B (B1 and B2); X.label be a vector of N specifying the original of each cell ('A1', 'A2', 'B1', or 'B2'); X.2d be the 2D embedding of the cells.

For marker detection of each DA region, let SeuratObj be a Seurat V2.3.0 ([tutorial](https://satijalab.org/seurat/v2.4/pbmc3k_tutorial.html)) object containing the merged datasets (N cells), after proper preprocessing steps (normalization, scaling, etc.).

~~~~
X.da.cells <- getDAcells(
  X = X, 
@@ -36,6 +38,13 @@ X.da.regions <- getDAregion(
  labels.2 = c("B1","B2"), 
  plot.embedding = X.2d
)

X.da.markers <- findMarkersForDAregion(
  cell.idx = X.da.cells$da.cell.idx,
  da.region.label = X.da.regions$cluster.res,
  obj = SeuratObj,
  only.pos = T, min.pct = 0.1, min.diff.pct = 0.09
)
~~~~