Calculating genetic interaction strengths with qfa R package.

GISFunctions.R

f(!("qfa" %in% rownames(installed.packages()))) install.packages("qfa", repos="http://R-Forge.R-project.org")
library(qfa)

GISReport=function(qfolder,cfolder,qTrtMed,cTrtMed,colName="MDRMDP",GISdirn=NULL,normalised=FALSE,wctest=FALSE,qStripRep=0,cStripRep=0,bootstrap=NULL,froot="QFA_EXPERIMENTS",makeImBrowser=FALSE,imTime=2.0,hpath=NULL){
        cat("###### WORKING DIRECTORY MUST BE LOGS(2) ######\n")
        cat("Query:",qfolder,qTrtMed,"Control:",cfolder,cTrtMed,"Directory:",GISdirn,"\n")

        #################################################################
        if(normalised){
                quer=file.path(froot,qfolder,"ANALYSISOUT",paste(qfolder,"NORM_FIT.txt",sep="_"))
                con=file.path(froot,cfolder,"ANALYSISOUT",paste(cfolder,"NORM_FIT.txt",sep="_"))
        }else{
                quer=file.path(froot,qfolder,"ANALYSISOUT",paste(qfolder,"FIT.txt",sep="_"))
                con=file.path(froot,cfolder,"ANALYSISOUT",paste(cfolder,"FIT.txt",sep="_"))
        }

        quer.fit=do.call("rbind",lapply(quer, read.delim,header=TRUE,sep="\t",stringsAsFactors=FALSE))
        control.fit=do.call("rbind",lapply(con, read.delim,header=TRUE,sep="\t",stringsAsFactors=FALSE))

        quer.fit$fit=quer.fit[[colName]]
        control.fit$fit=control.fit[[colName]]

        quer.fit=quer.fit[!quer.fit$RepQuad%in%qStripRep,]
        control.fit=control.fit[!control.fit$RepQuad%in%cStripRep,]

        # Add a TrtMed column, but only if it is missing...
        if(!"TrtMed"%in%colnames(quer.fit)) quer.fit$TrtMed=paste(quer.fit$Treatment,quer.fit$Medium,sep="_")
        if(!"TrtMed"%in%colnames(control.fit)) control.fit$TrtMed=paste(control.fit$Treatment,control.fit$Medium,sep="_")

        if(is.null(GISdirn)) GISdirn=file.path(froot,sort(qfolder)[1],"ANALYSISOUT")
        if(is.null(hpath)) GISdirn=file.path(froot,sort(qfolder)[1],"IMAGEBROWSER")

        # tell R where to store the comparison results
        QUER = quer.fit[quer.fit$ScreenID==qfolder,]
        CONT = control.fit[control.fit$ScreenID==cfolder,]

        a = QUER[QUER$TrtMed==qTrtMed,]
        b = CONT[CONT$TrtMed==cTrtMed,]
        clab = paste(unique(b$ScreenID),unique(b$Inoc),unique(b$Library),unique(b$User),unique(b$Screen.Name),unique(b$ExptDate),unique(b$TrtMed),collapse=" ")
        qlab = paste(unique(a$ScreenID),unique(a$Inoc),unique(a$Library),unique(a$User),unique(a$Screen.Name),unique(a$ExptDate),unique(a$TrtMed),collapse=" ")
        root = paste(unique(QUER$Client),qfolder[1],unique(QUER$Screen.Name),qTrtMed,"vs",cfolder[1],unique(CONT$Screen.Name),cTrtMed,sep="_")

        pdf(file.path(froot,sort(qfolder)[1],"ANALYSISOUT",paste(root,"_",colName[1],"_t-test.pdf",sep="")))

                # Calculate genetic interactions and produce epistasis plot
                epi=qfa.epi(a,b,0.05,plot=FALSE,wctest=wctest)
                if(wctest) mmain=paste("Normalised=",normalised[1]," ",colName[1],"Fitness plot (Wilcoxon test)")
                if(!wctest) mmain=paste("Normalised=",normalised[1]," ",colName[1],"Fitness plot (t-test)")
                if(!is.null(bootstrap)) mmain=paste("Normalised=",normalised[1]," ",colName[1],"Fitness plot (bootstrap)")

                qfa.epiplot(epi$Results,0.05,xxlab=clab,yylab=qlab,mmain=mmain,fmax=0)

                report.epi(epi$Results,file=file.path(GISdirn,paste(root,"_",colName[1],"_GIS.txt",sep="")))
                report.epi(epi$Results,file=file.path(froot,sort(qfolder)[1],"ANALYSISOUT",paste(root,"_",colName[1],"_GIS.txt",sep="")))

        dev.off()

        if(makeImBrowser){
                print("Making image browser...")
                ImBrowser(qfolder,cfolder,qTrtMed,cTrtMed,froot,imTime,hpath)
        }
}

ImBrowser=function(expt1,expt2,trtMed1,trtMed2,froot="QFA_EXPERIMENTS",imtime=2.0,hpath=NULL){
        trtMed1=gsub(" ","#~#",trtMed1)
        trtMed2=gsub(" ","#~#",trtMed2)
        string=paste("python /home/yeastimages/LOGS3/Scripts/makeHTML_GIS.py",expt1,expt2,trtMed1,trtMed2,froot,as.character(imtime),hpath)
        print(string)
        system(string)
}

getTreatMed=function(folder,normalised=TRUE){
        cat("###### WORKING DIRECTORY MUST BE LOGS(2) ######\n")
        res=list()
        for(qfolder in folder){
                if(normalised){
                        quer=paste("QFA_EXPERIMENTS/",qfolder,"/ANALYSISOUT/",qfolder,"_NORM_FIT.txt",sep="")
                }else{
                        quer=paste("QFA_EXPERIMENTS/",qfolder,"/ANALYSISOUT/",qfolder,"_FIT.txt",sep="")
                }
                quer.fit=read.delim(quer,header=TRUE,sep="\t",stringsAsFactors=FALSE)
                # Add a TrtMed column, but only if it is missing...
                if(!"TrtMed"%in%colnames(quer.fit)) quer.fit$TrtMed=paste(quer.fit$Treatment,quer.fit$Medium)
                res[[qfolder]]=list(Treatment=unique(quer.fit$Treatment), Medium=unique(quer.fit$Medium),TrtMed=unique(quer.fit$TrtMed))
        }
        res$nAUC=names(quer.fit)[grepl("nAUC.*",names(quer.fit))]
        res$nSTP=names(quer.fit)[grepl("nSTP.*",names(quer.fit))]
        return(res)
}
~
~
"Scripts/GISFunctions.R" [dos] 89L, 4602C