Visualization of gene expression with Nebulosa

# Data pre-processing
seurat.all <- readRDS("FinalMergedData-downsampled.rds")

# Data normalization
# Let’s use SCTransform to stabilize the variance of the data 
# by regressing out the effect of the sequencing depth from each cell
seurat.all <- SCTransform(seurat.all, verbose = FALSE)

# Visualize data with Nebulosa
# For usability, it resembles the FeaturePlot function from Seurat
# Let’s plot the kernel density estimate for ELAVL4 as follows
plot_density(seurat.all, "ELAVL4") + 
    facet_grid(seurat.all$Mt~seurat.all$Age)

Subclustering your CellType-Specific Seurat Object

# Specify your cell type 
my_celltype <- "Ast"
# Specify time point of interest
my_timept <- "2mo"

# Change the filepath in quotes "" below to match where you 
# have your celltype data saved
# The expected filepath is "~/AlzheimersDS/student-notebooks/FinalMergedData-my_celltype.rds"
# (where `my-celltype` is your celltype code)
my_celltype_data <- readRDS("/data/Alzheimers_DS/FinalMergedData-Ast.rds")

my_seurat <- subset(x = my_celltype_data, subset = Age == my_timept)

# Necessary data pre-processing steps for clustering:
# Use raw counts data for input to PCA
DefaultAssay(my_seurat) <- "RNA"
# create normalzied data slot from raw count matrix
my_seurat <- NormalizeData(my_seurat,normalization.method = "LogNormalize")
# calculate most variable genes for this assay
my_seurat <- FindVariableFeatures(my_seurat, selection.method = "vst", nfeatures = 2000)
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : pseudoinverse used at -2.2786
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : neighborhood radius 0.50287
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : reciprocal condition number 1.0681e-15
## Warning in simpleLoess(y, x, w, span, degree = degree, parametric =
## parametric, : There are other near singularities as well. 0.090619
# create scaled data slot from raw count matrix, to use for input to PCA
my_seurat <- ScaleData(my_seurat)
## Centering and scaling data matrix
my_seurat <- RunPCA(my_seurat)
## PC_ 1 
## Positive:  GNAS, JUND, UBE2M, INA, MAPT, TPT1, PCLO, CALY, TTC9B, IER2 
##     NEFM, VGF, UNC13A, SEC31A, TCEAL2, PKIA, NEAT1, HSP90AB1, XIST, BICD1 
##     SRRM4, RTN1, TCEAL5, MACF1, ARL4C, NRP2, AKAP9, CHGB, PEG10, SRPR 
## Negative:  TMSB4XP8, TMSB4X, RPL41, NTRK2, GPM6B, CXCR4, NPC2, NR2F1, RPL3P4, B2M 
##     ARID5B, CRISPLD1, MGST1, HES1, C1orf192, IL13RA2, ENPP2, PI15, PLP1, PLS3 
##     CRNDE, EDNRB, HIST1H4C, MIR4435-1HG, PALMD, C5orf38, FGFBP3, MINOS1P3, AC004540.4, FLRT3 
## PC_ 2 
## Positive:  MEIS2, STMN2, SOX4, DCX, GRIA2, DLX6-AS1, NRXN1, CELF4, NSG2, ANK3 
##     DLX5, MT-CO3, KIAA1107, MYT1L, SYT1, SCN2A, RUNX1T1, MAP1B, TTC9B, GABRA2 
##     ATP1B1, MT-CO1, SYT10, RAB3B, NEFL, MT-CO2, SBK1, PTPRT, RTN1, NCAM1 
## Negative:  CD63, FTL, MRPL17, ANXA5, SERPINH1, SAT1, P4HB, ALDOA, SLC2A1, NPC2 
##     DAD1, PLTP, TRMT112, SERF2, ANXA2, ARL1, KRT10, NANS, FTH1, PTGDS 
##     CRYAB, PGM3, TUBA4A, CIB1, SEC13, TTYH1, VIM, KDELR2, NTRK2, MRPL55 
## PC_ 3 
## Positive:  ARF4, RAB3A, CYCS, CDK5R2, KISS1R, NSG2, HIST1H4C, PHLDA2, STMN2, NEFM 
##     TMSB4XP8, SEC24D, INA, SYT5, UNC13A, SCG5, MZT1, SPINT2, RPL41, MFSD2A 
##     MIR7-3HG, CRYAB, FAM126B, TTC9B, TMC4, BAG5, TAF13, SCG3, RP11-727F15.9, POLE3 
## Negative:  FABP7, VIM, PTPRZ1, SFRP1, PTN, MOXD1, METRN, HOPX, C14orf23, TNC 
##     FAM181B, CD82, TTYH1, C1orf61, PLAT, CD81, EFNB2, PEA15, LAMP5, SLC1A3 
##     SCRG1, NKAIN4, FOXG1, FAM171A2, RORB, RPS12, LHX2, CDO1, MT3, SYNE2 
## PC_ 4 
## Positive:  ENO1, HMGB2, EIF4EBP1, C14orf23, RPS16, DDIT4, SLC2A1, LRRC59, TRIB3, EMX2 
##     INSIG2, ALDOA, SFRP1, SLC2A3, IGFBP2, ARF1, SLC16A1, KDELR2, CDCA7L, TRMT112 
##     C11orf24, PFKFB4, GSX2, LOX, HK2, P4HB, RN7SK, SLC1A5, RGS16, SSSCA1 
## Negative:  BCAN, TFPI, PI15, NOVA1, GPM6B, SPOCK2, NR2F1, CPE, GJA1, B2M 
##     PTGDS, NTRK2, TAGLN2, TPM2, NR2F2, FILIP1L, MT-ND3, MT-CO2, MT-ND4, PCSK1N 
##     MT-CYB, MT-CO1, LINC00152, SPOCK1, ARID5B, SPARCL1, PRSS23, NCKAP5, OLFML2B, MGST1 
## PC_ 5 
## Positive:  NEAT1, PLOD2, WSB1, RP3-368A4.6, IRX2, CTC-482H14.5, HILPDA, IER5L, JUND, SEMA5B 
##     P4HA1, ZMYM2, HRK, BHLHE40, HK2, NR2F2, RNF38, ZNF576, ZBED6, SNHG1 
##     VEGFA, FAM132A, TAF1D, CREBRF, TET1, G0S2, HRH3, SERPINE1, SHC2, STC2 
## Negative:  TPBG, RSPO2, RSPO3, AIF1, RSPO1, SERF2, CNTNAP2, BEX1, MT-ND2, TCEAL2 
##     TMSB10, CCDC80, MT-ND1, BARHL2, PRADC1, MT-ND4, CLU, TMEM38B, PON2, GCNT1 
##     ATP1B1, CLDN5, TPM1, NMU, RTN1, LEPREL4, DMRT3, S100A6, LSMD1, SCG5
pca_res <- my_seurat[['pca']]

my_seurat <- JackStraw(my_seurat,dims=20)
# Identify the most significant PCs to keep
my_seurat <- ScoreJackStraw(my_seurat, dims=1:20)
#  JackStraw Returns a Seurat object where 
# `object[['pca']]@jackstraw
# contains the PC scoring results with significant scores 
js_res <- my_seurat[['pca']]@jackstraw

PC_sig_thresh <- 1e-05

PC_scoring_results<- as.data.frame(js_res@overall.p.values)
PCs_to_keep <- PC_scoring_results %>% filter(Score <= PC_sig_thresh) %>% arrange(Score)

# Use the significant PCs identified from JackStraw analysis
# as input to clustering.
PCs <- PCs_to_keep$PC
# Run the 2-step clustering analysis
my_seurat <- FindNeighbors(my_seurat, dims = PCs)
## Computing nearest neighbor graph
## Computing SNN
my_seurat<- FindClusters(my_seurat, 
                         #resolution denotes bias towards number of communities
                         #algorithm should find, ranging from fewer(0) to more(1)
                         resolution = 0.5,
                          n.start = 10,
                          n.iter = 10,
                         # random seed for reproducibility
                         random.seed = 0)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 179
## Number of edges: 4548
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.7098
## Number of communities: 4
## Elapsed time: 0 seconds
# Set active identities of my subset seurat to the new 
# "seurat_clusters" attribute 
Idents(my_seurat) <- "seurat_clusters"

Interacive 3D plot of your CellType at a specific timepoint

# Before you plot, set the ranges of the axis you desire. This set axis range will be 
# present across all clusters, and plotly will not adjust for axis length anymore
# this axis length will persist even when selecting some clusters

# xaxis
axx <- list(
  nticks = 4,
  range = c(-10,10) #select range of xaxis
)

# yaxis
axy <- list(
  nticks = 4,
  range = c(-10,10) #select range of yaxis
)

#zaxis
axz <- list(
  nticks = 4,
  range = c(-10,10) #select range of zaxis
)


# Prepare a dataframe for cell plotting
plot.data <- FetchData(object = my_seurat, vars = c("PC_1", "PC_2", "PC_3", "seurat_clusters"))

# Make a column of row name identities (these will be your cell/barcode names)
plot.data$label <- paste(rownames(plot.data))


fig <- plot_ly(data = plot.data, 
        x = ~PC_1, y = ~PC_2, z = ~PC_3, 
        color = ~seurat_clusters, 
        colors = c("lightseagreen",
                   "gray50",
                   "darkgreen",
                   "red4",
                   "red",
                   "turquoise4",
                   "black",
                   "yellow4",
                   "royalblue1",
                   "lightcyan3",
                   "peachpuff3",
                   "khaki3",
                   "gray20",
                   "orange2",
                   "royalblue4",
                   "yellow3",
                   "gray80",
                   "darkorchid1",
                   "lawngreen",
                   "plum2",
                   "darkmagenta"),
        type = "scatter3d", 
        mode = "markers", 
        marker = list(size = 5, width=2), # controls size of points
        text=~label, #This is that extra column we made earlier for which we will use for cell ID
        hoverinfo="text") #When you visualize your plotly object, hovering your mouse pointer over a point shows cell names

fig_cube <- fig %>% layout(scene = list(xaxis=axx,yaxis=axy,zaxis=axz, aspectmode='cube')) # To maintain cubic aspect
fig_cube

Interacive gene-expression 3D plot

DefaultAssay(object = my_seurat) <- "RNA"

# create a dataframe
plot_data <- FetchData(object = my_seurat, vars = c("PC_1", "PC_2", "PC_3", "PTGDS"), slot = 'data')

# Say you want change the scale, so that every cell having an expression >1 will be one color
# Basically, you are re-adjusting the scale here, so that any cell having a certain expression will light up on your 3D plot

# First make another column in your dataframe, where all values above 1 are re-assigned a value of 1
# This information is stored in the 'changed' column of your dataframe
plot_data$changed <- ifelse(test = plot_data$PTGDS <1, yes = plot_data$PTGDS, no = 1)

# Add the label column, so that now the column has 'cellname-its expression value'
plot_data$label <- paste(rownames(plot_data)," - ", plot_data$PTGDS, sep="")


plot_ly(data = plot_data, 
        x = ~PC_1, y = ~PC_2, z = ~PC_3, 
        color = ~changed,
        opacity = .5,
        colors = c('darkgreen', 'red'), 
        type = "scatter3d", 
        mode = "markers",
        marker = list(size = 5, width=2), 
        text=~label,
        hoverinfo="text"
)

SCORPIUS: an unsupervised approach for inferring linear developmental chronologies from single-cell RNA sequencing data

Normalise with Seurat

my_seurat <- NormalizeData(my_seurat)

Fetch the expression data from the Seurat object as follows.

expression <- t(as.matrix(my_seurat@assays$RNA@data))

Also fetch some metadata from the Seurat object. Change group_name to whatever column in srt@meta.data you are interested in.

group_name <- my_seurat@meta.data$seurat_clusters
group_name <- as.factor(group_name)

Reduce dimensionality of the dataset

SCORPIUS uses Torgerson multi-dimensional scaling to reduce the dataset to three dimensions. This technique attempts to place the cells in a space such that the distance between any two points in that space approximates the original distance between the two cells as well as possible.

The distance between any two samples is defined as their correlation distance, namely 1 - (cor(x, y)+1)/2. The reduced space is constructed as follows:

space <- reduce_dimensionality(expression, dist = "spearman", ndim = 3)

The new space is a matrix that can be visualized with or without colouring of the different cell types.

draw_trajectory_plot(space, progression_group = group_name, contour = TRUE)

Inferring a trajectory through the cells

The main goal of SCORPIUS is to infer a trajectory through the cells, and orden the cells according to the inferred timeline.

SCORPIUS infers a trajectory through several intermediate steps, which are all executed as follows:

traj <- infer_trajectory(space)

The result is a list containing the final trajectory path and the inferred timeline for each sample time.

The trajectory can be visualized with respect to the samples by passing it to draw_trajectory_plot:

draw_trajectory_plot(
  space, 
  progression_group = group_name,
  path = traj$path,
  contour = TRUE
)

Finding candidate marker genes

We search for genes whose expression is seems to be a function of the trajectory timeline that was inferred, as such genes might be good candidate marker genes for dendritic cell maturation.

gimp <- gene_importances(expression, traj$time, num_permutations = 0, num_threads = 8)
gene_sel <- gimp[1:50,]
expr_sel <- expression[,gene_sel$gene]

To visualise the expression of the selected genes, use the draw_trajectory_heatmap function.

draw_trajectory_heatmap(expr_sel, traj$time, group_name)

Finally, these genes can also be grouped into modules as follows:

modules <- extract_modules(scale_quantile(expr_sel), traj$time, verbose = FALSE)
draw_trajectory_heatmap(expr_sel, traj$time, group_name, modules)

module <- modules %>%
  dplyr::select(-orig_index)
knitr::kable(module)
feature module within_module_ordering
NTRK2 1 0.0000000
GPM6B 1 0.0909091
CXCR4 1 0.1818182
NPC2 1 0.2727273
TMSB4XP8 1 0.3636364
TMSB4X 1 0.3636364
RPL41 1 0.5454545
RP11-234A1.1 1 0.6363636
RP11-742N3.1 1 0.7272727
AC007969.5 1 0.8181818
RPL13P12 1 0.9090909
RPL3P4 1 1.0000000
RPL36A 2 0.0000000
RPL23A 2 0.0000000
RPL27A 2 0.0952381
RPL31 2 0.1428571
RPL13A 2 0.1904762
RPL26 2 0.2380952
RPS29 2 0.2857143
RPL21 2 0.3333333
RPS27 2 0.3809524
RPL15 2 0.4285714
RPL7 2 0.4761905
RPL34 2 0.5238095
RPL39 2 0.5714286
RPS25 2 0.5714286
RPS3A 2 0.5714286
RPS18 2 0.7142857
RPL35 2 0.7619048
TMA7 2 0.8095238
USMG5 2 0.8571429
RPS10 2 0.9047619
LAMTOR5 2 0.9523810
SHFM1 2 1.0000000
MT-CO1 3 0.0000000
MT-CO2 3 0.1428571
MT-ATP6 3 0.2857143
MT-CO3 3 0.4285714
MT-CYB 3 0.5714286
MT-ND5 3 0.7142857
MT-RNR2 3 0.8571429
MT-RNR1 3 1.0000000
STMN2 4 0.0000000
AKAP9 4 0.1428571
ARGLU1 4 0.2857143
DDX17 4 0.4285714
BPTF 4 0.5714286
JUND 4 0.7142857
TPT1 4 0.8571429
GNAS 4 1.0000000

Timepoint-Specific Clustering

# Load ExDp1 data
Ast_data <- readRDS("/data/Alzheimers_DS/FinalMergedData-Ast.rds")

# Necessary data pre-processing steps for clustering:
# Use raw counts data for input to PCA
DefaultAssay(Ast_data) <- "RNA"
# create normalzied data slot from raw count matrix
Ast_data <- NormalizeData(Ast_data,normalization.method = "LogNormalize")
# calculate most variable genes for this assay
Ast_data <- FindVariableFeatures(Ast_data, selection.method = "vst", nfeatures = 2000)
VariableFeaturePlot(Ast_data)

# create scaled data slot from raw count matrix, to use for input to PCA
Ast_data <- ScaleData(Ast_data)
## Centering and scaling data matrix
Ast_data <- RunPCA(Ast_data, features = VariableFeatures(object = Ast_data))
## PC_ 1 
## Positive:  STMN2, SOX4, SOX11, NSG2, DCX, NSG1, SYT1, BCL11A, PCLO, TUBB3 
##     ELAVL4, INA, TMSB15A, RAB3A, NEFL, INSM1, TMSB4XP8, TTC9B, BBC3, GADD45G 
##     MYT1L, SRRM4, ATCAY, MEF2C, NEUROD2, EIF4EBP1, SPINT2, NEUROD6, SYT4, CELF4 
## Negative:  CLU, AQP4, APOE, BCAN, HEPN1, CTSH, PLTP, B2M, METRN, ATP1A2 
##     F3, MT3, ZFP36, MLC1, ITM2B, ITM2C, CSPG5, ADD3, PLA2G16, SPARCL1 
##     GEM, MGST1, GJA1, C21orf62, S100A10, CD81, PON2, HEPACAM, SEPW1, HADHB 
## PC_ 2 
## Positive:  NRXN1, HOPX, MEIS2, GRIA2, SLC1A2, PTPRZ1, RTN1, RORB, FABP7, SFRP1 
##     DCX, MOXD1, APC, NCAN, ROBO2, SLC1A3, MGAT4C, NNAT, MT-ATP6, PMP2 
##     MT-CO1, MT-CYB, SCRG1, CITED1, KAL1, F3, C1orf61, DBI, CHRDL1, FABP5 
## Negative:  FTL, SLC2A1, SERPINH1, STC2, DDIT3, CTSD, RP3-412A9.11, FTH1, ANXA2, HSPA9 
##     CRABP2, S100A11, P4HA1, P4HB, NUPR1, SLIT2, NPC2, SVIP, ALDOA, FAM162A 
##     SCD, TRIB3, PRSS23, LINC00152, ENO1, C5orf38, MIR4435-1HG, MGST1, UPP1, INSIG2 
## PC_ 3 
## Positive:  CNTNAP2, ATP1B1, RTN1, STMN2, SIX3, TCEAL2, ARL4C, S100A10, RP11-89K21.1, EMX2 
##     KIAA1598, CTSH, ZMYND10, LEMD1, CALB1, INA, BCL11A, TTC9B, GADD45G, NSG1 
##     TCTEX1D1, WDR16, FEZF2, SRRM4, NSG2, ZIC4, NRN1, PCP4, DCX, SPINT2 
## Negative:  GPM6B, RAB31, PI15, NTRK2, AC004540.4, NR2F1, PTN, NR2F2, S100A13, DBI 
##     FAM89A, PDPN, RAMP1, S100A16, ENPP2, CRNDE, PMP2, IRX1, KCNIP1, TFPI 
##     TTYH1, MFGE8, LIX1, ARID5B, BCAN, FGFR3, BBOX1, C5orf38, PPP1R14C, CST3 
## PC_ 4 
## Positive:  SFRP1, MOXD1, FZD8, CHRDL1, LHX2, MOB3B, GLI3, SLC16A1, C8orf4, FABP7 
##     AKAP12, HOPX, NRG1, AMMECR1, CROT, SGCG, PGK1, ALDOA, TNC, LRRTM3 
##     LAMP5, BNIP3, PTN, TRIB3, PTH2R, IQGAP2, GRB14, LOX, HSPB1, FTL 
## Negative:  STMN2, PCSK1N, SCG3, NELL2, MT-ND2, NR2F1, MT-ND4, NSG1, GABBR2, SMOC1 
##     SYT1, DCX, SOX4, MT-CO1, MT-CYB, GRIA2, SPARCL1, BBOX1, TTC9B, NSG2 
##     PI15, MT-ATP6, SNAP25, MT-RNR1, MGST1, SH3BP5, MT-ND5, NEUROD6, PDE1A, PLP1 
## PC_ 5 
## Positive:  C1orf192, ZMYND10, TMSB4XP8, TCTEX1D1, FAM183A, C1orf194, CAPSL, WDR16, PIFO, C7orf57 
##     C5orf49, LGALS1, RP11-89K21.1, CTSH, DTHD1, ZIC4, RP11-71N10.1, CCDC39, ARMC3, NKX6-2 
##     SIX3, EFEMP1, RSPH1, MORN5, KIF9, TEKT1, CALB1, SIX3-AS1, RP11-415A20.1, RP11-159K7.2 
## Negative:  CAMK2N1, PCLO, BBC3, TXNRD1, IER2, SYT1, PKIA, AKAP12, INA, FARP1 
##     JUND, HSPA9, KCNQ1OT1, NEAT1, BCL11A, SLC38A1, SNAP25, CDKN1A, GADD45G, TTC9B 
##     MDM2, SEZ6L2, PLK2, TRIB3, SRRM4, ATF5, MAP1LC3B, CTTNBP2, TKT, SH3GL2
DimHeatmap(Ast_data, dims = 1, cells = 500, balanced = TRUE)

Ast_data <- JackStraw(Ast_data, num.replicate = 100)
# Identify the most significant PCs to keep
Ast_data_data <- ScoreJackStraw(Ast_data, dims=1:20)
ElbowPlot(Ast_data_data)

# Run the 2-step clustering analysis
Ast_data_data <- FindNeighbors(Ast_data_data, dims = 1:15)
## Computing nearest neighbor graph
## Computing SNN
Ast_data_data <- FindClusters(Ast_data_data, 
                         #resolution denotes bias towards number of communities
                         #algorithm should find, ranging from fewer(0) to more(1)
                         resolution = 0.5)
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 681
## Number of edges: 18779
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8473
## Number of communities: 8
## Elapsed time: 0 seconds
# Set active identities of my subset seurat to the new 
# "seurat_clusters" attribute 
head(Idents(Ast_data_data), 5)
## GTGTTAGTCGTACGGC_1_1 TGGCCAGTCCCTAATT_1_1 GCGACCAGTACAAGTA_1_1 
##                    1                    1                    1 
## GTCGTAATCACGCGGT_1_1 GGCTCGACACATCTTT_1_1 
##                    1                    1 
## Levels: 0 1 2 3 4 5 6 7
Idents(Ast_data_data) <- "seurat_clusters"

my_Ast_data_data <- subset(x =Ast_data_data, subset = Age == "2mo")

featuredist <- VlnPlot(my_Ast_data_data, features = c("nFeature_RNA"), split.by = "Mt") +
  ggtitle("Distributions of gene read counts in subclusters of Ast cells at 2mo") +
  xlab("Subcluster") + theme(plot.title = element_text(size=12))
## The default behaviour of split.by has changed.
## Separate violin plots are now plotted side-by-side.
## To restore the old behaviour of a single split violin,
## set split.plot = TRUE.
##       
## This message will be shown once per session.
my_Ast_data_data_2 <- subset(x =Ast_data_data, subset = Age == "4mo")

featuredist_2 <- VlnPlot(my_Ast_data_data_2, features = c("nFeature_RNA"), split.by = "Mt") +
  ggtitle("Distributions of gene read counts in subclusters of Ast cells at 4mo") +
  xlab("Subcluster") + theme(plot.title = element_text(size=12))

my_Ast_data_data_3 <- subset(x =Ast_data_data, subset = Age == "6mo")

featuredist_3 <- VlnPlot(my_Ast_data_data_3, features = c("nFeature_RNA"), split.by = "Mt") +
  ggtitle("Distributions of gene read counts in subclusters of Ast cells at 6mo") +
  xlab("Subcluster") + theme(plot.title = element_text(size=12))

featuredist

featuredist_2

featuredist_3

# find markers for every cluster compared to all remaining cells, report
# only the positive ones
Ast.markers <- FindAllMarkers(object = Ast_data_data, only.pos = TRUE, min.pct = 0.25, 
    thresh.use = 0.25)
## Calculating cluster 0
## Calculating cluster 1
## Calculating cluster 2
## Calculating cluster 3
## Calculating cluster 4
## Calculating cluster 5
## Calculating cluster 6
## Calculating cluster 7
result <- Ast.markers %>% group_by(cluster) %>% top_n(2, avg_log2FC) %>% dplyr::select(cluster, gene)
marker.result <- aggregate(gene ~ cluster, data = result, paste, collapse = ",")
colnames(marker.result) <- c("Cluster ID", "Marker genes")
knitr::kable(marker.result)
Cluster ID Marker genes
0 TMSB4XP8,CXCR4
1 SFRP1,C14orf23
2 STMN2,SOX4
3 CST3,PTGDS
4 CTSH,APOE
5 EVL,CRMP1
6 HOPX,SFRP1
7 RSPO2,RSPO3