jgarces02 / flowct Goto Github PK

fsom.clustering <- function(fcs.SE, assay.i = "normalized", scale.data= TRUE, markers.to.use = rownames(fcs.SE), markers.to.plot = NULL, k.metaclustering = 40, metaclustering.name = NULL){
  require(FlowSOM)
  set.seed(1234)
  data <- as.flowSet.SE(fcs.SE, assay.i)
  
  ## FSOM clustering
  cat("Calculating SOM clustering...\n")
  fsom <- suppressMessages(ReadInput(data, transform = FALSE, scale = scale.data)) #read data
  fsom <- suppressMessages(BuildSOM(fsom, colsToUse = markers.to.use)) #build SOM
  
  cat("Building MST...\n")
  fsom <- suppressMessages(BuildMST(fsom, tSNE = TRUE, silent = T)) #build MST for visualization of clustering
  
  #plot the MST to evaluate the marker fluorescence (or general tree) intensity for each SOM
  if(!is.null(markers.to.plot)){
    if(markers.to.plot == "tree"){
      PlotStars(fsom)
    }else{
      for(marker in markers.to.plot) PlotMarker(fsom, marker)
    }
  }
  
  ## Metaclustering
  if(!is.null(k.metaclustering)){
    cat("Metaclustering...\n")
    if(!is.null(metaclustering.name)){
      mc <- suppressMessages(ConsensusClusterPlus::ConsensusClusterPlus(t(fsom$map$codes), maxK = k.metaclustering, reps = 100,
                                                                        pItem = 0.9, pFeature = 1, title = metaclustering.name, plot = "pdf",
                                                                        clusterAlg = "hc", innerLinkage = "average", finalLinkage = "average",
                                                                        distance = "euclidean", seed = 333, verbose = F))
    }else{
      mc <- suppressMessages(ConsensusClusterPlus::ConsensusClusterPlus(t(fsom$map$codes), maxK = k.metaclustering, reps = 100,
                                                                        pItem = 0.9, pFeature = 1, title = "consensus_plots", plot = "pdf",
                                                                        clusterAlg = "hc", innerLinkage = "average", finalLinkage = "average",
                                                                        distance = "euclidean", seed = 333, verbose = F))
      unlink("consensus_plots", recursive = TRUE)
    }
    
    #get cluster ids for each cell
    code_clustering1 <- mc[[k.metaclustering]]$consensusClass %>% as.factor()
    cell_clustering1 <- code_clustering1[fsom$map$mapping[,1]]
    
    #add clustering to original MST and color by cluster colors
    PlotStars(fsom, backgroundValues = code_clustering1, backgroundColor = alpha(div.colors(length(code_clustering1)), alpha = 0.7))
    
    return(list(fsom = fsom, metaclusters = cell_clustering1, plotStars_value = code_clustering1))
  }else{
    return(fsom)
  }
}

      pb <- txtProgressBar(min = 0, max = nrow(gene.sum), 
        style = 3)
      for (i in 1:nrow(gene.sum)) {
        prot.dat = all.prot.dat[Hugo_Symbol %in% gene.sum[i, 
          Hugo_Symbol]]
        syn.res = rbind(syn.res, cluster_prot(prot.dat = prot.dat, 
          gene = gene.sum[i, Hugo_Symbol], th = gene.sum[i, 
            th], protLen = gene.sum[i, aa.length]))
        setTxtProgressBar(pb, i)
      }

df<-as.data.frame(t(assay(fcs1000,"logcounts")))
df_expr<-cbind(df, SOM=fcs1000$SOM)
#or if only selected clusters
df_expr<-df_expr[df_expr$SOM==c("1","2","3","4"),]

#   prepare the dotplot
pmain<-ggplot(df_expr, aes(x=CD4, y=CD8, color=SOM)) +
  geom_point(size=0.5)+
  scale_color_manual(values = color_clusters) +
  theme_bw()
###if density plot needed
#ggplot(df_expr, aes(x=CD4, y=CD8)) +
  # geom_point(size=0.5)+theme_bw() +geom_hex(bins = 200) +
  # scale_fill_continuous(type = "viridis") 

library(cowplot) 

# histogram along x axis
xdens <- axis_canvas(pmain, axis = "x")+
  scale_fill_manual(values = color_clusters) +
  geom_density(data = df_expr, aes(x = CD4, fill = SOM),
               alpha = 0.7, size = 0.2)
# histogram along y axis
ydens <- axis_canvas(pmain, axis = "y", coord_flip = TRUE)+
  scale_fill_manual(values = color_clusters) +
  geom_density(data = df_expr, aes(x = CD8, fill = SOM),
               alpha = 0.7, size = 0.2)+
  coord_flip()
p1 <- insert_xaxis_grob(pmain, xdens, grid::unit(.2, "null"), position = "top")
p2<- insert_yaxis_grob(p1, ydens, grid::unit(.2, "null"), position = "right")
ggdraw(p2)

multidensity <- function(fcs.SCE, assay.i, show.markers = "all", color.by = NULL, subsampling = NULL, interactive = F, ridgeline.lim = 15, colors = NULL){
  if(show.markers == "all") show.markers <- rownames(fcs.SCE)
  if(!is.null(subsampling)) suppressMessages(fcs.SCE <- sub.samples(fcs.SCE, subsampling = subsampling))

  data <- t(assay(fcs.SCE, i = assay.i))
  data2 <- cbind(data, colData(fcs.SCE))

  # prepare tables: for plotting and with median values for each marker
  median_df <- data.frame(antigen = show.markers, median = apply(data[,show.markers], 2, median))
  ggdf <- as.data.frame(melt(as.data.table(data2), measure.vars = show.markers, value.name = "expression", variable.name = "antigen"))

  if(is.null(colors)) colors <- div.colors(unique(length(ggdf[,color.by])))

  if(length(unique(fcs.SCE$filename)) > ridgeline.lim){
    g <- ggplot(data = ggdf[grepl(paste0(show.markers, collapse = "|"), ggdf$antigen),],
                aes_string(x = "expression", color = color.by, group = "filename")) +
      # geom_density(size = 0.5) +
      stat_density(geom = "line", position = "identity", size = 0.5) +
      facet_wrap(~ antigen, scales = "free") +
      geom_vline(data = median_df, aes(xintercept = median), linetype = 2, color = "gray55") +
      scale_color_manual(name = color.by, values = colors) +
      theme_minimal() + theme(axis.text.x = element_text(angle = 90, hjust = 1),
                              strip.text = element_text(size = 7), axis.text = element_text(size = 5))

    if(interactive) ggplotly(g) else print(g)
  }else{
    suppressMessages(print(ggplot(ggdf, aes_string(x = "expression", y = "filename")) +
                             geom_density_ridges(alpha = 0.7) +
                             facet_wrap(~ antigen, scales = "free") +
                             geom_vline(data = median_df, aes(xintercept = median), linetype = 2, color = "gray55") +
                             theme_minimal() + theme(axis.text.x = element_text(angle = 90, hjust = 1),
                                                     strip.text = element_text(size = 7), axis.text = element_text(size = 7))))
  }
}

barplot.cell.pops_m <- function (cell_clusters, metadata, colname_sampleID, return_table)
{
  tab <- table(cell_clusters, metadata[,colname_sampleID])
  prop_tableL <- prop.table(tab, margin = 2) * 100
 
  ggdfL <- data.table::melt(prop_tableL, value.name = "proportion")
  colnames(ggdfL)[2] <- colname_sampleID
  mmL <- match(ggdfL[, colname_sampleID], metadata[, colname_sampleID])
  ggdfL <- data.frame(metadata[mmL, ], ggdfL)
  g <- ggplot(ggdfL, aes_string(x = "sample_id", y = "proportion",
                                fill = "cell_clusters")) + geom_bar(stat = "identity") +
    facet_wrap(~condition, scales = "free_x") + theme_minimal() +
    theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
    scale_fill_manual(values = colors_palette)
  if ("try-error" %in% class(suppressWarnings(try(x11(), silent = T)))) {
    cat("X11 is not active, boxplot is saved in -> barPlot_cell_prop.",
        deparse(substitute(cell_clusters)), ".jpg\n", sep = "")
    suppressMessages(ggsave(paste0("barPlot_cell_prop.",
                                   deparse(substitute(cell_clusters)), ".jpg"), device = "jpeg",
                            plot = g))
  }
  else {
    print(g)
  }
  if(return_table == "counts"){
    return(tab)
  }else{ #add "percentage" option and block any other possibilities
    return(prop_tableL)
  }
 
}

In data.table::melt(table(mnames)) :
  The melt generic in data.table has been passed a table and will attempt to redirect to the relevant reshape2 method; please note that reshape2 is deprecated, and this redirection is now deprecated as well. To continue using melt methods from reshape2 while both libraries are attached, e.g. melt.list, you can prepend the namespace like reshape2::melt(table(mnames)). In the next version, this warning will become an error.

> fcs <- normalization.flw(fcs.SCE = fcs, marker.to.norm = c("CCR6", "CCR4"),
+                          norm.method = "harmony", var.to.use = "patient_id")

error: Mat::init(): requested size is too large

 *** caught segfault ***
address (nil), cause 'memory not mapped'

Traceback:
 1: .External(list(name = "CppMethod__invoke_void", address = <pointer: 0x117fd830>,     dll = list(name = "Rcpp", path = "/home/cdasilvam/FlowCT.v2_environment1/renv/library/R-4.0/x86_64-pc-linux-gnu/Rcpp/libs/Rcpp.so",         dynamicLookup = TRUE, handle = <pointer: 0x1b6c8e90>,         info = <pointer: 0x19cfe20>), numParameters = -1L), <pointer: 0x29d889b0>,     <pointer: 0x361e00c0>, .pointer, ...)
 2: harmonyObj$setup(data_mat, phi, phi_moe, Pr_b, sigma, theta,     max.iter.cluster, epsilon.cluster, epsilon.harmony, nclust,     tau, block.size, lambda_mat, verbose)
 3: HarmonyMatrix(data_mat = assay(fcs.SCE, assay.i), meta_data = colData(fcs.SCE),     vars_use = var.to.use, do_pca = F, verbose = F)
 4: normalization.flw(fcs.SCE = fcs, marker.to.norm = c("CCR6", "CCR4"),     norm.method = "harmony", var.to.use = "patient_id")

Possible actions:
1: abort (with core dump, if enabled)
2: normal R exit
3: exit R without saving workspace
4: exit R saving workspace

>qc.and.removeDoublets(reduction_suffix="")
Processing: Screening_013543_TubeT.fcs
Processing: Screening_014090_TubeT.fcs
Processing: Screening_014184_TubeT.fcs
Processing: Screening_014408_TubeT.fcs
Processing: Screening_014497_TubeT.fcs
Processing: Screening_014786_TubeT.fcs
Processing: Screening_015308_TubeT.fcs
Processing: Screening_015879_TubeT.fcs
Processing: Screening_016286_TubeT.fcs
Processing: Screening_016469_TubeT.fcs
Processing: Screening_016692_TubeT.fcs
Processing: Screening_017425_TubeT.fcs
Processing: Screening_017499_TubeT.fcs
Processing: Screening_017529_TubeT.fcs
Processing: Screening_017557_TubeT.fcs
Processing: Screening_017759_TubeT.fcs
Processing: Screening_017789_TubeT.fcs
Processing: Screening_018108_TubeT.fcs
Processing: Screening_018769_TubeT.fcs
Processing: Screening_018796_TubeT.fcs
Processing: Screening_019278_TubeT.fcs
Processing: Screening_019649_TubeT.fcs
Processing: Screening_020616_TubeT.fcs
Processing: Screening_020653_TubeT.fcs
Processing: Screening_021098_TubeT.fcs
Error in cat("WARNING! >", i, "has lost some much cells (more that 30%) in the QC and doublets removal steps, consider to review it!",  :
  argument is missing, with no default
In addition: Warning message:
In dir.create(output_folder) : 'results_preprocessing' already exists

scale.data <- function(fcs.SE, assay.i = "normalized",scaled.matrix.name = "scaled"){
  data <- t(assay(fcs.SE, i = assay.i))
  rng <- matrixStats::colQuantiles(data, probs = c(0.01, 0.99))
  expr01 <- t((t(data) - rng[, 1]) / (rng[, 2] - rng[, 1]))
  expr01[expr01 < 0] <- 0
  expr01[expr01 > 1] <- 1
  SummarizedExperiment::assay(fcs.SE, i = scaled.matrix.name) <-t(expr01)
  return(fcs.SE)
  }

Error in `colnames<-`(`*tmp*`, value = new_names) : 
  attempt to set 'colnames' on an object with less than two dimensions

div.colors <- function(n, set.seed = 333){
  require(RColorBrewer)
  require(colorspace)
  
  if(n < 74){ 
    col<-qualitative_hcl(n, palette = "Dark3")  
  }else{
    qual_col_pals = sample(brewer.pal.info)
    col_vector = sample(unlist(mapply(brewer.pal, qual_col_pals$maxcolors, rownames(qual_col_pals))))
    col = sample(col_vector, 333)
    set.seed(set.seed); col <- sample(col, n, replace = T)
  }
  return(col)
}

#to be replaced
fcsL <- fcs1000[,fcs1000$SOM_named == "lymphocytes"]
metadata(fcsL)$subclustering <- "lymphocytes"

###new code to be included in the script
fcs$SOM_named <- clusters.rename(fcs$SOM, cluster = replacedata$original_cluster, name = replacedata$new_cluster)
fcsL<-subset(fcs, ,SOM_named=="lymphocytes")