View on GitHub

OpenCyto


A Robust BioConductor Framework for Automated Flow Data Analysis

Analysis of SDY207 Automated Gating Results

Analysis of SDY207 Automated Gating Results

We load some required libraries an define a few functions that we'll need later.

require(openCyto)
require(pheatmap)
require(COMPASS)
require(Kmisc)
require(reshape2)
require(reshape)
require(plyr)
require(ggplot2)
require(scales)
OUTPATH<-"./" #output path
WORKDIR<-"/Users/gfinak/Dropbox/GoTeam/Projects/Paper-OpenCytoPipeline/SDY207/reproducible/" #working directory, should be customized.
setwd(WORKDIR)

#hyperbolic arcsine transform for ggplot2
asinh_trans<-function(c){
  trans_new(name = 'asinh', transform = function(x) asinh(x*c), 
            inverse = function(x) sinh(x)/c)
}

#Parses boolean gate names
parse_name <- function(s) {
    split <- unlist(strsplit(sub('^\\d+\\.','',s),'&'))
    n = c(1,1,1,1,1,1,1,1,1)
    n[grep('\\!',split)] <- 0
    return(n)
}

#formats the population  statistics into a matrix for a heatmap
.reformat<-function(x){
  ss<-list()
  for(j in 1:ncol(x)){
    s<-NULL
    for(i in c("TN","TEF","TEM","TCM")){
      ii<-grep(paste0(i,"$"),rownames(x))
      n_cyto<-9
      s<-cbind(s,x[ii:(ii+2**n_cyto),j])
    }
    colnames(s)<-c("TN","TEF","TEM","TCM")
    names<-strsplit(sub('^\\d+\\.','',do.call(c,lapply(strsplit(rownames(s),"/"),function(x)x[length(x)]))),"&")
    s<-cbind(s,t(apply(do.call(rbind,names),1,parse_name)))
    colnames(s)<-c(colnames(s)[1:4],gsub("\\!","",names[[2]]))
    ss<-c(ss,list(s))
  }
  ss
}

Next we load the OpenCyto gated data that we generated here

gs <- load_gs('./gating/open_cyt_gs.dat')
## loading R object...
## loading tree object...
## Done

#extract stats
foo<-getPopStats(gs)

# format into matrices
ss<-.reformat(foo)
level <- (0.01)

#some filtering for populations that are expressed at least 1% of the time across any subset in any sample
idxs<-do.call(`|`,lapply(ss,function(x)rowSums(x[,1:4]>level)>0))
idxs[1] <- FALSE

# pull out information for annotating the heatmap.
a <- data.frame(ss[[1]][idxs,c(5:11)])
a$TNFa <- factor(a$TNFa == 1, labels=c('-','+'))
a$MIP1a <- factor(a$MIP1a == 1, labels=c('-','+')) 
a$IL2 <- factor(a$IL2 == 1, labels=c('-','+'))
a$CD107 <- factor(a$CD107==1, labels=c('-','+'))
a$IFNgamma <- factor(a$IFNgamma == 1, labels=c('-','+'))
a$MIP1b <- factor(a$MIP1b == 1, labels=c('-','+'))
a$GMCSF <- factor(a$GMCSF == 1, labels=c('-','+'))
cytokine_ann<-data.frame(matrix(as.integer(a=="+"),ncol=ncol(a),nrow=nrow(a)))
cytokine_ann<-as.data.frame(lapply(cytokine_ann,factor))
rownames(cytokine_ann)<-rownames(a)
colnames(cytokine_ann)<-colnames(a)

# We'll plot the average expression for each subset across samples
d<-do.call('+',lapply(ss,function(x)x[idxs,1:4]))/2
breaks <- unique(quantile(asinh(d*1000),seq(0,1,l=101)))
colors=colorRampPalette(RColorBrewer::brewer.pal(9,name="Blues"))(length(breaks)-1)

# These are the scale_labels
scale.labels<-sinh(1:7)/1000

Heatmap

A heatmap of the average proportion of cells (relative to the parent T-cell maturational subset, TN, TEF, TCM, or TEM) expressing each of the cytokine combinations in the legend at least at the 1% level in any sample.

COMPASS:::pheatmap(asinh(t(d*1000)), cluster_cols=F, cluster_rows=T, breaks=breaks, show_colnames=F,color=colors,cytokine_annotation=cytokine_ann,legend_breaks=1:7,legend_labels=signif(scale.labels,2))

plot of chunk unnamed-chunk-3

## pdf 
##   2
## pdf 
##   2

Degree of functionality

The distribution of cells within each T-cell maturational subset with different degrees of functionality. Data are rescaled so that each T-cell subset integrates to 1.

degree<-apply(cytokine_ann,1,function(x)sum(x==1))
ggplot(melt(t(apply(do.call(rbind,by((d),degree,function(x)colMeans(x))),2,function(x)x/sum(x)))))+geom_line(aes(x=X2,y=value,group=X1,lty=X1,col=X1),lwd=1)+scale_x_discrete("Degree of Functionality")+theme_bw()+scale_y_continuous(wrap("Fraction of Cells",25))+scale_linetype_discrete("Maturational Subset")+scale_color_discrete("Maturational Subset")+coord_trans(ytrans=GFMisc::asinh_trans(3))+theme(legend.position=c(0.8,0.63))

plot of chunk unnamed-chunk-4

## pdf 
##   2
## pdf 
##   2