## ----DmelSGIstyle, echo=FALSE, results="asis"------------------------------ BiocStyle::latex(bibstyle="apalike") ## ----DmelSGIopts,include=FALSE--------------------------------------------- library(knitr) opts_chunk$set(concordance=TRUE, resize.width="0.5\\textwidth", dev=c('pdf','png'), cache=TRUE, tidy = FALSE) ## ----DmelSGIinstallation,eval=FALSE---------------------------------------- # if (!requireNamespace("BiocManager", quietly=TRUE)) # install.packages("BiocManager") # BiocManager::install("DmelSGI") ## ----DmelSGIfunctions------------------------------------------------------ fpath <- function(d) { file.path(opts_knit$get("output.dir"), "result", d,"") } ## ----dataAcceess1,results='hide'------------------------------------------- library("DmelSGI") data("Interactions", package="DmelSGI") ## ----dataAccess2,eval=FALSE------------------------------------------------ # ? Interactions ## ----ReanalysisOfHornEtALLoadLibrary, results='hide', message=FALSE-------- library("DmelSGI") library("RNAinteractMAPK") basedir = getBaseDir() resultdir = file.path( basedir, "result", "ReanalysisOfHornEtAl") dir.create(resultdir, recursive = TRUE,showWarnings=FALSE) ## ----ReanalysisOfHornEtALLoadData------------------------------------------ data("Dmel2PPMAPK", package="RNAinteractMAPK") print(Dmel2PPMAPK) PI <- getData(Dmel2PPMAPK, type="pi", format="targetMatrix", screen="mean", withoutgroups = c("pos", "neg"))[,,1,] ## ----ReanalysisOfHornEtALNormalize----------------------------------------- for (j in 1:dim(PI)[2]) { for (k in 1:dim(PI)[3]) { PI[,j,k] = PI[,j,k] / (sqrt(sum(PI[,j,k] * PI[,j,k]) / (dim(PI)[2]-1))) } } ## ----ReanalysisOfHornEtALMaxVarianceSelection,eval=FALSE------------------- # Selected = c() # Selected = c(1,2,3) # R = 1:dim(PI)[1] # Res = PI # openVar = rep(-1,dim(PI)[1]+1) # openVar[1] = sum(Res * Res) / (dim(PI)[1]*(dim(PI)[2]-1)*dim(PI)[3]) # for (i in 1:dim(PI)[1]) { # H = rep(100000000.0,length(R)) # for (j in 1:length(R)) { # cat("i=",i," j=",j,"\n") # k=1:3 # A = PI[,c(Selected[seq_len(i-1)],R[j]),k,drop=FALSE] # dim(A) = c(dim(A)[1],prod(dim(A)[2:3])) # B = PI[,-c(Selected[seq_len(i-1)],R[j]),k,drop=FALSE] # dim(B) = c(dim(B)[1],prod(dim(B)[2:3])) # # Res = matrix(0.0, nr=dim(PI)[1],nc=ncol(B)) # for (z in 1:ncol(B)) { # model = lm(B[,z]~A+0) # Res[,z] = model$residuals # } # H[j] = sum(Res * Res) / (dim(PI)[1]*(dim(PI)[2]-1)*dim(PI)[3]) # } # M = which.min(H) # cat("selected: ",R[M],"\n") # Selected = c(Selected, R[M]) # openVar[i+1] = H[M] # R = R[-M] # } ## ----ReanalysisOfHornEtAL-------------------------------------------------- openVar = c(1, 0.584295886914632, 0.49354448724904, 0.440095163032832, 0.37969110256306, 0.330693818887106, 0.28896777328302, 0.26144276377077, 0.24550380797587, 0.212282252772014, 0.19041097617251, 0.16974901306481, 0.15642204582756, 0.141467140253324, 0.12781027389229, 0.11609596000734, 0.10374891651534, 0.093268306952119, 0.08446425055463, 0.07404659630757, 0.06599890651265, 0.057244319680828, 0.04944008500553, 0.04161924747819, 0.03515950952616, 0.028667487889006, 0.02313772533424, 0.01727915218118, 0.01282727545013, 0.007910401967279, 0.00357968641756, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 ) Selected = c(1, 2, 3, 16, 47, 9, 48, 63, 22, 74, 77, 53, 31, 27, 60, 6, 15, 93, 5, 82, 67, 45, 91, 7, 30, 25, 59, 13, 55, 61, 54, 35, 84, 4, 1, 2, 3, 8, 10, 11, 12, 14, 17, 18, 19, 20, 21, 23, 24, 26, 28, 29, 32, 33, 34, 36, 37, 38, 39, 40, 41, 42, 43, 44, 46, 49, 50, 51, 52, 56, 57, 58, 62, 64, 65, 66, 68, 69, 70, 71, 72, 73, 75, 76, 78, 79, 80, 81, 83, 85, 86, 87, 88, 89, 90, 92) ## ----ReanalysisOfHornEtALbarplot,fig.height=5------------------------------ N = 1:dim(PI)[1] bp = barplot(100.0*(1-openVar[N+1]),ylim=c(0,100), ylab="explained variance [in %]",xlab="number query genes", cex.axis=1.5,cex.lab=1.5) axis(side=1,at=bp[N %% 10 == 0],labels=(N)[N %% 10 == 0],cex.axis=1.5) ## ----QueryGeneSelectionLibrary,message=FALSE,results='hide'---------------- library(DmelSGI) basedir = getBaseDir() resultdir = file.path( basedir, "result", "QueryGeneSelection") dir.create(resultdir, recursive = TRUE,showWarnings=FALSE) data("datamatrix", package="DmelSGI") ## ----QueryGeneSelectionLoadLibrary, results='hide',message=FALSE----------- library("DmelSGI") data("SKDdata",package="DmelSGI") data("datamatrix",package="DmelSGI") ## ----QueryGeneSelectionPCA------------------------------------------------- D = apply(SKDdata$D[ ,,1,], c(1,3), mean, na.rm=TRUE) PCA = princomp(D) ## ----QueryGeneSelectionCol------------------------------------------------- col = ifelse(datamatrix$Anno$target$TID %in% datamatrix$Anno$query$TID, "red","gray80") I = order(datamatrix$Anno$target$TID %in% datamatrix$Anno$query$TID) S = PCA$scores S = S[I,] col = col[I] ## ----QueryGeneSelectionPairs----------------------------------------------- par(mar=c(0.2,0.2,0.2,0.2)) pairs(S[,1:5],pch=20,cex=0.7,col=col) ## ----ImageProcessingLibrary,message=FALSE,results='hide'------------------- library("DmelSGI") library("RColorBrewer") basedir = getBaseDir() resultdir = file.path( basedir, "result", "ImageProcessing") dir.create(resultdir, recursive = TRUE,showWarnings=FALSE) data("datamatrix",package="DmelSGI") data("Features",package="DmelSGI") ## ----ImageProcessingGlog,fig.width=4.2,fig.height=4.2---------------------- px = seq(-1.5, 9, length.out=200) trsf = list( log = function(x) log(ifelse(x>0, x, NA_real_)), glog = function(x, c=1) log( (x+sqrt(x^2+c^2))/2 )) colores = c("#202020", "RoyalBlue") matplot(px, sapply(trsf, do.call, list(px)), type="l", lty=c(2,1), col=colores, lwd=2.5, ylab="f(x)", xlab="x") legend("bottomright", fill=colores, legend=names(trsf)) ## ----ImageProcessingBarchartsMainEffectsExampleImagesData------------------ data("mainEffects",package="DmelSGI") ## ----ImageProcessingBarchartsMainEffectsExampleImages---------------------- Main = apply(mainEffects$target,c(1,4),mean,na.rm=TRUE) m = apply(Main,2,mad,center=0.0) for (i in 1:dim(Main)[2]) { Main[,i] = Main[,i] / m[i] } Main = rbind(Main,Fluc=c(0.0,0.0,0.0)) ## ----ImageProcessingExamplePhenotypesSingleKDcol--------------------------- ylim = range(Main[c("Fluc","ida","stg","Arpc1"),1:3]) col = brewer.pal(4,"Pastel1") ## ----ImageProcessingExamplePhenotypesSingleKD,resize.width="0.17\\textwidth",fig.show='hold',fig.width=1.5,fig.height=4---- par(mar=c(0.2,2,2,0.2)) barplot(Main["Fluc",c(1:3,13)],main="Fluc", col=col,ylim=ylim) abline(h=0.0) barplot(Main["ida",c(1:3,13)],main="ida", col=col,ylim=ylim) abline(h=0.0) barplot(Main["stg",c(1:3,13)],main="stg", col=col,ylim=ylim) abline(h=0.0) barplot(Main["Arpc1",c(1:3,13)],main="Arpc1", col=col,ylim=ylim) abline(h=0.0) plot(-10000,xaxt="n",yaxt="n",bty="n",xlim=c(0,1),ylim=c(0,1),xlab="",ylab="") legend("topleft",c("nr cells","MI","area","eccent."),fill=brewer.pal(4,"Pastel1")) ## ----ImageProcessingExamplePhenotypesSingleKDRasGAP1,resize.width="0.2\\textwidth",fig.width=2,fig.height=3---- barplot(Main[c("Fluc","RasGAP1"),1],col=brewer.pal(3,"Pastel1")[1], ylab=c("cell number","z-score"),yaxp=c(0,1,2),las=2) ## ----QualityControlLoadLibrary, results='hide', message=FALSE-------------- library("beeswarm") library("RColorBrewer") library("DmelSGI") library("hwriter") library("xtable") basedir = getBaseDir() resultdir = file.path( basedir, "result", "QualityControl") dir.create(resultdir, recursive = TRUE,showWarnings=FALSE) data("datamatrix", package="DmelSGI") ## ----QualityControlFeatureLoadData----------------------------------------- data("qualityControlFeature", package="DmelSGI") data("Features", package="DmelSGI") ## ----QualityControlFeatureCor---------------------------------------------- Fcor = qualityControlFeature$correlation Fcor = Fcor[order(-Fcor)] ## ----QualityControlFeaturePlot--------------------------------------------- par(mar=c(4.1,4.1,1,1)) plot(Fcor,pch=19,xlab="features", ylab=c("correlation of phenotype","between replicates"), ylim = range(qualityControlFeature$correlation,finite=TRUE), xlim=c(0,sum(is.finite(Fcor))),cex.lab=1.3, cex.axis=1.3) abline(h=0.6) ## ----QualityControlFeatureTableData---------------------------------------- data("Features", package="DmelSGI") data("qualityControlFeature", package="DmelSGI") ## ----QualityControlFeatureTable-------------------------------------------- Features = cbind(Features, QC=ifelse(qualityControlFeature$passed, "passed", "failed"), name = hrNames(row.names(Features))) write.table(Features, file=file.path(resultdir,"FeatureTable.txt"), sep="\t", quote=FALSE) ## ----QualityControlFeatureTableWebpage,echo=FALSE,results='hide'----------- file.copy(system.file("images/hwriter.css",package="hwriter"),resultdir) page=openPage(file.path(resultdir,"FeatureTable.html"), link.css="hwriter.css") hwrite("Features extracted from images", heading=1, page=page) hwrite("[Download as text file]", link="FeatureTable.txt", page=page) hwrite(Features, page=page) closePage(page, splash=FALSE) size = 110 Page = ceiling(seq_len(nrow(Features)) / size) for (p in 1:max(Page)) { XT = xtable(Features[Page==p,], caption=sprintf("Features extracted from images (Part %d/%d)", p,max(Page))) if (p==1) { con = file(file.path(resultdir,"FeatureTable.tex")) label(XT) = "TabFeature" } else { con = file(file.path(resultdir,"FeatureTable.tex"),open="a") writeLines("\\addtocounter{table}{-1}", con=con) } print(XT, file=con, size="tiny", caption.placement="top") close(con) } ## ----QualityControlFeatureTablePrint,echo=FALSE,results='asis'------------- XT = xtable(Features[1:7,], caption="Features extracted from images") print(XT,caption.placement="top") ## ----QualityControlGeneLoadData-------------------------------------------- data("qualityControlGene", package="DmelSGI") ## ----QualityControlGenePlot------------------------------------------------ par(mar=c(4.1,4.5,1,1)) Sample = which(qualityControlGene$Annotation$group == "sample") corGene = qualityControlGene$correlation[Sample] corGene = corGene[order(-corGene)] plot(corGene, pch=19, xlab="targeted genes", ylab=c("cor of phenotypic profile","between dsRNA designs"), cex.lab=1.3, cex.axis=1.3) abline(h=0.7) ## ----QualityControlGeneTableOfPassedGenesTxt------------------------------- data("qualityControlGene", package="DmelSGI") PassedSamples = which((qualityControlGene$Annotation$group == "sample") & (qualityControlGene$passed)) A = qualityControlGene$Annotation A$cor = qualityControlGene$correlation A = A[,c("TID", "Symbol", "cor", "Name")] A = A[PassedSamples,] A = A[order(A$cor),] A$cor = sprintf("%0.2f", A$cor) write.table(A, file=file.path(resultdir,"PassedGenes.txt"), sep="\t", quote=FALSE,row.names=FALSE) ## ----QualityControlGeneTableOfPassedGenesWebpage,echo=FALSE,results='hide'---- file.copy(system.file("images/hwriter.css",package="hwriter"),resultdir) page=openPage(file.path(resultdir,"PassedGenes.html"), link.css="hwriter.css") hwrite("List of genes that passed the quality control", heading=1, page=page) hwrite("[Download as text file]", link="PassedGenes.txt", page=page) hwrite(A, page=page) closePage(page, splash=FALSE) XT = xtable(A,caption="List of genes that passed the quality control") con = file(file.path(resultdir,"PassedGenes.tex")) label(XT) = "TabPassedGenes" print(XT, file=con, size="footnotesize", caption.placement="top") close(con) ## ----QualityControlGeneTableOfPassedGenesPrint,echo=FALSE,results='asis'---- XT = xtable(A[1:7,],caption="List of genes that passed the quality control") print(XT,caption.placement="top") ## ----QualityControlComparisonToRohnEtAl------------------------------------ data("SKDdata", package="DmelSGI") data("mainEffects", package="DmelSGI") data("RohnEtAl", package="DmelSGI") RohnEtAl = RohnEtAl[which(RohnEtAl$Computed.Target %in% SKDdata$Anno$target$TID),] I = match(RohnEtAl$Computed.Target, SKDdata$Anno$target$TID) RohnEtAlanno = RohnEtAl[,1:3] RohnEtAl = as.matrix(RohnEtAl[4:29]) D = apply(SKDdata$D[I,,,],c(1,4), mean, na.rm=TRUE) i=3; j=2 anova(lm(D[,i] ~ RohnEtAl[,j])) beeswarm(D[,i] ~ RohnEtAl[,j],pch=20, xlab=c(colnames(RohnEtAl)[j],"Rohn et al."), ylab=dimnames(D)[[2]][i]) i=13; j=5 anova(lm(D[,i] ~ RohnEtAl[,j])) beeswarm(D[,i] ~ RohnEtAl[,j],pch=20, xlab=c(colnames(RohnEtAl)[j],"Rohn et al."), ylab=dimnames(D)[[2]][i]) ## ----FeatureSelectionLoadLibrary, results='hide', message=FALSE------------ library("DmelSGI") library("RColorBrewer") library("hwriter") library("xtable") basedir = getBaseDir() resultdir = file.path( basedir, "result", "FeatureSelection") dir.create(resultdir, recursive = TRUE,showWarnings=FALSE) ## ----FeatureSelectionLoadData,results='hide'------------------------------- data("stabilitySelection", package="DmelSGI") ## ----FeatureSelectionPlotCor, resize.width="0.7\\textwidth",fig.width=11,fig.height=6---- par(mar=c(12,5,0.5,0.5)) barplot(stabilitySelection$correlation, names.arg=hrNames(stabilitySelection$selected),las=2, col=ifelse(stabilitySelection$ratioPositive > 0.5, brewer.pal(3,"Pastel1")[2], brewer.pal(3,"Pastel1")[1]), ylab = "correlation", cex.lab=1.3) ## ----FeatureSelectionPlotRatio, resize.width="0.7\\textwidth",fig.width=11,fig.height=6---- par(mar=c(12,5,0.5,0.5)) barplot(stabilitySelection$ratioPositive-0.5, names.arg=hrNames(stabilitySelection$selected),las=2, col=ifelse(stabilitySelection$ratioPositive > 0.5, brewer.pal(3,"Pastel1")[2], brewer.pal(3,"Pastel1")[1]), ylab = "ratio positive cor.", cex.lab=1.3, offset=0.5) ## ----FeatureSelectionTableOfFeaturesData----------------------------------- data("stabilitySelection", package="DmelSGI") data("datamatrix", package="DmelSGI") ## ----FeatureSelectionTableOfFeaturesTxt------------------------------------ df = as.data.frame(stabilitySelection[c("selected","correlation","ratioPositive")], stringsAsFactors=FALSE) row.names(df) = 1:nrow(df) df = cbind(df, Name=hrNames(stabilitySelection$selected), Selected = ifelse(stabilitySelection$ratioPositive > 0.5,"Selected",""), stringsAsFactors=FALSE) df = df[,c(1,4,2,3,5)] colnames(df) = c("ID","Name","Correlation","RatioPositive","Selected") write.table(df, file=file.path(resultdir,"StabilitySelectedFeatures.txt"), sep="\t", quote=FALSE,row.names=FALSE) ## ----FeatureSelectionTableOfFeaturesWebpage,echo=FALSE,results='hide'------ file.copy(system.file("images/hwriter.css",package="hwriter"),resultdir) page=openPage(file.path(resultdir,"StabilitySelectedFeatures.html"), link.css="hwriter.css") hwrite("Features selected by stability", heading=1, page=page) hwrite("[Download as text file]", link="StabilitySelectedFeatures.txt", page=page) hwrite(df, page=page) closePage(page, splash=FALSE) XT = xtable(df,caption="List of features selected by stability") label(XT) = "TabStabilitySelection" print(XT, file=file.path(resultdir,"StabilitySelectedFeatures.tex"), caption.placement="top") ## ----FeatureSelectionTable,results='asis'---------------------------------- XT = xtable(df[1:7,],caption="List of features selected by stability") label(XT) = "TabStabilitySelection" print(XT,caption.placement="top") ## ----PIscoresLibrary,message=FALSE,results='hide'-------------------------- library("DmelSGI") library("RColorBrewer") basedir = getBaseDir() resultdir = file.path( basedir, "result", "PairwiseInteractionScores") dir.create(resultdir, recursive = TRUE,showWarnings=FALSE) data("datamatrix", package="DmelSGI") ## ----PIscoresExamplePhenotypesDoubleKDdata--------------------------------- data("mainEffects",package="DmelSGI") data("pimatrix",package="DmelSGI") ## ----PIscoresExamplePhenotypesDoubleKD------------------------------------- examples = data.frame( ph = c("area","mitoticIndex"), targetGenes = c("FBgn0014020","FBgn0033029"), queryGenes=c("FBgn0030276","FBgn0261456"), stringsAsFactors=FALSE) Effects = array(0, dim=c(nrow(examples), 4, nrow(pimatrix$Anno$phenotype))) for (i in seq_len(nrow(examples))) { I = match(examples$targetGenes[i], pimatrix$Anno$target$TID) J = match(examples$queryGenes[i], pimatrix$Anno$query$TID) B = pimatrix$Anno$query$Batch[J] TP = pimatrix$Anno$target$TargetPlate[I] Effects[i,1,] = apply(mainEffects$target[I,,B,],2, mean, na.rm=TRUE) Effects[i,2,] = apply(mainEffects$query[J,,TP,], 2, mean, na.rm=TRUE) Effects[i,3,] = Effects[i,1,] + Effects[i,2,] Effects[i,4,] = apply(pimatrix$D[I,,J,,],3,mean) + Effects[i,3,] } for (i in seq_len(nrow(examples))) { tg = pimatrix$Anno$target$Symbol[match(examples$targetGenes[i],pimatrix$Anno$target$TID)] qg = pimatrix$Anno$query$Symbol[match(examples$queryGenes[i], pimatrix$Anno$query$TID)] pdf(file.path(resultdir,sprintf("ExamplePhenotypes-doubleKD-%s-%s.pdf", tg,qg)),height=8) par(mfrow=c(1,2),mar=c(18,3,3,0.2)) K=1 bp = barplot(Effects[i,,K],main="number cells", col=brewer.pal(3,"Pastel1")[1],cex.main=2,cex.axis=1.5) abline(h=0.0) lines(c(bp[4]-0.5,bp[4]+0.5),c(Effects[i,3,K],Effects[i,3,K])) arrows(x0=bp[4],Effects[i,3,K],x1=bp[4],Effects[i,4,K],code=3, length=min(0.25,abs(Effects[i,4,K] - Effects[i,3,K]))) axis(side=1,at=bp,labels=c(sprintf("%s+ctrl.",tg), sprintf("%s+ctrl.",qg),"expected", sprintf("%s+%s",tg,qg)),col=NA,cex.axis=1.8,las=2) if (examples$ph[i] == "area") { K = 5 } else { K = 2 } bp = barplot(Effects[i,,K],main=ifelse(K==2,"mitotic index","nuclear area"), col=brewer.pal(3,"Pastel1")[2],cex.main=2,cex.axis=1.5) abline(h=0.0) lines(c(bp[4]-0.5,bp[4]+0.5),c(Effects[i,3,K],Effects[i,3,K])) arrows(x0=bp[4],Effects[i,3,K],x1=bp[4],Effects[i,4,K],code=3) axis(side=1,at=bp,labels=c(sprintf("%s + ctrl.",tg), sprintf("%s + ctrl.",qg),"expected", sprintf("%s + %s",tg,qg)),col=NA, cex.axis=1.8,las=2) dev.off() } ## ----PIscoresExamplePhenotypesDoubleKDSWISNFdata--------------------------- data("Interactions", package="DmelSGI") f = "4x.count" names = c("mor","brm","Bap60","Snr1","osa") ## ----PIscoresExamplePhenotypesDoubleKDSWISNF,resize.width="0.4\\textwidth",fig.width=6,fig.height=4.5---- bp=barplot(Interactions$piscore[names, "RasGAP1", f], ylim=c(-0.6,0.6),ylab=sprintf("pi-score (%s)",hrNames(f)), las=2,cex.axis=1.5,cex.names=1.5,cex.lab=1.5,yaxp=c(-0.5,0.5,2)) abline(h=0) ## ----PIscoresMainResultTable1---------------------------------------------- library("DmelSGI") data("Interactions",package="DmelSGI") ## ----PIscoresMainResultTable2---------------------------------------------- PI = Interactions$piscore PADJ = Interactions$padj PI[is.na(PADJ)] = NA ## ----PIscoresMainResultTable3---------------------------------------------- dim(PI) = c(prod(dim(PI)[1:2]),dim(PI)[3]) dim(PADJ) = c(prod(dim(PADJ)[1:2]),dim(PADJ)[3]) ## ----PIscoresMainResultTable4---------------------------------------------- V = cbind(PI, PADJ) V = V[,rep(seq_len(dim(PI)[2]),each=2)+rep(c(0,dim(PI)[2]),times=dim(PI)[2])] colnames(V) = sprintf("%s.%s",rep(c("pi-score","padj"),times=dim(PI)[2]), rep(hrNames(Interactions$Anno$phenotype$phenotype), each=2)) ## ----PIscoresMainResultTable5---------------------------------------------- target = rep(Interactions$Anno$target$Symbol, times=dim(Interactions$piscore)[2]) query = rep(Interactions$Anno$query$Symbol, each=dim(Interactions$piscore)[1]) df = data.frame(targetGene=target, queryGene=query, V) write.table(df, file=file.path(resultdir,"interactions.txt"),sep="\t", row.names=FALSE,quote=FALSE) ## ----PIscoresGIQualityControlMainEffectsAcrossPhenotypes------------------- data("mainEffects", package="DmelSGI") D = apply(mainEffects$target, c(1,4), mean, na.rm=TRUE) colnames(D) = hrNames(colnames(D)) ## ----PIscoresGIQualityControl-mainEffects-Features1-11,resize.width="0.7\\textwidth",fig.width=11.5,fig.height=11.5---- par(mar=c(0.2,0.2,0.2,0.2)) pairs(D[,1:11],pch=20,cex=0.5) ## ----PIscoresGIQualityControl-mainEffects-Features12-21,resize.width="0.7\\textwidth",fig.width=11.5,fig.height=11.5---- par(mar=c(0.2,0.2,0.2,0.2)) pairs(D[,c(1,12:21)],pch=20,cex=0.5) ## ----PIscoresGIQualityControlMainEffectsAcrossBatches---------------------- data("mainEffects", package="DmelSGI") data("pimatrix", package="DmelSGI") D = apply(mainEffects$target[,,,2], c(1,3), mean, na.rm=TRUE) ## ----PIscoresGIQualityControl-mainEffect-batches1to12,resize.width="0.7\\textwidth",fig.width=11.5,fig.height=11.5---- par(mar=c(0.2,0.2,0.2,0.2)) pairs(D,pch=20,cex=0.5) ## ----PIscoresGIQualityControl-mainEffect-batches1and12--------------------- par(mar=c(4.5,4.5,1,1)) plot(D[,1],D[,12],pch=20, xlab="mitotic index [glog] in batch 1", ylab="mitotic index [glog] in batch 12", main="",cex=1.5,cex.lab=2,cex.axis=2,cex.main=2) text(x=1,y=-2.0,sprintf("cor = %0.2f",cor(D[,3],D[,4])),cex=2) ## ----PIscoresGIQualityControlPIscoreAcrossPhenotypes1---------------------- data("Interactions", package="DmelSGI") D = Interactions$piscore dim(D) = c(prod(dim(D)[1:2]),dim(D)[3]) colnames(D) = hrNames(Interactions$Anno$phenotype$phenotype) ## ----PIscoresGIQualityControlPIscoreAcrossPhenotypes2---------------------- set.seed(1043289201) S = sample(1:dim(D)[1],1000) D1 = D[S,] ## ----PIscoresGIQualityControl-piscore-Features1-11,resize.width="0.7\\textwidth",fig.width=11.5,fig.height=11.5---- par(mar=c(0.2,0.2,0.2,0.2)) pairs(D1[,1:11],pch=20,cex=0.5) ## ----PIscoresGIQualityControl-piscore-Features12-21,resize.width="0.7\\textwidth",fig.width=11.5,fig.height=11.5---- par(mar=c(0.2,0.2,0.2,0.2)) pairs(D1[,c(1,12:21)],pch=20,cex=0.5) ## ----PIscoresGIQualityControlPIscoreAcrossPhenotypes3---------------------- set.seed(1043289201) S = sample(1:dim(D)[1],2000) D1 = D[S,] ## ----PIscoresGIQualityControlPIscoreAcrossPhenotypes4---------------------- colnames(D1) = hrNames(colnames(D1)) for (i in c(2,16)) { X = D1[,c(1,i)] pdf(file=file.path(resultdir,sprintf("GeneticInteractionQC-piscore-cellnumber-%s.pdf", gsub("[ ,()]","",colnames(X)[2])))) s = mad(X[,2],na.rm=TRUE) / mad(X[,1],na.rm=TRUE) r = max(abs(X[,1]),na.rm=TRUE) r = c(-r,r) X[which(X[,2] > s*r[2]),2] = s*r[2] X[which(X[,2] < s*r[1]),2] = s*r[1] par(mar=c(4.5,4.5,1,1)) plot(X[,2],X[,1],pch=20, xlab=sprintf("pi-score (%s)",colnames(X)[2]), ylab=sprintf("pi-score (%s)",colnames(X)[1]), main="",cex=1.5,cex.lab=2,cex.axis=2,cex.main=2,xlim=s*r,ylim=r) dev.off() cat("correlation nrCells - ",dimnames(D1)[[2]][i]," = ", cor(X[,2],X[,1],use="pairwise.complete"),"\n") } ## ----PIscoresGIQualityControlPIscoreAcrossPhenotypes5---------------------- data("pimatrix", package="DmelSGI") D = pimatrix$D D1 = (D[,1,,1,] + D[,1,,2,]) / 2.0 D2 = (D[,2,,1,] + D[,2,,2,]) / 2.0 dim(D1) = c(prod(dim(D1)[1:2]),dim(D1)[3]) dim(D2) = c(prod(dim(D2)[1:2]),dim(D2)[3]) colnames(D1) = colnames(D2) = hrNames(pimatrix$Anno$phenotype$phenotype) for (i in c(1,2,16)) { cc = cor(D1[,i],D2[,i],use="pairwise.complete") cat("correlation between replicates ",dimnames(D1)[[2]][i]," = ",cc,"\n") } ## ----PIscoresGIQualityControlPIscoreNumberInteractions1-------------------- data("Interactions", package="DmelSGI") p.value.cutoff = 0.01 N = matrix(NA_integer_, nr=nrow(Interactions$Anno$phenotype), nc=2) colnames(N) = c("pos","neg") for (i in 1:nrow(Interactions$Anno$phenotype)) { PI = Interactions$piscore[,,i] N[i,2] = sum(PI[Interactions$padj[,,i] <= p.value.cutoff] > 0) N[i,1] = sum(PI[Interactions$padj[,,i] <= p.value.cutoff] < 0) } ## ----PIscoresGIQualityControlPIscoreNumberInteractions2-------------------- N = N / prod(dim(Interactions$piscore)[1:2]) ## ----PIscoresGIQualityControl-nrOfInteractions,fig.width=7,fig.height=6---- par(mar=c(15,5,0.5,0.5)) barplot(t(N),col=c("cornflowerblue","yellow"), names.arg=hrNames(Interactions$Anno$phenotype$phenotype), las=2, ylab ="fraction of interactions", cex.names=1.2,cex.axis=1.5,cex.lab=1.5) ## ----PIscoresGIQualityControlPIscoreNumberInteractions3-------------------- isinteraction = rep(FALSE, prod(dim(Interactions$piscore)[1:2])) Ncum = rep(NA_integer_, nrow(Interactions$Anno$phenotype)) for (i in 1:nrow(Interactions$Anno$phenotype)) { isinteraction[Interactions$padj[,,i] <= p.value.cutoff] = TRUE Ncum[i] = sum(isinteraction) / prod(dim(Interactions$piscore)[1:2]) } ## ----PIscoresGIQualityControl-nrOfInteractionsCumulative,fig.width=7,fig.height=6---- par(mar=c(15,5,0.5,0.5),xpd=NA) bp=barplot(Ncum, col=brewer.pal(3,"Pastel1")[2], ylab="fraction of interactions",las=2, names.arg = rep("",length(Ncum)),cex.axis=1.5,cex.lab=1.5) text(bp,-0.01,hrNames(Interactions$Anno$phenotype$phenotype),adj=c(1,0.5),srt=38,cex=1.2) ## ----GeneticInteractionCubeLoadLibrary,results='hide',message=FALSE-------- library("DmelSGI") library("grid") library("RColorBrewer") library("gplots") library("beeswarm") basedir = getBaseDir() resultdir = file.path( basedir, "result", "GeneticInteractionCube") dir.create(resultdir, recursive = TRUE,showWarnings=FALSE) dir.create(file.path("Figures","GeneticInteractionCube"), showWarnings=FALSE) ## ----GeneticInteractionCubeData-------------------------------------------- data("Interactions", package="DmelSGI") ## ----GeneticInteractionCubeNormalize1,eval=FALSE--------------------------- # data("pimatrix", package="DmelSGI") # D = pimatrix$D # D2 = aperm(D, c(1,3,2,4,5)) # dim(D2) = c(prod(dim(D2)[1:2]),prod(dim(D2)[3:4]),dim(D2)[5]) # # SD = apply(D2,c(1,3),sd, na.rm=TRUE) # MSD = apply(SD, 2, function(x) { median(x,na.rm=TRUE) } ) ## ----GeneticInteractionCubeNormalize2-------------------------------------- MSD = c(0.0833528975114521, 0.134136618342975, 0.0498996012784751, 0.204772216536139, 0.0142975582945938, 0.0428299793772605, 0.0576235314621808, 0.0833934805305705, 0.0328437541652814, 0.147643254412127, 0.0866394118952878, 0.140840565863283, 0.0154131573539473, 0.0286467941877466, 0.0496616658001497, 0.0164694485385577, 0.233130597062897, 0.222961290060361, 0.00228512594775289, 0.0773453995034531, 0.0892678802977647) D = Interactions$piscore for (i in 1:dim(D)[3]) { D[,,i] = D[,,i] / MSD[i] } ## ----GeneticInteractionCubeCut--------------------------------------------- cuts = c(-Inf, seq(-6, -2, length.out=(length(DmelSGI:::colBY)-3)/2), 0.0, seq(2, 6, length.out=(length(DmelSGI:::colBY)-3)/2), +Inf) ## ----GeneticInteractionCubeHClust------------------------------------------ ordTarget = orderDim(D, 1) ordQuery = orderDim(D, 2) ordFeat = orderDim(D, 3) Ph = c("4x.intNucH7","4x.count","4x.LCD2", "4x.ratioMitotic", "10x.meanNonmitotic.nucleus.DAPI.m.majoraxis", "4x.areaNucAll", "10x.meanNonmitotic.cell.Tub.m.eccentricity", "4x.areapH3All", "4x.intH3pH4") D1 = D[ordTarget, ordQuery, Ph] dimnames(D1)[[3]] = hrNames(dimnames(D1)[[3]]) ## ----GeneticInteractionCubeDPiMlibrary------------------------------------- library("DmelSGI") library("RColorBrewer") library("gplots") library("grid") data("Interactions",package="DmelSGI") data("FBgn2anno",package="DmelSGI") data("DPiM",package="DmelSGI") ## ----GeneticInteractionCubeDPiMcor----------------------------------------- PI = Interactions$piscore dim(PI) = c(dim(PI)[1],prod(dim(PI)[2:3])) C = cor(t(PI)) row.names(C) = colnames(C) = Interactions$Anno$target$TID ## ----GeneticInteractionCubeDPiMsets---------------------------------------- m1 = match(DPiM$interactions$Interactor_1, Interactions$Anno$target$TID) m2 = match(DPiM$interactions$Interactor_2, Interactions$Anno$target$TID) I = which(!is.na(m1) & !is.na(m2) & (m1 != m2)) I = cbind(m1[I],m2[I]) ccDPiM = C[I] Iall = upper.tri(C) Iall[I] = FALSE Iall[cbind(I[,2],I[,1])] = FALSE ccall = C[Iall] ## ----GeneticInteractionCubeDPiM,fig.width=10,fig.height=4------------------ par(mar=c(5,5,0.2,5),xpd=NA) d1 = density(ccDPiM, from=-1,to=1) d2 = density(ccall, from=-1,to=1) plot(d2, main="", col=brewer.pal(9, "Set1")[2],lwd=5, xlab="correlation of interaction profiles",ylim=c(0,1.5), cex.lab=1.8,cex.axis=1.8) lines(d1$x,d1$y,col=brewer.pal(9, "Set1")[1],lwd=5) legend("topleft",c("co-purified in DPiM","not captured by DPiM"),bty="n", fill = brewer.pal(9, "Set1")[1:2],cex=1.8) ## ----GIlandscapeLoadLibrary,results='hide',message=FALSE------------------- library("DmelSGI") library("igraph") library("RSVGTipsDevice") library("hwriter") library("grid") basedir = getBaseDir() resultdir = file.path( basedir, "result", "GeneticInteractionLandscape") dir.create(resultdir, recursive = TRUE,showWarnings=FALSE) ## ----GIlandscapeLoadData--------------------------------------------------- data("Interactions", package="DmelSGI") ## ----GIlandscapeNormalization1,eval=FALSE---------------------------------- # data("pimatrix", package="DmelSGI") # D = pimatrix$D # D2 = aperm(D, c(1,3,2,4,5)) # dim(D2) = c(prod(dim(D2)[1:2]),prod(dim(D2)[3:4]),dim(D2)[5]) # # SD = apply(D2,c(1,3),sd, na.rm=TRUE) # MSD = apply(SD, 2, function(x) { median(x,na.rm=TRUE) } ) ## ----GIlandscapeNormalization2--------------------------------------------- Sel = 1:1293 MSD = c(0.0833528975114521, 0.134136618342975, 0.0498996012784751, 0.204772216536139, 0.0142975582945938, 0.0428299793772605, 0.0576235314621808, 0.0833934805305705, 0.0328437541652814, 0.147643254412127, 0.0866394118952878, 0.140840565863283, 0.0154131573539473, 0.0286467941877466, 0.0496616658001497, 0.0164694485385577, 0.233130597062897, 0.222961290060361, 0.00228512594775289, 0.0773453995034531, 0.0892678802977647) ## ----GIlandscapeNormalization3--------------------------------------------- PI = Interactions$piscore for (i in 1:dim(PI)[3]) { PI[,,i] = PI[,,i] / MSD[i] } dim(PI) = c(dim(PI)[1],prod(dim(PI)[2:3])) ## ----GIlandscapeNrInteractions--------------------------------------------- nint = apply(Interactions$padj <= 0.01,1,sum) Nint = matrix(nint,nrow=length(nint),ncol=length(nint)) Nint[t(Nint) < Nint] = t(Nint)[t(Nint) < Nint] row.names(Nint) = colnames(Nint) = Interactions$Anno$target$Symbol ## ----GIlandscapeLoadGeneClusters------------------------------------------- data("SelectedClusters", package="DmelSGI") SelectedClusters = c(SelectedClusters, list(allOthers = Interactions$Anno$target$Symbol[ !( Interactions$Anno$target$Symbol %in% unlist(SelectedClusters))])) Labels = SelectedClusters ## ----GIlandscapeCol-------------------------------------------------------- set.seed(26323) Col = c("SWI/SNF" = "#ED1C24", "Condensin / Cohesin" = "#B8D433", "Cytokinesis" = "#67BF6B", "SAC" = "#3C58A8", "DREAM complex" = "#B64F9D", "Centrosome / Mitotic spindle" = "#EF4123", "CCT" = "#8CC63F", "Sequence-specific TFs" = "#67C4A4", "26S Proteasome" = "#3A51A3", "CSN" = "#F17A22", "RNA helicase" = "#6CBE45", "APC/C" = "#66C9D7", "Ribosomal biogenesis" = "#4A50A2", "Condensin / Cohesin (2)" = "#ED127A", "SAGA & Mediator" = "#F1B61C", "Cell-cell signalling" = "#61BB46", "Vesicle trafficking and cytoskeleton" = "#2CB4E8", "DNA repair and apoptosis" = "#6950A1", "ARP2/3 complex" = "#ED1940", "Tor signalling" = "#D54097", "Ras / MAPK signalling" = "#65BC46", "RNA PolII" = "#4074BA", "Wnt signalling" = "#8E4F9F") col = rep("gray80",nrow(Interactions$Anno$target)) names(col) = Interactions$Anno$target$Symbol for (i in 1:length(Labels)) { col[Labels[[i]]] = Col[i] } ## ----GIlandscapePCA-------------------------------------------------------- dimPCA = 25 PCA = prcomp(PI) X = sweep(PI,2,PCA$center) %*% PCA$rotation[,1:dimPCA] X = sweep(X,2,apply(X,2,sd), FUN="/") theCorrPCA = cor(t(X), use = "pairwise.complete.obs") row.names(theCorrPCA) = colnames(theCorrPCA) = Interactions$Anno$target$Symbol ## ----GIlandscapeDist------------------------------------------------------- D = 2 - 2*theCorrPCA D[lower.tri(D,diag=TRUE)] = NA ## ----GIlandscapeEdgeList--------------------------------------------------- thresholdDist = 0.8 wedges = data.frame(V1 = rep(Interactions$Anno$target$Symbol,times=dim(D)[1]), V2 = rep(Interactions$Anno$target$Symbol,each=dim(D)[1]), nint = as.vector(Nint), weight = as.vector(D),stringsAsFactors=FALSE) wedges = wedges[which(wedges$weight <= thresholdDist),] ## ----GIlandscapeigraph----------------------------------------------------- g <- graph.data.frame(wedges, directed=FALSE) V(g)$color = col[V(g)$name] V(g)$frame.color = ifelse(V(g)$name %in% SelectedClusters$allOthers, "#666666","#000000") V(g)$size = 1.5 V(g)$size[!(V(g)$name %in% Labels$allOthers)] = 2.5 V(g)$label = rep("",length(V(g)$name)) E(g)$color <- "#e7e7e7" E(g)$color[E(g)$nint > 5] <- "#cccccc" ## ----GIlandscapeFruchtermanReingold---------------------------------------- set.seed(234816) a = 0.07 b = 2.0 co <- layout.fruchterman.reingold(graph=g, params=list(weights=(thresholdDist - E(g)$weight), area=a*vcount(g)^2, repulserad=b*a*vcount(g)^2*vcount(g))) ## ----GIlandscapeUnitBox---------------------------------------------------- co[,1] = 2*(co[,1] - min(co[,1])) / diff(range(co[,1]))-1 co[,2] = 2*(co[,2] - min(co[,2])) / diff(range(co[,2]))-1 row.names(co) = V(g)$name ## ----GIlandscapePermuteVertices-------------------------------------------- g = permute.vertices(graph=g, permutation = rank((!(V(g)$name %in% SelectedClusters$allOthers)), ties.method="random")) co = co[V(g)$name,] ## ----GIlandscape-graph,resize.width="0.9\\textwidth",fig.width=14,fig.height=10---- # par(mar=c(0.1,0.1,0.1,0.1)) plot(g, layout=co) plotHairballLabels(g, co, Labels[-length(Labels)], Col) ## ----GIlandscapeHairballSVG,echo=FALSE,results='hide',eval=FALSE----------- # # Plot the genetic interaction landscape as SVG # devSVGTips(file.path(resultdir,"GIlandscape-graph.svg"), toolTipMode=1, # title="Genetic interaction landscape", # width = 10, height = 10) # par(xpd=NA,mar=c(5,5,5,5)) # plot(0,type='n', # xlim=range(co, na.rm=TRUE), # ylim=range(co, na.rm=TRUE), # xlab='', ylab='', # main="Genetic interaction landscape", # xaxt="n",yaxt="n") # invisible(sapply(1:(dim(co)[1]), function(i) { # setSVGShapeToolTip(desc=V(g)$name[i]) # points(co[i,1], co[i,2],bg=col[V(g)$name[i]],cex=1.5,pch=21) # })) # dev.off() ## ----GIlandscapeHairballHwriter,echo=FALSE,results='hide'------------------ # Write a HTML webpage to access the SVG graphic Lpng = "GIlandscape-graph.png" Lpdf = "
[pdf] - [svg]
" M = sprintf("%s
%s",hwriteImage(Lpng,table=FALSE),Lpdf) file.copy(system.file(file.path("images","hwriter.css"),package="hwriter"),file.path(resultdir,"hwriter.css")) page = openPage(file.path(resultdir,"GIlandscape-index.html"),link.css="hwriter.css") hwrite(M, page=page) hwriteImage("GIlandscape-phenotypes.png",link="GIlandscape-phenotypes.pdf",page=page) closePage(page,splash=FALSE) ## ----GIlandscapeHairballZoominA1------------------------------------------- A = SelectedClusters[c("APC/C","SAC","Centrosome / Mitotic spindle", "Condensin / Cohesin","26S Proteasome")] genesA = c("Arp10","fzy","vih","Klp61F","polo") gsubA = induced.subgraph(g, which(V(g)$name %in% c(unlist(A), genesA))) ## ----GIlandscapeHairballZoominA2------------------------------------------- data("SelectedClustersComplexes",package="DmelSGI") V(gsubA)$color[which(V(gsubA)$name %in% SelectedClustersComplexes$gammaTuRC)] = "orange2" V(gsubA)$color[which(V(gsubA)$name %in% SelectedClustersComplexes$'Dynein/Dynactin')] = "yellow2" V(gsubA)$label = V(gsubA)$name V(gsubA)$size = 8 V(gsubA)$label.cex = 0.4 V(gsubA)$label.color = "#222222" E(gsubA)$color = "#777777" ## ----GIlandscape-graph-zoomin-A-------------------------------------------- set.seed(38383) LA = layout.fruchterman.reingold(gsubA,params=list(area=121^2*40)) plot(gsubA,layout=LA) ## ----GIlandscapeHairballZoominB1------------------------------------------- B = SelectedClusters[c("Tor signalling","Ribosomal biogenesis","RNA helicase", "Ribosomal biogenesis","RNA helicase", "DNA repair and apoptosis","Cell-cell signalling", "Vesicle trafficking and cytoskeleton")] B$'Tor signalling' = B$'Tor signalling'[!(B$'Tor signalling' %in% c("trc","InR","gig","Tsc1","Pten"))] genesB = c("Dbp45A","hpo","14-3-3epsilon", "CG32344", "CG9630", "kz", "pit", "Rs1","CG8545","twin") gsubB = induced.subgraph(g, which(V(g)$name %in% c(unlist(B), genesB))) ## ----GIlandscapeHairballZoominB2------------------------------------------- PolI = c("CG3756","l(2)37Cg","RpI1","RpI12","RpI135","Tif-IA") PolIII = c("Sin","RpIII128","CG5380","CG12267","CG33051","CG7339") V(gsubB)$color[which(V(gsubB)$name %in% PolI)] = "orange2" V(gsubB)$color[which(V(gsubB)$name %in% PolIII)] = "yellow2" V(gsubB)$label = V(gsubB)$name V(gsubB)$size = 8 V(gsubB)$label.cex = 0.4 V(gsubB)$label.color = "#222222" E(gsubB)$color = "#777777" ## ----GIlandscape-graph-zoomin-B-------------------------------------------- set.seed(122138) LB = layout.fruchterman.reingold(gsubB,params=list(area=121^2*40)) plot(gsubB,layout=LB) ## ----DirectionalInteractionsLibrary,results='hide',message=FALSE----------- library("DmelSGI") library("abind") library("igraph") basedir = getBaseDir() resultdir = file.path( basedir, "result", "DirectionalInteractions") dir.create(resultdir, recursive = TRUE,showWarnings=FALSE) data("Interactions", package="DmelSGI") data("pimatrix", package="DmelSGI") data("mainEffects", package="DmelSGI") data("SelectedClustersComplexes", package="DmelSGI") ## ----DirectionalInteractionsMainEffects------------------------------------ pi = pimatrix$D mt = mainEffects$target mq = mainEffects$query ## ----DirectionalInteractionsNormalize-------------------------------------- myrange = function(x) diff(quantile(x, probs=c(0.01, 0.99), na.rm=TRUE)) featureSD = apply(pi, 5, myrange) for(k in seq(along=featureSD)){ pi[,,,,k] = pi[,,,,k]/featureSD[k] mt[,,,k] = mt[,,,k]/featureSD[k] mq[,,,k] = mq[,,,k]/featureSD[k] } ## ----DirectionalInteractionsRearrage--------------------------------------- data = array(NA_real_, dim=c(21, 3, dim(pi)[1:4])) for(it in seq_len(dim(pi)[1])) { targetplate = pimatrix$Anno$target$TargetPlate[ it ] for(dt in seq_len(dim(pi)[2])) { for(iq in seq_len(dim(pi)[3])) { batch = pimatrix$Anno$query$Batch[iq] xt = mt[it, dt, batch, ] for(dq in seq_len(dim(pi)[4])) { y = pi[it, dt, iq, dq, ] xq = mq[iq, dq, targetplate, ] nay = sum(is.na(y)) naxq = sum(is.na(xq)) if((nay>1)||(naxq>1)) { } else { data[, 1, it, dt, iq, dq] = xt data[, 2, it, dt, iq, dq] = xq data[, 3, it, dt, iq, dq] = y } # else } } } } dimnames(data) = list(pimatrix$Anno$phenotype$phenotype, c("xt","xq","pi"), pimatrix$Anno$target$Symbol, 1:2, pimatrix$Anno$query$Symbol, 1:2) ## ----DirectionalInteractionsFit,eval=FALSE--------------------------------- # resCoef = array(NA_real_, dim=c(3, dim(pi)[1:4])) # resSq = array(NA_real_, dim=c(3, dim(pi)[1:4])) # resPV = array(NA_real_, dim=c(2, dim(pi)[1:4])) # for(it in seq_len(dim(pi)[1])) { # for(dt in seq_len(dim(pi)[2])) { # for(iq in seq_len(dim(pi)[3])) { # for(dq in seq_len(dim(pi)[4])) { # if (all(is.finite(data[, , it, dt, iq, dq]))) { # model = lm(data[,3,it,dt,iq,dq] ~ # data[,1,it,dt,iq,dq]+data[,2,it,dt,iq,dq]) # a = anova(model) # resCoef[, it, dt, iq, dq] = model$coefficients # resSq[1, it, dt, iq, dq] = a[1,2] # resSq[2, it, dt, iq, dq] = a[2,2] # resSq[3, it, dt, iq, dq] = a[3,2] # resPV[1, it, dt, iq, dq] = a[1,5] # resPV[2, it, dt, iq, dq] = a[2,5] # } # else # } # } # } # } # # dimnames(resCoef) = list(c("const","xt","xq"), # pimatrix$Anno$target$Symbol, # 1:2, # pimatrix$Anno$query$Symbol, # 1:2) # # dimnames(resSq) = list(c("xt","xq","res"), # pimatrix$Anno$target$Symbol, # 1:2, # pimatrix$Anno$query$Symbol, # 1:2) # # dimnames(resPV) = list(c("xt","xq"), # pimatrix$Anno$target$Symbol, # 1:2, # pimatrix$Anno$query$Symbol, # 1:2) # # fitepistasis = list(Coef = resCoef, Sq = resSq) # # save(fitepistasis, file="fitepistasis.rda") ## ----DirectionalInteractionsLoadFit,results='hide'------------------------- data("fitepistasis", package="DmelSGI") ## ----DirectionalInteractions-example2pheno-sti-Cdc23,fig.width=5*1.3,fig.height=4*1.3---- plot2Phenotypes(data,"sti","Cdc23",1,5,lwd=3,cex.axis=1.3,cex.lab=1.3,cex=1.3) ## ----DirectionalInteractions-example21pheno-sti-Cdc23,fig.width=8,fig.height=4---- plotPIdata(data,"sti","Cdc23",cex.axis=1.3,cex.lab=1.3) ## ----DirectionalInteractionsResSQ------------------------------------------ resSq = fitepistasis$Sq x = apply(resSq,2:5,sum, na.rm=TRUE) resSq[1,,,,] = resSq[1,,,,] / x resSq[2,,,,] = resSq[2,,,,] / x resSq[3,,,,] = resSq[3,,,,] / x ## ----DirectionalInteractionsTH--------------------------------------------- NSIG = (apply(Interactions$padj <= 0.01,1:2,sum,na.rm=TRUE)) SQ = apply(resSq[,,,,], c(1,2,4),mean,na.rm=TRUE) Coef = apply(fitepistasis$Coef[,,,,], c(1,2,4),mean,na.rm=TRUE) thT = quantile(as.vector(SQ["xt",,]),probs=c(0.1,0.95)) thQ = quantile(as.vector(SQ["xq",,]),probs=c(0.1,0.95)) ## ----DirectionalInteractionsEdges------------------------------------------ IT = which((NSIG >= 1) & (SQ["xt",,] > thT[2]) & (SQ["xq",,] < thQ[1]),arr.ind=TRUE) IQ = which((NSIG >= 1) & (SQ["xq",,] > thQ[2]) & (SQ["xt",,] < thT[1]),arr.ind=TRUE) IT = cbind(Interactions$Anno$target$Symbol[IT[,1]], Interactions$Anno$query$Symbol[IT[,2]]) IQ = cbind(Interactions$Anno$target$Symbol[IQ[,1]], Interactions$Anno$query$Symbol[IQ[,2]]) ET = data.frame(geneFrom = IT[,2], geneTo = IT[,1], sign = sign((Coef["xt",,])[IT]), coef = (Coef["xt",,])[IT], coefRev = NA, sqFrom = (SQ["xq",,])[IT], sqTo = (SQ["xt",,])[IT], mode = "target", stringsAsFactors=FALSE) ITREV = cbind(IT[,2],IT[,1]) B = (ITREV[,1] %in% dimnames(Coef)[[2]]) & (ITREV[,2] %in% dimnames(Coef)[[3]]) ET$coefRev[B] = (Coef["xq",,])[ITREV[B,]] EQ = data.frame(geneFrom = IQ[,1], geneTo = IQ[,2], sign = sign((Coef["xq",,])[IQ]), coef = (Coef["xq",,])[IQ], coefRev = NA, sqFrom = (SQ["xt",,])[IQ], sqTo = (SQ["xq",,])[IQ], mode = "query", stringsAsFactors=FALSE) IQREV = cbind(IQ[,2],IQ[,1]) B = (IQREV[,1] %in% dimnames(Coef)[[2]]) & (IQREV[,2] %in% dimnames(Coef)[[3]]) EQ$coefRev[B] = (Coef["xt",,])[IQREV[B,]] edges = rbind(ET,EQ) edges$color = ifelse(edges$sign<0, "dodgerblue", "crimson") ## ----DirectionalInteractionsKey-------------------------------------------- key = sprintf("%s__%s",edges$geneFrom,edges$geneTo) k = key[duplicated(key)] g = tapply(edges$sign[key %in% k],key[key %in% k], function(x) { length(unique(x)) > 1}) if (any(g)) { cat(sum(key %in% names(which(g))), " edges with contraditing sign are removed.\n") edges = edges[!(key %in% names(which(g))),] } key = sprintf("%s__%s",edges$geneFrom,edges$geneTo) edges = edges[!duplicated(key),] key = sprintf("%s__%s",edges$geneFrom,edges$geneTo) key2 = sprintf("%s__%s",edges$geneTo,edges$geneFrom) if (any(key %in% key2)) { cat(sum(key %in% key2)," edges with contraditing direction are removed.\n") edges = edges[!((key %in% key2) | (key2 %in% key)),] } ## ----DirectionalInteractionsTable------------------------------------------ E = data.frame(edges[,c("geneFrom","geneTo")], effect=ifelse( edges$sign == 1, "aggravating", "alleviating")) write.table(E, file=file.path(resultdir,"DirectionalEpistaticInteractions.txt"), sep="\t",row.names=FALSE,quote=FALSE) ## ----DirectionalInteractionsSelectedClusters------------------------------- SelectedClustersComplexes = SelectedClustersComplexes[c("Cytokinesis","Condensin/Cohesin", "SAC","Apc/C","gammaTuRC","Dynein/Dynactin")] SelectedClustersComplexes$'Dynein/Dynactin' = c(SelectedClustersComplexes$'Dynein/Dynactin', "Klp61F") SelectedClustersComplexes$polo = "polo" SelectedClustersComplexes[["vih"]] = "vih" SelectedClustersComplexes[["Elongin-B"]] = "Elongin-B" SelectedClustersComplexes[["Skp2"]] = "Skp2" QG = unlist(SelectedClustersComplexes) edges2 = edges edges2 = edges2[(edges2$geneFrom %in% QG) & (edges2$geneTo %in% QG),] nodes2 = list() for (cl in names(SelectedClustersComplexes)) { nodes2[[cl]] = unique(c(edges2$geneFrom, edges2$geneTo)) nodes2[[cl]] = nodes2[[cl]][nodes2[[cl]] %in% SelectedClustersComplexes[[cl]] ] } cat(sprintf("Writing graph with %d nodes and %d edges.\n", length(nodes2), nrow(edges2))) out = c("digraph DirectionalInteractions {", paste("graph [size=\"10,10\" ratio=0.35 mode=major outputorder=edgesfirst overlap=false", "rankdir = \"LR\"];",sep=" "), "graph [splines=true];") ## ----DirectionalInteractionsDF--------------------------------------------- de = data.frame( from=c("CG31687","xxx","CG31687","xxx","Gl","xxx4","xxx4","rod","abcd1","polo","Grip75", "xxx2","Bub3","xxx3","Grip84","feo","xxx5","Apc10"), to =c("xxx","cenB1A","xxx","Elongin-B","xxx4","Cdc23","vih","polo","polo","feo","xxx2", "Bub3","xxx3","Klp61F","polo","xxx5","glu","sti"), stringsAsFactors=FALSE) for (i in seq_len(nrow(de))) { out <- c(out, rep(sprintf("\"%s\" -> \"%s\" [style=invis];",de$from[i],de$to[i]),5)) } dg = c(de$from,de$to) dg = unique(dg[!(dg %in% Interactions$Anno$target$Symbol)]) out <- c(out, sprintf("%s [style=invis]\n",dg)) ## ----DirectionalInteractions-Mitosis--------------------------------------- out <- c(out, # "newrank=true;", # "1 -> 2 -> 3 -> 4 -> 5 -> 6 -> 7 -> 8 -> 9 [style=invis];", sprintf("\"%s\" -> \"%s\" [color=\"%s\" penwidth=3 arrowsize=2];", edges2$geneFrom, edges2$geneTo, edges2$color)) for (cl in seq_along(SelectedClustersComplexes)) { out = c(out,sprintf("subgraph cluster%s {",LETTERS[cl])) out = c(out, sprintf(paste("\"%s\" [label=\"%s\" shape=ellipse style=filled fillcolor=\"%s\"", "labelfontsize=%d margin=\"0.02,0.02\" tooltip=\"%s\"];"), nodes2[[cl]], nodes2[[cl]], "#F5ECE5", 20, nodes2[[cl]])) out = c(out, sprintf("label=\"%s\";",names(SelectedClustersComplexes)[cl])) out = c(out," }") } out = c(out,"}") file = file.path(resultdir,"DirectionalInteractions-Mitosis") writeLines(out, con=sprintf("%s.dot",file) ) ## ----DirectionalInteractions-Mitosis-dot, eval=FALSE----------------------- # system(sprintf("dot %s.dot -Tpdf -o%s.pdf",file,file)) ## ----DirectionalInteractions-Mitosis-plots,resize.width="0.24\\textwidth",fig.show='hold',fig.width=4,fig.height=4---- plot2Phenotypes(data,"polo","Grip84",2,3,lwd=1,cex.axis=0.5,cex.lab=0.5,cex=0.5,length = 0.1) plot2Phenotypes(data,"BubR1","Elongin-B",5,19,lwd=1,cex.axis=0.5,cex.lab=0.5,cex=0.5,length = 0.1) plot2Phenotypes(data,"rod","SMC2",1,4,lwd=1,cex.axis=0.5,cex.lab=0.5,cex=0.5,length = 0.1) plot2Phenotypes(data,"Cdc23","SMC2",1,4,lwd=1,cex.axis=0.5,cex.lab=0.5,cex=0.5,length = 0.1) plot2Phenotypes(data,"Klp61F","Cdc16",1,2,lwd=1,cex.axis=0.5,cex.lab=0.5,cex=0.5,length = 0.1) plot2Phenotypes(data,"ald","l(1)dd4",4,19,lwd=1,cex.axis=0.5,cex.lab=0.5,cex=0.5,length = 0.1) plot2Phenotypes(data,"ald","Dlic",2,19,lwd=1,cex.axis=0.5,cex.lab=0.5,cex=0.5,length = 0.1) ## ----DirectionalInteractionsExamples,fig.show='hold',resize.width="0.3\\textwidth",fig.width=10,fig.height=10,eval=FALSE---- # Genes = list("sti"=c("sti","Cdc23"), # "RasGAP1" = c("RasGAP1", "dalao", "Snr1", "osa","brm", "mor", # "Bap60", "Dsor1", "Pvr", "Sos", "pnt")) # # for (g in seq_along(Genes)) { # QG = Genes[[g]] # # edges2 = edges # edges2 = edges2[(edges2$geneFrom %in% QG) & (edges2$geneTo %in% QG),] # cat(sprintf("Writing graph with %d nodes and %d edges.\n", length(nodes2), nrow(edges2))) # edges2$color[edges2$color == "crimson"] = "#DB1D3D" # edges2$color[edges2$color == "dodgerblue"] = "#4C86C6" # edges2$width = 5 # edges2$arrow.size = 2 # # genes = data.frame(gene = unique(c(edges2$geneFrom,edges2$geneTo)),stringsAsFactors=FALSE) # genes$color = rep("#F5ECE5",nrow(genes)) # genes$frame.color = NA # genes$label.color = "black" # genes$label.cex=1.3 # genes$size=50 # g = graph.data.frame(edges2,vertices = genes) # # set.seed(3122) # plot(g) # } ## ----DirectionalInteractions-SWISNF-plots,resize.width="0.24\\textwidth",fig.show='hold',fig.width=4,fig.height=4---- plot2Phenotypes(data,"Snr1","RasGAP1",1,5,lwd=2,cex.axis=1,cex.lab=1,cex=1,length = 0.2) plot2Phenotypes(data,"Snr1","RasGAP1",1,13,lwd=2,cex.axis=1,cex.lab=1,cex=1,length = 0.2) plot2Phenotypes(data,"brm","RasGAP1",1,5,lwd=2,cex.axis=1,cex.lab=1,cex=1,length = 0.2) plot2Phenotypes(data,"brm","RasGAP1",1,13,lwd=2,cex.axis=1,cex.lab=1,cex=1,length = 0.2) plot2Phenotypes(data,"mor","RasGAP1",1,5,lwd=2,cex.axis=1,cex.lab=1,cex=1,length = 0.2) plot2Phenotypes(data,"mor","RasGAP1",1,13,lwd=2,cex.axis=1,cex.lab=1,cex=1,length = 0.2) ## ----DirectionalInteractionsLoadDI,fig.show='hold',resize.width="0.29\\textwidth",fig.width=7,fig.height=7,eval=FALSE---- # data("TID2HUGO", package="DmelSGI") # data("Interactions", package="DmelSGI") # HugoNames = sapply(TID2HUGO, function(x) { # if(length(x) > 0) { # j=0 # for (i in 1:nchar(x[1])) { # if (length(unique(substr(x,1,i))) == 1) { # j=i # } # } # res = paste(substr(x,1,j)[1],paste(substr(x,j+1,nchar(x)),collapse="/"),sep="") # } else { # res = "" # } # res # }) # HugoNames = paste(Interactions$Anno$target$Symbol," (",HugoNames,")",sep="") # names(HugoNames) = Interactions$Anno$target$Symbol # # # data("Intogen", package="DmelSGI") # # SelCancer = sapply(TID2HUGO, function(x) { any(x %in% Intogen$symbol) }) # SelCancer = Interactions$Anno$target$Symbol[which(Interactions$Anno$target$TID %in% # names(which(SelCancer)))] # # Genes = list("Pten" = c("Pten","gig"), # "Arp3" = c("Arp3","Sos"), # "Myb" = c("Myb","mip120","mip130","polo","fzy","Elongin-B","CtBP", # "sti","pav","tum","feo","Rho1","dia","scra","SMC4"), # "nonC" = c("nonC","spen"), # "Nup75" = c("Nup75","Sin3A","CtBP","jumu","RecQ4")) # # for (g in seq_along(Genes)) { # QG = Genes[[g]] # # edges2 = edges # edges2 = edges2[(edges2$geneFrom %in% QG) & (edges2$geneTo %in% QG),] # edges2 = edges2[(edges2$geneFrom %in% QG[1]) | (edges2$geneTo %in% QG[1]),] # cat(sprintf("Writing graph with %d nodes and %d edges.\n", length(nodes2), nrow(edges2))) # edges2$color[edges2$color == "crimson"] = "#DB1D3D" # edges2$color[edges2$color == "dodgerblue"] = "#4C86C6" # edges2$width = 5 # edges2$arrow.size = 2 # # genes = data.frame(gene = unique(c(edges2$geneFrom,edges2$geneTo)),stringsAsFactors=FALSE) # genes$color = rep("#DDDDDC",nrow(genes)) # genes$color[genes$gene %in% SelCancer] = "#777777" # genes$frame.color = NA # genes$label.color = "black" # genes$label.cex=1.5 # genes$size=50 # g = graph.data.frame(edges2,vertices = genes) # V(g)$name = HugoNames[V(g)$name] # # set.seed(3122) # plot(g) # legend("bottomright", inset = c(-0.07,-0.07),fill=c("#777777", "#DDDDDC"), # c("genes recurrently mutated in cancer","not recurrently mutated"),cex=0.5) # } # ## ----DmelSGIsessioninfo---------------------------------------------------- sessionInfo()