Cytofkit

library(flowCore)
library(cytofkit)

# First run 
# cytofkit_GUI()
# compute the analysis
# then run the following command to view the results

cytofkitShinyAPP("C:/demo/190917-atelier/demo/cytofkit_demo.RData")

Cytofast

Cytofast has nice plots.

Read Cytofkit FCS export

library(flowCore)
library(cytofast)

cfData = readCytofkitFCS("demo/cytofkit_demo_analyzedFCS")
## Reading .FCS sample: 1
## Reading .FCS sample: 2
## Reading .FCS sample: 3
## Reading .FCS sample: 4
## Reading .FCS sample: 5
## Reading .FCS sample: 6
# cfData = readCytofkitFCS("demo/cytofkit_demo_analyzedFCS", clusterID = "Rphenograph_clusterIDs")

Optionnally, do some cleaning

# Remove unneeded channels
keep <- colnames(cfData@expr)  # all
keep <- keep[!grepl("^(Time|Event_length)", keep, ignore.case = TRUE)]
keep <- keep[!grepl("(_ADN1|_ADN2)", keep, ignore.case = TRUE)]
keep <- keep[!grepl("^X\\d+", keep) | grepl("_", keep)]  # unannotated X channel
cfData@expr <- cfData@expr[, keep]

# Rename channels
colnames(cfData@expr) <- gsub("(^X\\d+.+?_)", "", keep)  # remove metal tag

# Manual ordering of channels
first <- c("CD14", "CD19", "CD3", "CD4", "CD8", "CD56", "CD45RA", "CD95", "CD127")
keep <- colnames(cfData@expr)  # all
cfData@expr <- cfData@expr[, c(keep[1:2], first, setdiff(keep[-(1:2)], first))]

Read meta data that group FCS samples

# Create a template file to annotate samples
if (!file.exists("meta.csv")) {
  meta <- data.frame(cfData@samples, group = "")
  write.csv(meta, "meta.csv")
}

# Fill the template using Excel or so

# Import annotation
meta <- read.csv("meta.csv", row.names = 1)
meta <- meta[match(row.names(cfData@samples), meta[,"sampleID"]),] # match sampleID
cfData@samples <- cbind.data.frame(cfData@samples, meta[, -1, drop = FALSE])
cfData@samples
##                  sampleID               FCSfilename group
## HD1_Ungated   HD1_Ungated  cytofkit_HD1_Ungated.fcs    HD
## HD2_Ungated   HD2_Ungated  cytofkit_HD2_Ungated.fcs    HD
## HD3_Ungated   HD3_Ungated  cytofkit_HD3_Ungated.fcs    HD
## HIV1_Ungated HIV1_Ungated cytofkit_HIV1_Ungated.fcs   HIV
## HIV2_Ungated HIV2_Ungated cytofkit_HIV2_Ungated.fcs   HIV
## HIV3_Ungated HIV3_Ungated cytofkit_HIV3_Ungated.fcs   HIV

Transform expression values, in the same way you display marker intensity

# Transform expression with asinh( x / 5 )
cfData@expr[,-(1:2)] <- asinh(cfData@expr[,-(1:2)]/5)

Works with cell counts

# Retrieve cell counts
cfData <- cellCounts(cfData)
head(cfData@counts)
##                1    2   3   4  5   6   7   8   9  10  11  12   13 14  15
## HD1_Ungated   20  114 109  79 35  40 426 600 127   8 123  70  975 17   2
## HD2_Ungated    5  318 304 248 56  41 179 460 538   8 141 177  854 36  21
## HD3_Ungated    9  143 127  71 25  20 360 531  95   9 116  61  975 14   6
## HIV1_Ungated 324 3275 708 111 33 102   1 123 478 102  67 278  637 87 118
## HIV2_Ungated 172 1407 175 537 75 140  64 305 313  65 155 199 1130 30 147
## HIV3_Ungated 151 2069 325 752 95 176  43 210 277 136  83 110  487 21 120
##               16  17  18 19  20  21  22  23  24  25  26  27 28   29  30 31
## HD1_Ungated  391  46 211 61 193  59 486  93 156 235 102  38 71  115 200 13
## HD2_Ungated   50  11 253 26 191  23 750 151 119  63  35 194 71  652 325 10
## HD3_Ungated  321  90 174 51 183  47 520 109 145 246  84  17 84  127 190 31
## HIV1_Ungated  95   4 236 43  30 176 133 312 329  10  12 113 33   73 139 12
## HIV2_Ungated  82 144 247 18  55  55 430 404  86  37  14 146 48 1527 183 20
## HIV3_Ungated 129  10 240 44  55  14 634 178  98  66  26 325 20 1650  90 28
##              32 33 34   35 36   37   38  39  40
## HD1_Ungated  20 66 87 1195 42 1965 1207  66 137
## HD2_Ungated  49 28 27  966 38 1118 1242 146  76
## HD3_Ungated  35 32 31 1448 32 1889 1318  83 151
## HIV1_Ungated 42 24  7  583 97  548  397  18  90
## HIV2_Ungated 62 26 24  730 64  105  471  44  64
## HIV3_Ungated 18 17 14  435 77   90  576  24  87
cytoHeatmaps(cfData, group="group", legend=TRUE)

Cell counts need to be transformed

# Store raw counts
cellCountRaw <- cfData@counts

# Transform cell counts with log2( x + floor )
cfData@counts <- log2(cfData@counts+10)
# Center
cfData@counts <- sweep(cfData@counts, 2, colMeans(cfData@counts), "-")
round(head(cfData@counts))
##               1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
## HD1_Ungated  -1 -2 -1 -1  0 -1  2  1 -1 -1  0 -1  0 -1 -2  1  0  0  1  1
## HD2_Ungated  -2 -1  0  0  0 -1  1  0  1 -1  0  0  0  0 -1 -1 -1  0  0  1
## HD3_Ungated  -2 -2 -1 -1 -1 -1  2  1 -1 -1  0 -1  0 -1 -2  1  1  0  0  1
## HIV1_Ungated  2  2  2 -1  0  1 -3 -1  1  1 -1  1  0  1  1  0 -2  0  0 -1
## HIV2_Ungated  1  1  0  1  1  1 -1  0  0  1  0  1  0  0  2 -1  2  0 -1 -1
## HIV3_Ungated  1  2  0  2  1  1 -1 -1  0  2  0  0 -1  0  1  0 -1  0  0 -1
##              21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40
## HD1_Ungated   0  0 -1  0  2  1 -1  0 -2  0  0 -1  1  1  1  0  2  1  0  0
## HD2_Ungated  -1  1  0  0  0  0  1  0  1  1  0  0  0  0  0  0  1  1  1  0
## HD3_Ungated   0  0 -1  0  2  1 -2  1 -1  0  1  0  0  0  1 -1  2  1  1  1
## HIV1_Ungated  2 -2  1  1 -2 -1  0  0 -2  0  0  0  0 -1  0  1  0 -1 -1  0
## HIV2_Ungated  0  0  1 -1 -1 -1  1  0  2  0  0  1  0  0  0  0 -2 -1  0 -1
## HIV3_Ungated -1  1  0  0  0  0  2 -1  2 -1  0 -1 -1 -1 -1  0 -2  0 -1  0
cytoHeatmaps(cfData, group="group", legend=TRUE)

Let’s take HD as reference

# Center to HD
cfData@counts <- sweep(cfData@counts, 2, colMeans(cfData@counts[cfData@samples$group == "HD",]), "-")
round(head(cfData@counts)*10)
##               1  2  3  4  5  6   7   8  9 10 11 12 13 14 15  16  17 18 19
## HD1_Ungated   6 -6 -5 -5 -1  2   5   2 -5  0  0 -4  1 -2 -6  10   2  0  4
## HD2_Ungated  -4  8  9 11  5  3  -7  -2 15  0  1  9 -1  6  8 -17 -12  3 -6
## HD3_Ungated  -1 -3 -3 -6 -4 -5   2   0 -9  1 -1 -5  1 -4 -2   7  10 -3  2
## HIV1_Ungated 40 42 21  0 -1 14 -48 -20 13 26 -8 15 -5 16 28  -9 -18  2  0
## HIV2_Ungated 32 29  1 22  9 18 -21  -8  7 20  3 10  3  4 31 -11  17  2 -9
## HIV3_Ungated 30 35 10 26 12 21 -26 -13  5 30 -6  2 -9  0 28  -5 -13  2  0
##               20  21  22 23 24  25  26  27  28  29  30 31 32 33  34  35 36
## HD1_Ungated    0   4  -2 -3  2   6   5  -4  -1  -8  -2 -2 -5  6   9   0  1
## HD2_Ungated    0  -6   4  4 -2 -12  -8  17  -1  16   5 -4  5 -4  -5  -3  0
## HD3_Ungated    0   2  -1 -1  1   6   3 -12   1  -7  -3  6  1 -2  -4   3 -2
## HIV1_Ungated -23  19 -20 14 12 -31 -18   9 -10 -14  -7 -3  3 -5 -16 -10 12
## HIV2_Ungated -16   4  -4 17 -6 -18 -17  13  -6  28  -3  2  7 -5  -6  -7  7
## HIV3_Ungated -16 -11   1  6 -5 -11 -11  24 -15  29 -13  5 -6 -9 -11 -14  9
##               37  38  39 40
## HD1_Ungated    3  -1  -4  2
## HD2_Ungated   -5   0   6 -6
## HD3_Ungated    2   1  -2  3
## HIV1_Ungated -15 -16 -19 -3
## HIV2_Ungated -38 -14  -9 -8
## HIV3_Ungated -40 -11 -16 -4
cytoHeatmaps(cfData, group="group", legend=TRUE)

Box plots

# Detailed view of counts aka percentages
cytoBoxplots(cfData, group = "group")

Functional marker histograms

# Detailed of functional markers
msiPlot(cfData, markers = c("CD45", "CD4"), byGroup='group')

Do some statistics

# Add t.test
cfData <- cytottest(cfData, group="group", adjustMethod = "fdr")
head(cfData@results)
##   clusters       mean_HD  mean_HIV      pvalue   adjusted
## 1        1  0.000000e+00 3.3869024 0.001575980 0.06303922
## 2        2 -7.402390e-17 3.5342799 0.003460116 0.06920232
## 3        3  3.699840e-17 1.0399227 0.223951773 0.35832284
## 4        4  7.399680e-17 1.5870323 0.190784696 0.34688127
## 5        5 -1.849468e-17 0.6281237 0.261691958 0.40260301
## 6        6 -3.700743e-17 1.7840534 0.006007491 0.07469355
# Add some columns
# fold is a multiplicative coefficient applied to the reference
# fold is negative to stand for 1/fold
cfData@results <- within(cfData@results, {
  diff <- mean_HIV - mean_HD
  fold <- sign(diff) * round(2^abs(diff), 2)
  unselected <- abs(diff) < 0.5 | pvalue > 0.05
  label <- clusters
  label[unselected] <- NA
  })

View the Volcano plot

library(ggplot2)
## 
## Attaching package: 'ggplot2'
## The following object is masked from 'package:cytofast':
## 
##     expr
library(ggrepel)
set.seed(42)

p <- ggplot(cfData@results, aes(diff, -log10(pvalue))) + geom_point() + geom_hline(yintercept = -log10(c(0.05, 1))) + geom_vline(xintercept = c(-0.5, 0 , 0.5)) + ggtitle("Volcano Plot, diff = Log2 Fold Change of percentages")
p + geom_text_repel(aes(label = label))
## Warning: Removed 26 rows containing missing values (geom_text_repel).

p <- ggplot(cfData@results, aes(2^diff, -log10(pvalue))) + geom_point() + geom_hline(yintercept = -log10(c(0.05, 1))) + geom_vline(xintercept = c(1/1.5, 1 , 1.5)) + ggtitle("Volcano Plot, diff = Log2 Fold Change of percentages") + scale_x_continuous(trans = "log2")
p + geom_text_repel(aes(label = label))
## Warning: Removed 26 rows containing missing values (geom_text_repel).

and the table of counts

library(knitr)
kable(
  cbind(cfData@results[with(cfData@results, !unselected),], t(
    cellCountRaw[,with(cfData@results, !unselected)]))[,-c(1:3,6:7)]
)
pvalue adjusted fold diff HD1_Ungated HD2_Ungated HD3_Ungated HIV1_Ungated HIV2_Ungated HIV3_Ungated
1 0.0015760 0.0630392 10.46 3.3869024 20 5 9 324 172 151
2 0.0034601 0.0692023 11.59 3.5342799 114 318 143 3275 1407 2069
6 0.0060075 0.0746935 3.44 1.7840534 40 41 20 102 140 176
7 0.0484989 0.1385683 -8.91 -3.1550009 426 179 360 1 64 43
10 0.0112945 0.0746935 5.84 2.5460736 8 8 9 102 65 136
15 0.0148138 0.0746935 7.60 2.9259433 2 21 6 118 147 120
20 0.0149387 0.0746935 -3.60 -1.8473936 193 191 183 30 55 55
23 0.0470964 0.1385683 2.33 1.2222525 93 151 109 312 404 178
26 0.0415599 0.1385683 -2.92 -1.5464926 102 35 84 12 14 26
35 0.0199767 0.0887852 -2.06 -1.0447995 1195 966 1448 583 730 435
36 0.0123828 0.0746935 1.87 0.9053814 42 38 32 97 64 77
37 0.0478986 0.1385683 -8.70 -3.1215816 1965 1118 1889 548 105 90
38 0.0088101 0.0746935 -2.60 -1.3802123 1207 1242 1318 397 471 576
39 0.0245895 0.0983579 -2.78 -1.4742611 66 146 83 18 44 24