RGLab / flowStats

flowStats: algorithms for flow cytometry data analysis using BioConductor tools
15 stars 10 forks source link

warpset cutting off low-end or high end data #48

Open The1stMartian opened 1 year ago

The1stMartian commented 1 year ago

Hi, I'm having an issue where warpset is changing values at the lowest end of certain samples, making them a single (lower) value. What is the cause of this and how do I prevent it from happening?

To see if some of the lower values were simply below some threshold, I added a single value to every cell in every channel, effectively shifting everything to the right, then I did the normalization again. This time the lower values were spared, while the upper values were cut off and moved to a single max value. Again, I'm not sure what I need to adjust to avoid this.

shiftComp

The command I've been using is "flowStats::warpSet(x = fixed, stains = fluorophores2, clipRange = 0)". (I set clipRange to zero thinking that this option could be to blame, but I essentially see the same thing when I leave it as default.)

mikejiang commented 1 year ago

reproducible example script with data, please

The1stMartian commented 1 year ago

reproducible example script with data, please

Hi Mike, I'm afraid I can't upload the data as it's not mine, but here's the script.

inputFolder = "/input" outputFolder = "/output" flowSetFolder = "/FlowSet" metaDataFile = "FCMetadata.csv" normFolder = "NormalizationData"

suppressMessages(library(tidyverse)) suppressMessages(library(flowCore)) suppressMessages(library(flowAI)) suppressMessages(library(flowAssist))

dir.create(file.path(outputFolder, flowSetFolder), recursive = TRUE, showWarnings = FALSE) dir.create(file.path(outputFolder, normFolder), recursive = TRUE, showWarnings = FALSE)

myfiles <- list.files(path=inputFolder, pattern = ".FCS", ignore.case = TRUE) input.FCS <- read.flowSet(myfiles, path=inputFolder, truncate_max_range = FALSE)

input.FCS.comp <-compensate(input.FCS, spillover(input.FCS[[1]])$SPILL) rm(input.FCS)

input.FCS.comp.cleaned <- flow_auto_qc(input.FCS.comp) rm(input.FCS.comp)

trans <- estimateLogicle(input.FCS.comp.cleaned[[1]], fluorophores) input.FCS.comp.cleaned.trans <- transform(input.FCS.comp.cleaned, trans) rm(input.FCS.comp.cleaned)

removeList = c("FSC-A", "FSC-H", "FSC-W", "SSC-A", "SSC-H","SSC-W", "Time") fluorophores<-as.character(colnames(exprs(input.FCS.comp.cleaned.trans[[1]]))) fluorophores<-fluorophores[!(fluorophores %in% removeList)]

normalizedData = flowStats::warpSet(x = input.FCS.comp.cleaned.trans, stains = fluorophores, clipRange = 0)

png(file.path(outputFolder, normFolder, "FluorDensity_Pre.png"), width = 1600, height = 1600) densityplot(~ ., input.FCS.comp.cleaned.trans) dev.off()

png(file.path(outputFolder, normFolder, "FluorDensity_Post.png"), width = 1600, height = 1600) densityplot(~ ., normalizedData) dev.off()

multi_gs <- GatingSet(normalizedData)

marker1 = "FSC-A" protein1 = "Forward Scatter" marker2 = "SSC-A" protein2 = "Side Scatter" currentPopName = "non_debris"

fs_data <- gs_pop_get_data(multi_gs)

non_debris <- fsApply(fs_data, function(fr) openCyto:::.flowClust.2d(fr, channels = c("FSC-A","SSC-A")))

non_debris <- rectangleGate(filterId="Non_debris_manual", "FSC-A" = c(0.25E5,3E5), "SSC-A" = c(10,3E5))

gs_pop_add(multi_gs, non_debris, parent = "root", name="non_debris")

gs_pop_remove(multi_gs, "non_debris")

recompute(multi_gs)

plt = autoplot(multi_gs, x=marker1, y=marker2, currentPopName, bins=2012) + labs(x= marker1, y = marker2) + theme(text=element_text(size = 16)) plt_path = file.path(inputFolder, gatingFolder, paste0(currentPopName,".png")) ggsave(plt_path, plt, width = 30, height = 30) plt

marker1 = "SSC-A" marker2 = "SSC-H" currentPopName = "singlets" parentPopName = "non_debris" currentPopName = "singlets"

fs_data <- gs_pop_get_data(multi_gs, parentPopName) singlets <- fsApply(fs_data, function(fr) openCyto:::.singletGate(fr, channels =c(marker1, marker2))) gs_pop_add(multi_gs, singlets, parent = parentPopName, name = currentPopName) recompute(multi_gs)

plt = autoplot(multi_gs, x = marker1, y = marker2, currentPopName, bins = 256) +labs(x= marker1, y = marker2) plt_path = file.path(outputFolder, paste0(currentPopName, ".png")) ggsave(plt_path, plt, width = 30, height = 30) print("Finished saving singlets plot...") plt

marker1 = "SSC-A" protein1 = "Side Scatter" marker2 = "BV570-A" protein2 = "Vital Stain" currentPopName = "live" parentPopName = "singlets"

fs_data <- gs_pop_get_data(multi_gs, parentPopName)

x_boundary = c(10,2E5) y_boundary = c(0,3.5) live <- rectangleGate(filterId="livegateID", "SSC-A" = x_boundary, "BV570-A" = y_boundary)

gs_pop_add(multi_gs, live, parent = parentPopName, name=currentPopName)

recompute(multi_gs)

plt = autoplot(multi_gs, x=marker1, y=marker2, currentPopName, bins=512) + labs(x = paste(marker1, protein1), y = paste(marker2, protein2)) + theme(text=element_text(size = 16)) plt_path = file.path(outputFolder, paste0(currentPopName,".png")) ggsave(plt_path, plt, width = 30, height = 30) print("Finished saving live/dead plot...")

marker1 = "BV570-A" protein1 = "Vital Stain" marker2 = "BUV496-A" protein2 = "CD45" parentPopName = "live" currentPopName = "CD45+"

fs_data <- gs_pop_get_data(multi_gs, parentPopName)

x_boundary = c(1, 4) y_boundary = c(2.5,5) CD45pos <- rectangleGate(filterId="CD45", "BV570-A" = x_boundary, "BUV496-A" = y_boundary)

gs_pop_add(multi_gs, CD45pos, parent = parentPopName, name=currentPopName) recompute(multi_gs)

plt = autoplot(multi_gs, x=marker1, y=marker2, currentPopName, bins=512) + ggcyto_par_set(limits = list(x = c(0, 5), y= c(0, 5))) + labs(x = paste(marker1, protein1, "\n\nPopulation:", currentPopName), y = paste(marker2, protein2)) + theme(text=element_text(size = 16)) plt_path = file.path(outputFolder, paste0(currentPopName,".png")) ggsave(plt_path, plt, width = 30, height = 30) print("Finished saving CD45 plot...")

marker1 = "BUV805-A"
protein1 = "CD16" marker2 = "BUV395-A"
protein2 = "CD56" currentPopName = "CD16+" parentPopName = "CD45+"

CD16pos <- rectangleGate(filterId="CD16_neg", "BUV805-A" = c(2, 6), "BUV395-A" = c(0,5))

gs_pop_add(multi_gs, CD16pos, parent = parentPopName, name=currentPopName) recompute(multi_gs)

plt = autoplot(multi_gs, x=marker1, y=marker2, currentPopName, bins=512) + ggcyto_par_set(limits = list(x = c(0, 9), y= c(0, 9))) + labs(x = paste(marker1, protein1, "\n\nPopulation:", currentPopName), y = paste(marker2, protein2)) + theme(text=element_text(size = 16)) plt_path = file.path(outputFolder, paste0(currentPopName,".png")) ggsave(plt_path, plt, width = 30, height = 30)