### R code from vignette source 'vignettes/RchyOptimyx/inst/doc/RchyOptimyx.Rnw'

###################################################
### code chunk number 1: c10
###################################################
library(flowType)
data(HIVData)
data(HIVMetaData)
HIVMetaData <- HIVMetaData[which(HIVMetaData[,'Tube']==2),];


###################################################
### code chunk number 2: c11
###################################################
Labels=(HIVMetaData[,2]=='+')+1;


###################################################
### code chunk number 3: c12
###################################################
library(flowCore);
library(RchyOptimyx);

##Markers for which cell proportions will be measured.
PropMarkers <- 5:10

##Markers for which MFIs will be measured.
MFIMarkers <- PropMarkers

##Marker Names
MarkerNames <- c('Time', 'FSC-A','FSC-H','SSC-A',
                 'IgG','CD38','CD19','CD3',
                 'CD27','CD20', 'NA', 'NA')

##Apply flowType
ResList <- fsApply(HIVData, 'flowType', PropMarkers, 
                   MFIMarkers, 'kmeans', MarkerNames);

##Extract phenotype names
phenotype.names=names(ResList[[1]]@CellFreqs)


###################################################
### code chunk number 4: c13
###################################################
all.proportions <- matrix(0,3^length(PropMarkers)
                          -1,length(HIVMetaData[,1]))
for (i in 1:length(ResList)){
  all.proportions[,i] = ResList[[i]]@CellFreqs
  all.proportions[,i] = all.proportions[,i] / 
    max(all.proportions
        [which(names(ResList[[i]]@CellFreqs)==''),i])
}


###################################################
### code chunk number 5: c14
###################################################
Pvals <- vector();
EffectSize <- vector();
for (i in 1:dim(all.proportions)[1]){

  ##If all of the cell proportions are 1 (i.e., the phenotype 
  ##with no gates) the p-value is 1.
  if (length(which(all.proportions[i,]!=1))==0){
    Pvals[i]=1;
    EffectSize[i]=0;
    next;
  }
  temp=t.test(all.proportions[i, Labels==1], 
    all.proportions[i, Labels==2])
  Pvals[i] <- temp$p.value
  EffectSize[i] <- abs(temp$statistic)  
}

Pvals[is.nan(Pvals)]=1
names(Pvals)=phenotype.names

##Bonferroni's correction
selected <- which(p.adjust(Pvals)<0.1);

print(names(selected))


###################################################
### code chunk number 6: c15
###################################################
library(sfsmisc)
Signs=t(digitsBase(1:(3^length(PropMarkers)-1),
  3,ndigits=length(PropMarkers)))
rownames(Signs)=phenotype.names;
colnames(Signs)=MarkerNames[PropMarkers]
head(Signs)


###################################################
### code chunk number 7: c16
###################################################
res<-RchyOptimyx(Signs, -log10(Pvals), 
                 paste(Signs['IgG-CD38-CD19-CD27+CD20-',], 
                       collapse=''), factorial(6), FALSE)
plot(res, phenotypeScores=-log10(Pvals), ylab='-log10(Pvalue)')


###################################################
### code chunk number 8: c17
###################################################
res<-RchyOptimyx(Signs, -log10(Pvals), 
                 paste(Signs['IgG-CD38-CD19-CD27+CD20-',], 
                       collapse=''), 15, FALSE)
plot(res, phenotypeScores=-log10(Pvals), ylab='-log10(Pvalue)')


###################################################
### code chunk number 9: loadlibs
###################################################
library(RchyOptimyx)
data(HIVData)


###################################################
### code chunk number 10: c0
###################################################
LogRankPvals['KI-67+CD4-CCR5+CD127-']


###################################################
### code chunk number 11: c01
###################################################
paste(Signs['KI-67+CD4-CCR5+CD127-',], collapse='')


###################################################
### code chunk number 12: c1
###################################################
par(mar=c(1,4,1,1))
res<-RchyOptimyx(Signs, -log10(LogRankPvals), 
                 '2111012110', 24,FALSE)
plot(res, phenotypeScores=-log10(LogRankPvals), 
     ylab='-log10(Pvalue)')


###################################################
### code chunk number 13: c2
###################################################
par(mar=c(1,4,1,1))
plot(res, phenotypeScores=-log10(LogRankPvals), ylab='-log10(Pvalue)', 
     uniformColors=TRUE, edgeWeights=FALSE, edgeLabels=FALSE, 
     nodeLabels=TRUE)


###################################################
### code chunk number 14: c3
###################################################
summary(res)


###################################################
### code chunk number 15: c4
###################################################
plot(ecdf(res@pathScores), verticals=TRUE)


###################################################
### code chunk number 16: c5
###################################################
par(mar=c(1,4,1,1))
res<-RchyOptimyx(Signs, -log10(LogRankPvals), '2111012110', 4,FALSE)
plot(res, phenotypeScores=-log10(LogRankPvals), ylab='-log10(Pvalue)')


###################################################
### code chunk number 17: c6
###################################################
LogRankPvals['KI-67+CCR5+']


###################################################
### code chunk number 18: c7
###################################################
res<-RchyOptimyx(Signs, OverlapScores, 
                 paste(Signs['CD28+CD45RO-CD57-CCR5-CD27+CCR7+',], 
                       collapse=''), 720, FALSE)
par(mar=c(1,4,1,1))
plot(res, phenotypeScores=OverlapScores, ylab='Purity',
     uniformColors=TRUE, edgeWeights=FALSE, edgeLabels=FALSE,
     nodeLabels=TRUE)


###################################################
### code chunk number 19: c8
###################################################
plot(ecdf(res@pathScores), verticals=TRUE)


###################################################
### code chunk number 20: c9
###################################################
res<-RchyOptimyx(Signs, OverlapScores, 
                 paste(Signs['CD28+CD45RO-CD57-CCR5-CD27+CCR7+',],
                       collapse=''), 5, FALSE)
par(mar=c(1,4,1,1))
plot(res, phenotypeScores=OverlapScores, ylab='Purity')


###################################################
### code chunk number 21: c10
###################################################
OverlapScores['CD45RO-CCR5-CCR7+']