-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgating.R
141 lines (115 loc) · 5.33 KB
/
gating.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
# Load packages
library(viridis)
library(tidyverse)
library(FCSplankton)
# set the path to the Project folder
setwd("PATH/TO/PROJECT/")
# Set analysis for unstained (TRUE) or stained samples (FALSE)
unstained <- TRUE
# load list of FCS files
if(unstained){folder <- "unstained"
}else{folder <- "stained"}
file_list <- list.files(paste0(folder,"/raw"), pattern = ".fcs$", full.names=T)
### Initialization ###
######################
# load one file
this_file <- file_list[1]
fcs <- read_influx(this_file, transformation = TRUE) # create a dataframe WITH log-amplified data)
# Rename PMTs
print(names(fcs))
id <- c(2,3,4,5) ## replace number by column indice of FSC, 692, 580 and 530 respectively
names(fcs)[id]
names_pmt <- names(fcs)[id] # original names of FSC, 692, 580 and 530 respectively
### Gating populations of interest ###
######################################
gating <- TRUE
summary_table <- NULL
# summary_table <- read_csv(paste0(folder, "/summary.csv")) # if you want to append results to an existing summary_table
# Path where to save Gating output
system(paste0("mkdir ",folder, "/gating"))
for (this_file in file_list[]){
# Read data
# this_file <- file_list[1]
print(paste("gating", this_file))
fcs <- read_influx(this_file, transformation = TRUE) # create a flowFrame WITH transformation (for log-amplified data)
# change header
id <- match(names_pmt, names(fcs))
names(fcs)[id] <- c("scatter", "red", "orange", "green")
# Load gating if exist
if(!gating){
previous <- sub("raw","gating",paste0(this_file, ".RData"))
if(file.exists(previous)) load(previous)
}
# Gates Beads
if(gating & folder == "unstained") gates.log <- set_gating_params(fcs, "beads", "scatter", "orange")
if(gating & folder == "stained") gates.log <- set_gating_params(fcs, "beads", "scatter", "red")
# Assign "beads" label to particles
fcs <- classify_fcs(fcs, gates.log[][1])
# Normalization
beads <- fcs[which(fcs$pop == "beads"),]
fcs$norm.scatter <- fcs$scatter / mean(beads$scatter)
fcs$norm.orange <- fcs$orange / mean(beads$orange)
fcs$norm.red <- fcs$red / mean(beads$red)
fcs$norm.green <- fcs$green / mean(beads$green)
# Gating population - WARNINGS: don't change population names!
if(gating & folder == "unstained"){
gates.log <- set_gating_params(fcs, "synecho", "norm.scatter", "norm.orange", gates.log)
gates.log <- set_gating_params(fcs, "prochloro", "norm.scatter", "norm.red", gates.log)
gates.log <- set_gating_params(fcs, "picoeuk", "norm.scatter", "norm.red", gates.log)
#gates.log <- set_gating_params(fcs, "small-picoeuk", "norm_scatter", "norm_red", gates.log)
#gates.log <- set_gating_params(fcs, "large-picoeuk", "norm_scatter", "norm_red", gates.log)
}
if(gating & folder == "stained"){
gates.log <- set_gating_params(fcs, "bacteria", "norm.scatter", "norm.green", gates.log)
}
# Apply gates and label particles according to 'gates.log'
fcs <- classify_fcs(fcs, gates.log)
# Save Gating
if(gating){
save(gates.log, file=sub("raw","gating",paste0(this_file, ".RData")))
}
# Save plot
if(gating & folder == "stained"){
png(sub("raw","gating",paste0(this_file, ".png")),width=12, height=12, unit="in", res=200)
par(mfrow=c(2,2), pty="s", cex=1.2, oma=c(0,0,1,0), mar=c(5,5,1,1))
plot_vct_cytogram(fcs, "norm.scatter","norm.red", ylab="red\n (normalized to beads)", xlab="scatter\n (normalized to beads)")
plot_vct_cytogram(fcs, "norm.scatter","norm.orange", ylab="orange\n (normalized to beads)", xlab="scatter\n (normalized to beads)")
plot_vct_cytogram(fcs, "norm.scatter","norm.green", ylab="green\n (normalized to beads)", xlab="scatter\n (normalized to beads)")
dev.off()
}
if(gating & folder == "unstained"){
png(sub("raw","gating",paste0(this_file, ".png")),width=12, height=6, unit="in", res=200)
par(mfrow=c(1,2), pty="s", cex=1.2, oma=c(0,0,1,0), mar=c(5,5,1,1))
plot_vct_cytogram(fcs, "norm.scatter","norm.red", ylab="red\n (normalized to beads)", xlab="scatter\n (normalized to beads)")
plot_vct_cytogram(fcs, "norm.scatter","norm.orange", ylab="orange\n (normalized to beads)", xlab="scatter\n (normalized to beads)")
dev.off()
}
# Aggregate statistics
stat_table <- NULL
for(population in unique(fcs$pop)){
#print(i)
p <- subset(fcs, pop == population)
count <- nrow(p)
if(count == 0) {
scatter <- 0
red <- 0
}else{
scatter <- round(mean(p$norm.scatter),6)
red <- round(mean(p$norm.red),6)
orange <- round(mean(p$norm.orange),6)
if(folder == "stained"){green <- round(mean(p$norm.green),6)}
else{green <- NA}
}
var <- cbind(population,count,scatter,red,orange,green)
stat_table <- rbind(stat_table, var)
}
table <- data.frame(cbind(stat_table, file=basename(this_file)))
# remove entries that already exist
id <- which(!is.na(match(summary_table$file,unique(table$file))))
if(length(id) > 0) summary_table <- summary_table[-id,]
# add replace with new entries
summary_table <- rbind(summary_table, table)
}
### save summary table ###
##########################
write_csv(summary_table,path=paste0(folder, "/summary.csv"))