cytofCore.updateFlowFrameKeywords = function(flowFrame){
params = parameters(flowFrame)
pdata = pData(params)
for (i in 1:ncol(flowFrame)){
s = paste("$P",i,"S",sep="");
n = paste("$P",i,"N",sep="");
r = paste("$P",i,"R",sep="");
b = paste("$P",i,"B",sep="");
e = paste("$P",i,"E",sep="");
keyval=list();
keyval[[s]] = colnames(flowFrame)[i];
keyval[[n]] = colnames(flowFrame)[i];
keyval[[r]] = ceiling(max(exprs(flowFrame)[,i])-min(exprs(flowFrame)[,i]))
keyval[[b]] = 32;
keyval[[e]] = "0,0";
keyword(flowFrame) = keyval;
pdata[i,"minRange"]=min(exprs(flowFrame)[,i])
pdata[i,"maxRange"]=max(exprs(flowFrame)[,i])
}
pData(params)=pdata
parameters(flowFrame)=params
# keyval[["$DATATYPE"]] <- "F"
return(flowFrame)
}
cytofCore.combineChannels = function(flowFrame,channelList,newName=NULL){
if (!all(channelList %in% colnames(flowFrame))){
offending = which(!(channelList %in% colnames(flowFrame)))
stop( paste("Channel(s)",paste(channelList[offending],collapse=","),"not found in frame.") )
}
# Add data in combined channels
combinedValues = apply(exprs(flowFrame)[,channelList],1,sum)
channelNumber = ncol(flowFrame)+1
channelId = paste("$P",channelNumber,sep="");
if (is.null(newName)){
newName = paste("channel",paste(channelList,collapse="-"),"combined",sep="_");
}
params = parameters(flowFrame)
pdata = pData(params)
plist = matrix(c(newName,newName,ceiling(max(combinedValues)-min(combinedValues)),floor(min(combinedValues)),ceiling(max(combinedValues))));
rownames(plist) = c("name","desc","range","minRange","maxRange");
colnames(plist) = c(channelId);
pdata = rbind(pdata,t(plist))
pData(params) = pdata
return(cytofCore.updateFlowFrameKeywords(flowFrame(exprs=cbind(exprs(flowFrame), `colnames<-`(cbind(combinedValues), newName)),parameters=params,description=description(flowFrame))))
}
cytofCore.subtract = function(flowFrame,value=100,exclude=c(1,2)){
if (length(exclude)>0){
if (!is.numeric(exclude)){
if (!all(exclude %in% colnames(flowFrame))){
offending = which(!(exclude %in% colnames(flowFrame)))
stop( paste("Channel(s)",paste(exclude[offending],collapse=","),"not found in frame.") )
}
} else {
exclude = colnames(flowFrame)[exclude];
}
}
subtractCols = setdiff(colnames(flowFrame),exclude)
exprs(flowFrame)[,subtractCols] = exprs(flowFrame)[,subtractCols]-value
return(cytofCore.updateFlowFrameKeywords(flowFrame))
}
cytofCore.concatenateDirectoryFiles = function(inputDir,outputDir=NULL,pattern=NULL,overwrite=F,timeParam="time"){
currentwd = getwd();
if (is.null(outputDir)){
outputDir = inputDir;
}
if (!file.exists(inputDir)){
stop(paste("Input Directory",inputDir,"not found."))
}
if (!file.exists(outputDir)){
stop(paste("Output Directory",outputDir,"not found."))
}
setwd(inputDir);
if (is.null(pattern)){
outputFileName="combined.fcs";
} else {
outputFileName = paste(paste(pattern,collapse="-"),"_combined.fcs",sep="")
}
if (file.exists(outputFileName) && !overwrite){
stop(paste("File ",outputFileName," already exists in dir ",outputDir,". Set overwrite argument to TRUE or remove file first.",sep=""))
}
fcsFiles = list.files(inputDir, pattern=".fcs$", ignore.case=TRUE)
matched = c(1:length(fcsFiles));
for (match in pattern){
matched=intersect(matched,grep(match,fcsFiles))
}
combineFiles = fcsFiles[matched]
begin=T
for (file in combineFiles){
cat(paste("Adding",file,"\n"));
if (begin){
fcs = read.FCS(file)
combinedData = exprs(fcs);
columns = colnames(combinedData)
timeCol = which(tolower(columns) == "time")
if (!length(timeCol)){
stop(paste("Time Column",timeParam,"not found for file",file));
}
fcsDescription = description(fcs)
channels = parameters(fcs)$desc
begin=F
} else {
tmpData = exprs(read.FCS(file));
# check that column names are the same before concatenating
if (any(colnames(tmpData) != columns)) {
stop(paste("File", file, "does not have matching panel."));
}
tmpData[,timeCol] = tmpData[,timeCol]+(max(combinedData[,timeCol])+round((max(combinedData[,timeCol])-min(combinedData[,timeCol]))/nrow(combinedData)))
combinedData = rbind(combinedData,tmpData)
}
}
# copy the channel descriptions from the first file
cytofCore.write.FCS(combinedData, outputFileName, what="numeric", channelDescription=channels,
referenceDescription=fcsDescription, oldDescription=fcsDescription)
gc()
setwd(currentwd);
return(paste(outputDir,"->",outputFileName))
}
cytofCore.concatenateFiles = function(fcsFiles,timeCol="Time"){
begin=T
for (file in fcsFiles){
if (!file.exists(file)){
stop(paste(file,"Not Found! Stopping."))
}
cat(paste("Adding",file,"\n"));
if (begin){
combinedData = exprs(read.FCS(file));
if (!(timeCol %in% colnames(combinedData))){
stop(paste("Time Column",timeCol,"not found for file",file));
}
begin=F
} else {
tmpData = exprs(read.FCS(file));
if (!all(colnames(tmpData)==colnames(combinedData))){
stop(paste("File names in files do not match.",colnames(tmpData),colnames(combinedData),sep="\n"))
}
if (!(timeCol %in% colnames(tmpData))){
stop(paste("Time Column",timeCol,"not found for file",file));
}
tmpData[,timeCol] = tmpData[,timeCol]+(max(combinedData[,timeCol])+round((max(combinedData[,timeCol])-min(combinedData[,timeCol]))/nrow(combinedData)))
combinedData = rbind(combinedData,tmpData)
}
}
outputFileName = paste(dirname(fcsFiles[1]),"/concatenated_",gsub("\\s+","_",date()),".fcs",sep="")
write.FCS(cytofCore.updateFlowFrameKeywords(flowFrame(exprs=combinedData)),outputFileName);
return(outputFileName)
}
cytofCore.updatePanel = function(templateFile=NULL, fcsFolder=NULL){
if (is.null(templateFile) || is.null(fcsFolder)) {
library("tcltk")
}
if (is.null(templateFile)) {
#select template file, which has list of metals you want to include and their desired marker names
Filters <- matrix(c("template", ".txt","template", ".TXT", "template",".fcs","template",".FCS","template",".csv","template","CSV"), 6, 2, byrow = TRUE)
templateFile<- tk_choose.files(caption="Select template file", multi=FALSE, filters=Filters)
}
if (!file.exists(templateFile)) {
stop(paste("Template file", templateFile, "not found."))
}
templateExt = substr(templateFile,nchar(templateFile)-2,nchar(templateFile))
#read in metal and marker list from template file
if (templateExt=="csv" || templateExt=="CSV") {
data =read.csv(templateFile, header=FALSE)
channels=as.character(data[,1])
markers=as.character(data[,2])
} else if (templateExt=="txt" || templateExt=="TXT") {
data=read.table(templateFile, header=FALSE)
channels=as.character(data[,1])
markers=as.character(data[,2])
} else if (templateExt=="fcs" || templateExt=="FCS") {
data=read.FCS(templateFile)
channels=as.character(parameters(data)$name)
markers=as.character(parameters(data)$desc)
}
#for simplicity, copy over the channel values when the markers are empty
for (i in 1:length(markers)) {
if ((nchar(markers[i]) == 0) || is.na(markers[i])) {
markers[i]=channels[i]
}
}
#select folder of FCS files to relabel
if (is.null(fcsFolder)) {
fcsFolder<-tk_choose.dir()
}
newFolder<-file.path(fcsFolder,"relabeled")
dir.create(newFolder)
fcsFiles<-list.files(path=fcsFolder, pattern=".fcs$", ignore.case = TRUE)
addZeroCol<-matrix(0,ncol=length(channels)) #keep track of whether to add missing column for all files
for (file in fcsFiles) {
oldFile=read.FCS(file.path(fcsFolder,file))
oldChannels=as.character(parameters(oldFile)$name)
oldData=exprs(oldFile)
newData=matrix(data=0,nrow=nrow(oldData),ncol=length(channels))
keepCols=matrix(TRUE,ncol=length(channels)) #keep track of columns to include in final matrix
for (i in 1:length(channels)) {
ind=grep(paste("\\Q",channels[i],"\\E",sep=""),oldChannels)
if (length(ind)==1) { #found this template metal in this fcs file
newData[,i]=oldData[,ind]
} else if (addZeroCol[i]==0) { #haven't decided yet whether to create dummy column
print(paste("Channel", channels[i], "was not found in",file,"!"))
yn=readline("Do you want to create a column of zeros in its place? y/n: ")
allone=readline("Do you want to apply this to all FCS files in this folder? y/n: ")
if (yn=="y" && allone=="y") { #create dummy column in this and all files
addZeroCol[i]=1
} else if (yn=="n" && allone=="y") { #don't create dummy column in any file
addZeroCol[i]=-1
keepCols[i]=FALSE
} else if (yn=="n") {#don't create dummy column but only in this file
keepCols[i]=FALSE
}
} else if (addZeroCol[i]==-1) { #previously decided not to create this dummy column
keepCols[i]=FALSE
}
}
newData=newData[,keepCols]
newChannels=channels[keepCols]
newMarkers=markers[keepCols]
colnames(newData) = newChannels
print(paste("Writing",file))
referenceDescription = NULL
if (templateExt=="fcs" || templateExt=="FCS") {
referenceDescription=description(data)
}
cytofCore.write.FCS(newData,file.path(newFolder,file),channelDescriptions=newMarkers,referenceDescription=referenceDescription,oldDescription=description(oldFile))
}
}
cytofCore.averageDdCoefficients = function(folder,masses) {
#average slopes and intercepts across a set of conf files for selected masses
fileList=list.files(path=folder, pattern=".conf$")
slopes=matrix(ncol=length(masses),nrow=length(fileList))
intercepts=matrix(ncol=length(masses),nrow=length(fileList))
for (i in 1:length(fileList)) {
conf=cytofCore.read.conf(file.path(folder,fileList[i]))
if (i==1) {
massesMeasured=conf$Mass
} else if (any(conf$Mass != massesMeasured)) {
stop("Panels of conf files are not consistent.")
}
slopes[i,]=approx(conf$Mass,conf$Slope,masses,method="linear",rule=2)$y
intercepts[i,]=approx(conf$Mass,conf$Intercept,masses,method="linear",rule=2)$y
}
results=rbind(colMeans(slopes),colMeans(intercepts))
colnames(results)=masses
rownames(results)=c("slopes","intercepts")
return(results)
}
cytofCore.rewriteImdCoeffs = function(imdFile,confFolder) {
#rewrite the slopes and intercepts in the dual calibration as the average from a set of conf files
library("XML")
# extract xml from tail of imd
imd <- file(imdFile, "r+b")
if (!isSeekable(imd)) {
stop("Cannot seek in specified file or connection")
}
endTag="</ExperimentSchema>"
seek(imd,where=-2*nchar(endTag),origin="end")
fileEndString=paste(readBin(imd, what="character", size=2, n=2*nchar(endTag), signed=FALSE),collapse="")
if (fileEndString != endTag) {
stop("The xml tail is either missing or irregular.")
}
#read in chunks until the beginning of the xml is reached
chunkLength=1024
seek(imd,where=-2*chunkLength,origin="current")
xmlChunk=c()
while (!(0 %in% xmlChunk)) {
xmlChunk=c(readBin(imd, what="integer", size=2, n=chunkLength, signed=FALSE),xmlChunk)
seek(imd,where=-4*chunkLength,origin="current")
if (length(xmlChunk)>2^17) {stop("Didn't find zeros preceding xml.")}
}
# return to the position at which the xmlChunk begins
seek(imd,where=2*chunkLength,origin="current")
# convert to char and trim extra before the xml
startTag="<ExperimentSchema"
matchInd=match(utf8ToInt("<"),xmlChunk)
if (is.na(matchInd)) {stop("XML start tag not found")}
imdString=sub(".*<ExperimentSchema","<ExperimentSchema",intToUtf8(xmlChunk[matchInd:length(xmlChunk)]))
if (substr(imdString,1,nchar(startTag)) != startTag) {stop("XML start not aligned.")}
# parse xml
xmlList=xmlToList(xmlInternalTreeParse(imdString))
# dual calibration info
dualList=xmlList[names(xmlList)=="DualAnalytesSnapshot"]
masses=c()
oldSlopes=c()
oldIntercepts=c()
for (metal in dualList) {
masses=c(masses,metal$Mass)
oldSlopes=c(oldSlopes,metal$DualSlope)
oldIntercepts=c(oldIntercepts,metal$DualIntercept)
#dualCalibration=rbind(dualCalibration,c(metal$Mass,metal$DualSlope, metal$DualIntercept))
}
#colnames(dualCalibration)=c("mass","slope","intercept")
#get new coeffs
averagedCoeffs=cytofCore.averageDdCoefficients(confFolder,as.numeric(masses))
newXml=imdString
for (i in 1:length(masses)) {
#make new slope string same length as old slope string
numChars=nchar(oldSlopes[i])
slopeString=as.character(averagedCoeffs["slopes",i])
charDiff=numChars-nchar(slopeString)
if (charDiff>0) {
slopeString=paste(slopeString,paste(rep("0",charDiff),sep="",collapse=""),sep="")
} else if (charDiff<0) { slopeString=substr(slopeString,1,numChars) }
newXml=sub(oldSlopes[i],slopeString,newXml,fixed=T)
#make new intercept string same length as old intercept string
numChars=nchar(oldIntercepts[i])
intString=as.character(averagedCoeffs["intercepts",i])
charDiff=numChars-nchar(intString)
if (charDiff>0) {
intString=paste(intString,paste(rep("0",charDiff),sep="",collapse=""),sep="")
} else if (charDiff<0) {
intString=substr(intString,1,numChars)
}
newXml=sub(oldIntercepts[i],intString,newXml,fixed=T)
}
# convert back to int and make new xmlChunk
newXmlInt=utf8ToInt(newXml)
if (length(newXmlInt) != (length(xmlChunk)-matchInd + 1)) {
stop("New xml is not same length as old xml.")
}
xmlChunk[matchInd:length(xmlChunk)]=newXmlInt
# overwrite the old xml with the new xml, given that the IMD file is still open at the same location
writeBin(xmlChunk,imd,size=2)
close(imd)
cat(imdFile, " has been rewritten.")
}
cytofCore.copyImdXml = function(sourceImd,targetImd) {
# Write an XML Tail on an IMD file that was truncated before the tail was written by copying it from another IMD file.
# Note that the other IMD file MUST have the same panel and dual slopes for this to be accurate.
# If you do not have another IMD file with the same panel and tail, create one with this panel.
# If the dual slopes have been adjusted since the target file was acquired, it is possible to acquire in Dd mode and then rewrite the slopes with cytofCore.rewriteImdCoeffs
# read XML from source IMD file and convert to integers
sourceXml=cytofCore.read.imd.xml(sourceImd)
imdInt=utf8ToInt(sourceXml$rawText)
# write XML to target file
imd <- file(targetImd, "ab")
if (!isSeekable(imd)) {
stop("Cannot seek in specified file or connection")
}
writeBin(imdInt,imd,size=2)
close(imd)
cat("XML from ",sourceImd, " has been appended to ",targetImd)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.