2_SPADE_10000.R

# Load the preprocessed data:
# See script_preprocessing.R
#    ff:            Compensated flowFrame
#    ff_t:          Compensated and logicle transformed flowFrame
#    manual:        Array with label for each cell
#    selected:      Array with TRUE/FALSE whether cell falls in single live 
#                   cells
#    gatingMatrix:  Matrix with rows corresponding to cells and a column for 
#                   each manual gate. Each column contains TRUE/FALSE values
#                   indicating whether the cells fall in the specific gate
#    colsToCluster: Columns to use for clustering
load("FR-FCM-ZZQY/21-10-15_Tube_028.Rdata")

# Load the plot settings
# See script_plotSettings.R
#    cellTypeColors,markerColors,markerNames,
#    circular_markerOrder,grid_markerOrder
load("plotSettings.Rdata")

# Load the SPADE library
library(spade)

# Set some parameters
SPADE_nClusters = 50

# Set seed for reproducable results
set.seed(42)
# Record start time
start <- Sys.time()
    # Run the SPADE algorithm
    out_dir <- "SPADE_10000"
    SPADE.driver("fcs/21-10-15_Tube_028_selected_10000.fcs",
                 out_dir=out_dir,
                 cluster_cols=colnames(ff)[colsToCluster], 
                 comp=FALSE, transforms=logicleTransform(),
                 k=SPADE_nClusters)
    
    # Plot the results
    mst_graph <- igraph:::read.graph(paste(out_dir,"mst.gml",
                                           sep=.Platform$file.sep),
                                     format="gml")
    layout <- read.table(paste(out_dir,"layout.table",sep=.Platform$file.sep))
    SPADE.plot.trees(mst_graph,out_dir,layout=as.matrix(layout),
                     out_dir=paste(out_dir,"pdf",sep=.Platform$file.sep))
# Record end time
t_SPADE_10000<- Sys.time() - start


# Fit the SPADE results in a FlowSOM object for the visualization 
# used in the paper.
library(FlowSOM)
spade_fcs <- read.FCS(paste0(out_dir,"/21-10-15_Tube_028_selected_10000.fcs.density.fcs.cluster.fcs"))
spade_res <- list("MST"=list(),"map"=list(),"data"=exprs(ff_t[selected,][1:10000,]))
spade_res$MST$graph <- mst_graph
spade_res$MST$l <- as.matrix(layout)
spade_res$map$mapping <- exprs(spade_fcs[,"cluster"])
vsize <- table(spade_res$map$mapping)/nrow(spade_res$map$mapping)
spade_res$MST$size <- (vsize/(max(vsize, na.rm = TRUE)) * 3 + 2) * 4 # same resizing as in SPADE library
spade_res$map$codes <- spade_res$map$medianValues <- t(sapply(seq_len(SPADE_nClusters), function(i) {
    apply(subset(spade_res$data, spade_res$map$mapping[, 1] == i), 
          2, median)
}))
spade_res$map$medianValues[is.nan(spade_res$map$medianValues)] <- 0
spade_res$map$grid <- matrix(1:SPADE_nClusters,nrow=SPADE_nClusters,ncol=1)


# Save results
save(t_SPADE_10000, spade_res, file="SPADE_10000.Rdata")

# Plot results
pdf("SPADE_10000.pdf",useDingbats = FALSE)
    PlotPies(UpdateNodeSize(spade_res,reset=T),manual[selected][1:10000],
         colorPalette = colorRampPalette(c("#FFFFFF",cellTypeColors)))
    
    # Plot individual markers
    par(mfrow=c(3,3))
    for(m in grid_markerOrder){
        channel <- colnames(ff)[m]
        PlotMarker(spade_res,channel,
                   colorPalette = colorRampPalette(c("#dddddd",markerColors[channel])),
                   main=markerNames[channel])
    }
    par(mfrow=c(1,1))
dev.off()