Commit d15ff486 authored by Philippe BORDRON's avatar Philippe BORDRON
Browse files

first try of interactive pca

parent 83131de9
#!/usr/bin/env RScript
library(ggplot2)
library(ggiraph)
library(grid)
library(optparse)
......@@ -49,9 +50,10 @@ ACP<-function(d,transpose=T,scale=F,center=T) {
acp2d<-function(pca, comp=1:2, group=NULL, plotVars = FALSE, pointSize=2, plotText=FALSE, fixedCoord=FALSE, main=NULL, ellipse=FALSE, color=NULL){
if(!require("ggplot2")) stop("You must install ggplot2");
if(!require("ggiraph")) stop("You must install ggiraph");
if(length(comp)!=2) stop("You must give a vector of 2 integer for comp parameter");
percentVar <- pca$percentVar
functPlot<-ifelse(plotText,geom_text,geom_point)
functPlot<-ifelse(plotText, geom_text_interactive, geom_point_interactive)
coor=ifelse(plotVars,"rotation","x")
ngroup <- 0
if(is.null(group)){
......@@ -81,7 +83,7 @@ acp2d<-function(pca, comp=1:2, group=NULL, plotVars = FALSE, pointSize=2, plotTe
override.aes = list(size=legend_text_size/(2*pointSize)) # adjust dot size into legend
))
graph <- graph + functPlot(size=pointSize)+
graph <- graph + functPlot(size=pointSize, tooltip = rownames(d)) +
xlab(paste0("PC",comp[1],": ",round(percentVar[comp[1]] * 100),"% variance")) +
ylab(paste0("PC",comp[2],": ",round(percentVar[comp[2]] * 100),"% variance"))
if(fixedCoord)graph <- graph + coord_fixed(ratio=percentVar[comp[2]]/percentVar[comp[1]])
......@@ -165,6 +167,7 @@ compo=length(which(acpT$percentVar>0.05))+1
exprDatT<-exprDatT[,c(rownames(sampleAnnot))]
pdf(file = paste(outputDir,"/NormAndPCA.pdf",sep=""),width=10,height=10)
acpTa <- NULL
if(batchEffect || plateBatch){
......@@ -197,9 +200,12 @@ if(batchEffect || plateBatch){
acp2d(acpTa,comp=c(1,2),plotText = TRUE,group = sampleAnnot[,batchEffectCol],main="PCA on normalized/transformed/adjusted data")
}
p <- acp2d(acpTa,comp=c(1,2),group = sampleAnnot[,condCol],main="PCA on normalized/transformed/adjusted data")
htmlwidgets::saveWidget(girafe(ggobj = p, pointsize = 12, width_svg = 6, height_svg = 4), paste(outputDir, "PCA.html", sep="/"))
png(paste(outputDir,"/PCA.png",sep=""),width=15,height=10,units="cm",res=300,pointsize=6)
acp2d(acpTa,comp=c(1,2),group = sampleAnnot[,condCol],main="PCA on normalized/transformed/adjusted data")
dev.off()
print(p)
dev.off()
plotAllACP(compo,acpTa,sampleAnnot,condCol)
......@@ -210,12 +216,16 @@ if(batchEffect || plateBatch){
barplot(acpT$percentVar*100,names.arg = round(acpT$percentVar*100,2),main = "Contribution of each component in the PCA")
plotAllACP(compo,acpT,sampleAnnot,condCol)
p <- acp2d(acpT,comp=c(1,2),group = sampleAnnot[,condCol],main="PCA on normalized/transformed data")
htmlwidgets::saveWidget(girafe(ggobj = p, pointsize = 12, width_svg = 6, height_svg = 4), paste(outputDir, "PCA.html", sep="/"))
png(paste(outputDir,"/PCA.png",sep=""),width=15,height=10,units="cm",res=300,pointsize=6)
acp2d(acpT,comp=c(1,2),group = sampleAnnot[,condCol],main="PCA on normalized/transformed data")
print(p)
dev.off()
}
dev.off()
save(acpT, acpTa, sampleAnnot, file=paste(outputDir, "/PCA.RData", sep=""))
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment