## Define a function that calculates the NRS per sample

NRS <- function(x, ncomp = 3){
  pr <- prcomp(x, center = TRUE, scale. = FALSE) 
  score <- rowSums(outer(rep(1, ncol(x)), 
    pr$sdev[1:ncomp]^2) * abs(pr$rotation[,1:ncomp]))
  return(score)
}

nrs_sample <- fsApply(fcs[, lineage_markers], NRS, use.exprs = TRUE)

nrs <- colMeans(nrs_sample, na.rm = TRUE)

lineage_markers_ord <- names(sort(nrs, decreasing = TRUE))

ggdf <- melt(nrs_sample, id.var = "sample_id", 
  value.name = "nrs", variable.name = "antigen")

ggdf$antigen <- factor(ggdf$antigen, levels = lineage_markers_ord)

ggplot(ggdf, aes(x = antigen, y = nrs)) +
  geom_point(aes(color = condition), alpha = 0.9, 
    position = position_jitter(width = 0.3, height = 0)) +
  geom_boxplot(outlier.color = NA, fill = NA) +
  stat_summary(fun.y = "mean", geom = "point", shape = 21, fill = "white") +
  theme_bw() +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  scale_color_manual(values = color_conditions)