Skip to content

Instantly share code, notes, and snippets.

@datad

datad/DataFrames.R

Last active June 8, 2018 20:18
Show Gist options
  • Save datad/39b9401a53e7f4b44e9bea4d584ac3e8 to your computer and use it in GitHub Desktop.
Save datad/39b9401a53e7f4b44e9bea4d584ac3e8 to your computer and use it in GitHub Desktop.
Download ZIP
_utile4disSuptyper
Raw
_utile
#' version 1 (1/31/2018)
#' short description
######### source this ----
rm(list = ls())
#' Get all needed data
LoadMyData <- funtion()
{
dataFile <- "temp/1-2v3data.RData"
if( file.exists(dataFile) ){
load(dataFile, env = globalenv())
} else {
#Get w
load(file="temp/100_NMFsurv.RData", env = globalenv())
#Get A
load("results/BRCA15.RData", env = globalenv())
A <<- Boolcp_3comp$A #patients vs components
rm(list=ls()[-which(ls() %in% c("w","A"))])
libs<-c("Packages.R", "permutation.R", "clustering.R", "Plots.R")
libs<-paste("https://gist.githubusercontent.com/datad/39b9401a53e7f4b44e9bea4d584ac3e8/raw/", libs,sep='')
sapply(libs, function(u) {source(u)})
onlyBaseLibs()
loadls("SNFtool survival ROntoTools org.Hs.eg.db")
loadlib("plyr", FALSE) #for function count
loadlib("lmtest", FALSE) #for comparing two models - likelihood
loadls("survival survAUC", FALSE) #for Cox proportional hazard model
#For non-negative matrix factorization vvv
install.packages("digest", repos="http://R-Forge.R-project.org") #required
loadls("plyr survival lmtest", FALSE)
loadls("stats registry", FALSE)
loadls("methods utils pkgmaker registry rngtools cluster", FALSE)
loadlib("NMF", FALSE)
#For non-negative matrix factorization ^^^^
loadlib("glmnet") #for lasso regression L1
A <<- A
save.image(file="temp/myData.RData")
}
}
######### run this ----
LoadMyData()
#' Title
#'
#' @section step 1
#'
#' @examples
foo <- function() {
}
#end
sess <- sessionInfo() #save session on variable
save.image(file="temp/x.RData")
Raw
clustering.R
#version 2 (1/1/2018)
computeAPN <- function(pathwMolecules, cluster)
{
measuresAPN <-0
## stability validation measures
for (del in seq_along(pathwMolecules)) {
molsDel <- pathwMolecules[-del]
clusterDel <- custerByMols(molsDel)
stabmeasAPN <- stabilityAPN(cluster, clusterDel)
measuresAPN <- measuresAPN + stabmeasAPN
}
measuresAPN <- measuresAPN/length(pathwMolecules)
measuresAPN
}
stabilityAPN <- function(cluster, clusterDel) {
nc1 <- length(cluster)
nc2 <- length(clusterDel)
## measure APN
## calculate a ncxnc matrix of proportion of non-overlaps in the two collection of nc clusters
overlap <- xtabs(~cluster + clusterDel)
## measure AD
## calculate a ncxnc matrix of average-distance in the two collection of nc clusters
dij <- matrix(rep(NA,nc1*nc2),nc1,nc2)
rs <- matrix(rowSums(overlap),nrow=nrow(overlap),ncol=ncol(overlap),byrow=FALSE)
#rs <- matrix(colSums(overlap),nrow=nrow(overlap),ncol=ncol(overlap),byrow=TRUE)
APN1 <- 1-sum(overlap^2/rs)/sum(overlap)
APN1
}
#an actual pathway
clusterByPathwaySNF2 <- function(pathway.i, surData){
require(ROntoTools)
genesOfInterest <- nodes(pathway.i)
listColVals <- clusterByGenesSNF2(genesOfInterest, surData)
listColVals
}
#genesOfInterest: are genes and mirs from pathways
clusterByGenesSNF2 <- function(genesOfInterest,surData){
require(survival)
listColVals <- list()
#Number of genes on pathway
listColVals$nGenes <- length(genesOfInterest)
genesOfInterest <- geneMIRSetFormat(genesOfInterest)
groups <- getGroupsFromGenesSNF2(genesOfInterest)
listColVals$groups <- groups
stabilityAPN <- computeAPN(pathwMolecules=genesOfInterest, cluster=groups)
listColVals$stabilityAPN <- stabilityAPN
# #begin PREVIOUS METRIC
# #Kaplan SNF
# fit <- survdiff(Surv(time = Survival, event = Death) ~ groups, data = surData)
# listColVals$kaplan <- fit
#
# #Cox SNF
# coxFit <- coxph(Surv(time = Survival, event = Death) ~ groups,
# data = surData, ties="exact")
# resCox <- mySumm(coxFit)
#
# listColVals$coxPval <- resCox$sctest[3]
# listColVals$coxScore <- resCox$sctest[1]
# listColVals$resCox <- resCox
# #end PREVIOUS METRIC
listColVals
}
#Global = geneEData
#genesOfInterest: are genes from pathways
getGroupsFromGenesSNF2 <- function(genesOfInterest, displayT=FALSE){
require(SNFtool)
GEofI <- intersect(genesOfInterest,colnames(geneEData))
MIRofI <- intersect(genesOfInterest,colnames(mirEData))
genePath <- geneEData[, GEofI]
mirPath <- mirEData[, MIRofI]
#2. Calculate the pair-wise distance
PSMgeneE = dist2(as.matrix(genePath),as.matrix(genePath));
PSMmirE = dist2(as.matrix(mirPath),as.matrix(mirPath));
#3. construct similarity graphs
W1 = affinityMatrix(PSMgeneE, K, alpha)
W2 = affinityMatrix(PSMmirE, K, alpha)
## These similarity graphs have complementary information about clusters.
#displayClusters(W1,truelabel);
#displayClusters(W2,truelabel);
#next, fuse all the graphs
W = SNF(list(W1,W2), K, T)
####With this unified graph W of size n x n, you can do either spectral clustering
#or Kernel NMF. If you need help with further clustering, please let us know.
#4. clustering
groupGE = spectralClustering(W,C);
if(displayT) displayClusters(W,groupGE);
#SNFNMI = Cal_NMI(group, truelabel)
###you can also find the concordance between each individual network and the fused network
#ConcordanceMatrix = Concordance_Network_NMI(list(W, W1,W2));
groups <- factor(groupGE)
names(groups) <- row.names(geneEData)
return(groups)
}
#Global var survivalData
thetaSNF2 <- function(x){
require(survival)
indicesOfIPgene <- x[which(x<totalNumberOfGenes)]
indicesOfIPmir <- x[which(x>totalNumberOfGenes)]
indicesOfIPmir <- indicesOfIPmir - totalNumberOfGenes
genesOfInterest <- allGenes[indicesOfIPgene]
mirsOfInterest <- allMIRs[indicesOfIPmir]
molsOfInterest <- c(genesOfInterest,mirsOfInterest)
groups <- getGroupsFromGenesSNF2(molsOfInterest)
##this was before
# #Kaplan SNF
# fit <- survdiff(Surv(time = Survival, event = Death) ~ groups, data = survivalData)
# #Cox SNF
# coxFit <- coxph(Surv(time = Survival, event = Death) ~ groups, data = survivalData,
# ties="exact")
# resCox <- mySumm(coxFit)
# coxPval <- resCox$sctest[3]
# coxPval
stabilityAPN <- computeAPN(molsOfInterest,groups)
stabilityAPN
}
#Global var survivalData
thetaSNF <- function(x){
require(survival)
genesOfInterest <- allGenes[x]
groups <- getGroupsFromGenesSNF(genesOfInterest)
#Kaplan SNF
fit <- survdiff(Surv(time = Survival, event = Death) ~ groups, data = survivalData)
#Cox SNF
coxFit <- coxph(Surv(time = Survival, event = Death) ~ groups, data = survivalData,
ties="exact")
resCox <- mySumm(coxFit)
coxPval <- resCox$sctest[3]
coxPval
}
#genesOfInterest: are genes from pathways
clusterByGenesSNF <- function(genesOfInterest,surData){
require(survival)
listColVals <- list()
#Number of genes on pathway
listColVals$nGenes <- length(genesOfInterest)
genesOfInterest <- geneSetFormat(genesOfInterest)
groups <- getGroupsFromGenesSNF(genesOfInterest)
#Kaplan SNF
fit <- survdiff(Surv(time = Survival, event = Death) ~ groups, data = surData)
listColVals$kaplan <- fit
#Cox SNF
coxFit <- coxph(Surv(time = Survival, event = Death) ~ groups,
data = surData, ties="exact")
resCox <- mySumm(coxFit)
listColVals$coxPval <- resCox$sctest[3]
listColVals$coxScore <- resCox$sctest[1]
listColVals$resCox <- resCox
listColVals$groups <- groups
listColVals
}
#an actual pathway
clusterByPathwaySNF <- function(pathway.i, surData){
require(ROntoTools)
genesOfInterest <- nodes(pathway.i)
listColVals <- clusterByGenesSNF(genesOfInterest, surData)
listColVals
}
#before
# genesOfInterest <- geneSetFormat(genesOfInterest)
#Global = geneEData
#genesOfInterest: are genes from pathways
getGroupsFromGenesSNF <- function(genesOfInterest){
require(SNFtool)
genePath <- geneEData[, genesOfInterest]
##################
#2. Calculate the pair-wise distance
PSMgeneE = dist2(as.matrix(genePath),as.matrix(genePath));
#3. construct similarity graphs
W1 = affinityMatrix(PSMgeneE, K, alpha)
#Groups with just geneExpression
groupGE = spectralClustering(W1,C);
groups <- factor(groupGE)
names(groups) <- row.names(geneEData)
return(groups)
}
########
#For k-means
########
#before
# genesOfInterest <- geneSetFormat(genesOfInterest)
#Global = geneEData
#genesOfInterest: are genes from pathways
getGroupsFromGenesK <- function(genesOfInterest){
genePath <- geneEData[, genesOfInterest]
kresults <- kmeans(genePath, C, iter.max = 10000, nstart = 50)
groupGE = kresults$cluster
groups <- factor(groupGE)
return(groups)
}
#Global var survivalData
thetaK <- function(x){
genesOfInterest <- allGenes[x]
groups <- getGroupsFromGenesK(genesOfInterest)
#Kaplan SNF
fit <- survdiff(Surv(time = Survival, event = Death) ~ groups, data = survivalData)
#Cox SNF
coxFit <- coxph(Surv(time = Survival, event = Death) ~ groups, data = survivalData,
ties="exact")
resCox <- mySumm(coxFit)
coxPval <- resCox$sctest[3]
coxPval
}
#genesOfInterest: are genes from pathways
clusterByGenesK <- function(genesOfInterest,surData){
listColVals <- list()
#Number of genes on pathway
listColVals$nGenes <- length(genesOfInterest)
genesOfInterest <- geneSetFormat(genesOfInterest)
groups <- getGroupsFromGenesK(genesOfInterest)
#Kaplan SNF
fit <- survdiff(Surv(time = Survival, event = Death) ~ groups, data = surData)
listColVals$kaplan <- fit
#Cox SNF
coxFit <- coxph(Surv(time = Survival, event = Death) ~ groups,
data = surData, ties="exact")
resCox <- mySumm(coxFit)
listColVals$coxPval <- resCox$sctest[3]
listColVals$coxScore <- resCox$sctest[1]
listColVals$resCox <- resCox
listColVals$groups <- groups
listColVals
}
#an actual pathway
clusterByPathwayK <- function(pathway.i, surData){
genesOfInterest <- nodes(pathway.i)
listColVals <- clusterByGenesK(genesOfInterest, surData)
listColVals
}
Raw
DataFrames.R
#data.frame
# DF=tT
# checkTypeOfDF(DF=controlDat)
# checkTypeOfDF(DF=diseaseDat)
# checkTypeOfDF(DF)
# head(DF)
#
checkTypeOfDF <- function(DF=controlDat) {
sapply(DF, mode)
sapply(DF, class)
#IF WANT TO CHANGE
#transform(d, char = as.numeric(char))
}
mergeByrowName <- function(controlDat, diseaseDat){
allSamples <- merge(controlDat, diseaseDat, by="row.names", stringsAsFactors=FALSE )
allSamples <- allSamples [,-1]
row.names(allSamples) <- row.names(controlDat)
allSamples
}
#
# #create data.frame
# #http://www.dummies.com/how-to/content/how-to-create-a-data-frame-from-scratch-in-r.html
#
# employee <- c('John Doe','Peter Gynn','Jolie Hope')
# salary <- c(21000, 23400, 26800)
# startdate <- as.Date(c('2010-11-1','2008-3-25','2007-3-14'))
#
#
#
#
# employ.data <- data.frame(employee, salary, startdate)
#
# #create empty data frame
# #http://stackoverflow.com/questions/10689055/create-an-empty-data-frame
# df <- data.frame(Date=as.Date(character()),
# File=character(),
# User=character(),
# stringsAsFactors=FALSE)
#
# newrow = c("1/1/2010","thefile","dd")
# df <- rbind(df,newrow)
#
# #read
#
#
# head(employ.data)
#
# employ.data[,c('employee','salary')]
#
#
# employ.data[1:5,c('employee','salary')]
# employ.data[1:length(employ.data[,1]),c('employee','salary')]
#
#
#
# #add a column with a list of two columns
# #in data table
# totalCounts[,pValue:= pVHyper(freqS,freqU) ]
#
# #in dataframe
# ##http://stackoverflow.com/questions/3651651/adding-a-column-to-a-dataframe-in-r
# #https://stat.ethz.ch/pipermail/r-help/2002-April/020399.html
#
# DF <- data.frame(start=c(1,3,5,7), end=c(2,6,7,9))
# DF$newcol <- apply(DF,1,function(row) mean( row[1] : row[2] ))
# data$Z <- data$X + 5 * data$Y
#
#
#
# #Add a Column to a Dataframe From a List of Values
# #http://stackoverflow.com/questions/11130037/add-a-column-to-a-dataframe-from-a-list-of-values
#
# columnsNames = LETTERS[1:10]
# nRows = 2
#
# num.iters <- 10
# l <- vector('list', num.iters)
# for (i in 1:num.iters) {
# l[[i]] <- rep("NA",nRows) # the column data
# names(l)[i] <- columnsNames[i] # the column name
# }
# do.call(cbind, l) # ... if your cols are the same datatype and you want a matrix
# data.frame(l) # otherwise
#
#
# #DIVIDE THE DATA FRAME IN TWO ONE WITH COLUMNS EMPTY
#
# x.sub <- subset(DF, interaction == '')
#
# #get a vector for data set
# #http://stackoverflow.com/questions/2545228/converting-a-dataframe-to-a-vector-by-rows
#
# test <- data.frame(x=c(26,21,20),y=c(34,29,28))
# vectorCOlum = as.vector(as.matrix(test))
# vectorRow =as.vector(t(test))
#
# vectorColum1 = as.vector(as.matrix(test[,'x']))
# vectorRow1 = as.vector(as.matrix(test[1,]))
#
#
# #Change from factor to char
# #http://stackoverflow.com/questions/2851015/convert-data-frame-columns-from-factors-to-characters
#
# employ.data[1,1]
#
# class(employ.data[1,1])
#
# employ.data <- data.frame(lapply(employ.data, as.character), stringsAsFactors=FALSE)
#
# class(employ.data[1,1])
#
#
# #change value of column
# #http://stackoverflow.com/questions/12888363/how-to-change-column-values-in-a-data-frame
# employ.data <- within(employ.data, salary <- salary+1000)
#
# employ.data[,2] = log(employ.data[,2])
#
# #with a new function declared by us
# #http://stackoverflow.com/questions/11330138/finding-columns-with-all-missing-values-in-r
# interactionsTab[,'interaction'] <- sapply(interactionsTab[,'interaction'], interaction2weight )
#
#
#
# #delete or drop columns
# #http://stackoverflow.com/questions/4605206/drop-columns-r-data-frame
# DF <- data.frame(
# x=1:10,
# y=10:1,
# z=rep(5,10),
# a=11:20
# )
# drops <- c("x","z")
# DF[,c("a","y")]
#
#
#
# #rename
# #http://www.cookbook-r.com/Manipulating_data/Renaming_columns_in_a_data_frame/
# names(DF) [names(DF)=="a"] <- "two"
#
#
# #add a row to the dataframe
# #http://stackoverflow.com/questions/11561856/add-new-row-to-dataframe
#
# #In a r possition
#
# existingDF <- as.data.frame(matrix(seq(20),nrow=5,ncol=4))
# r <- 3
# newrow <- seq(4)
# existingDF <- rbind(existingDF[1:r,],newrow,existingDF[-(1:r),])
#
#
# #at the end
# #it is the same for adding dataframes at the end
# newrow = c(7:10)
# existingDF = rbind(existingDF,newrow)
#
#
# #transpose data frame
# #http://stackoverflow.com/questions/6778908/r-transposing-a-data-frame
#
# DF<- as.data.frame(t(DF))
#
# #delete the name of the rows
# #http://r.789695.n4.nabble.com/Remove-row-names-column-in-dataframe-td856647.html
#
# row.names(DF) <- NULL
#
# #Extract or Replace Parts of a Data Frame
# #http://stat.ethz.ch/R-manual/R-devel/library/base/html/Extract.data.frame.html
#
#
#
# #combine columns
#
# m <- cbind(1, 1:7) # the '1' (= shorter vector) is recycled
#
#
#
# #merge columns
# #https://stat.ethz.ch/pipermail/r-help/2011-June/280275.html
#
#
# prefix <- c("cheap", "budget")
# roots <- c("car insurance", "auto insurance")
# suffix <- c("quote", "quotes")
#
# df1 <- data.frame(prefix, roots, suffix)
#
# newThing <- do.call(c, df1)
#
#
#
#
# class(newThing)
#
#
# dat <- paste(df1[,2], df1[,1], sep=" ")
# head(dat)
# dat[,c(1,2)] <- NA
# head(dat)
#
#
#
#
# #removing duplicated
# #http://stats.stackexchange.com/questions/6759/removing-duplicated-rows-from-r-data-frame
#
# DF <- data.frame(
# x=1:10,
# y=10:1,
# z=rep(5,10),
# a=11:20
# )
#
# DF = rbind(DF, c(1,2,3,4))
# DF = rbind(DF, c(1,2,3,4))
# DF = rbind(DF, c(1,2,3,4))
# DF = rbind(DF, c(1,2,3,4))
#
# duplicated(DF)
#
# DF[duplicated(DF), ]
#
# DF[!duplicated(DF), ]
#
#
# #remove column
# #http://stackoverflow.com/questions/9202413/how-do-you-delete-a-column-in-data-table
#
#
# #GET NUMBER OF ROWS LENGHT SIZE IN DATA TABLE
# #http://stats.stackexchange.com/questions/5253/how-do-i-get-the-number-of-rows-of-a-data-frame-in-r
#
# nrow(dataset)
#
# dim(dataset)
#
#
#
# #delete or remove a row
# #http://stackoverflow.com/questions/13520515/command-to-remove-row-from-a-data-frame
#
# eldNew2 <- tt[-1,]
# dim(eldNew)
#
#
#
#
# #To find whether a column exists in data frame or not
# if("d" %in% colnames(dat))
# {
# cat("Yep, it's in there!\n");
# }
#
#
# #split data frame
# #http://stackoverflow.com/questions/7069076/split-column-at-delimiter-in-data-frame
#
# df <- data.frame(ID=11:13, FOO=c('a|b','b|c','x|y'))
# foo <- data.frame(do.call('rbind', strsplit(as.character(df$FOO),'|',fixed=TRUE)))
#
#
Raw
Packages.R
#version 2 (1/1/2018)
#Install libraries automatically
#Clean screen and delete elements
# rm(list=ls())
# cat("\014")
# source("../Utils/utilitary/Packages.R")
# sourceAFolder("../Utils/utilitary/")
#What is this?
#source("http://bioconductor.org/workflows.R")
# workflowInstall("arrays")
loadlib <- function (libname, bioc=TRUE, verbosel = FALSE){
#newInstruction <- try( require(libname, character.only = TRUE,
# quietly = !verbosel , warn.conflicts = verbosel), silent=TRUE )
#if( inherits(newInstruction, "try-error") )
if(!require(libname, character.only = TRUE, quietly = !verbosel ,
warn.conflicts = verbosel ) ) #previous statement #dd 05/18
{
errorMsn <- paste("Could not load package '",libname,"'.")
if(readline(prompt = paste(errorMsn,
"\nDo you like to installed it from Bioconductor? (y/n)")) == 'y')
{
source("http://bioconductor.org/biocLite.R")
biocLite(libname)
#biocLite(c("GenomicFeatures", "AnnotationDbi"))
require(libname, character.only = TRUE) || stop(errorMsn)
} else if(readline(prompt = paste(errorMsn,
"\nDo you like to installed it from CRAN? (y/n)")) == 'y'){
install.packages(libname)
require(libname, character.only = TRUE) || stop(errorMsn)
}
else{
stop(errorMsn)
}
}
if(bioc)
{
source("http://bioconductor.org/biocLite.R")
message(biocValid())
message("Above is the output of biocValid()")
}
}
loadls <- function (libs, bioc=TRUE){
libs <- unlist(strsplit(libs, " "))
lapply(libs,loadlib,bioc)
}
Raw
Packages_Extra.R
#version 2 (1/1/2018)
#Install libraries automatically
#Clean screen and delete elements
# rm(list=ls())
# cat("\014")
# source("../Utils/utilitary/Packages.R")
# sourceAFolder("../Utils/utilitary/")
#What is this?
#source("http://bioconductor.org/workflows.R")
# workflowInstall("arrays")
loadlib <- function (libname, bioc=TRUE){
if(!require(libname, character.only = TRUE))
{
errorMsn <- paste("Could not load package '",libname,"'.")
if(readline(prompt = paste(errorMsn,
"\nDo you like to installed it from Bioconductor? (y/n)")) == 'y')
{
source("http://bioconductor.org/biocLite.R")
biocLite(libname)
#biocLite(c("GenomicFeatures", "AnnotationDbi"))
require(libname, character.only = TRUE) || stop(errorMsn)
} else if(readline(prompt = paste(errorMsn,
"\nDo you like to installed it from CRAN? (y/n)")) == 'y'){
install.packages(libname)
require(libname, character.only = TRUE) || stop(errorMsn)
}
else{
stop(errorMsn)
}
}
if(bioc)
{
source("http://bioconductor.org/biocLite.R")
message(biocValid())
message("Above is the output of biocValid()")
}
}
loadls <- function (libs, bioc=TRUE){
libs <- unlist(strsplit(libs, " "))
lapply(libs,loadlib,bioc)
}
lspack <- function(libname){
print(ls(paste("package:",libname, sep='')))
}
#
listLoadedLibs <- function()
{
print((.packages()))
(.packages())
}
detachlib <- function(libname){
try( detach(paste("package:",libname, sep=''), unload=TRUE, character.only=TRUE) , silent = TRUE)
}
detachlibs <- function(libs){
libs <- unlist(strsplit(libs, " "))
lapply(libs,detachlib)
}
onlyBaseLibs <- function(){
baseLibs <- c("stats", "graphics", "grDevices", "utils", "datasets", "methods", "base")
ldl <- listLoadedLibs()
detachlibs(ldl[-which(ldl %in% baseLibs)])
}
searchMethods <- function(methodName) {
loadlib("sos")
findFn(methodName, maxPages=2, sortby="MaxScore")
}
sourceAFolder <- function(folder) {
scripts <- dir(folder, pattern = "*.R")
scripts <- paste(folder,scripts,sep="/")
sapply(scripts,source)
}
printRVersion <- function(){
R.Version()
getRversion()
}
#update bioconductor after instar new R and new Rstudio
updateBiocon <- function(){
remove.packages("BiocInstaller")
source("https://bioconductor.org/biocLite.R")
biocLite("BiocInstaller")
biocLite()
}
updateAllCRANPackages <- function() {
update.packages()
inst <- packageStatus()$inst
inst[inst$Status != "ok", c("Package", "Version", "Status")]
loadedPackages <- (.packages())
}
#########
#
unistallPackages <- function(packageNames){
#Option B: remove.packages2() (function from reposTools).
libs <- unlist(strsplit(packageNames, " "))
lapply(libs,function(x){remove.packages(x)})
}
#
reinstallPackages <- function(thePackages){
unistallPackages(thePackages)
loadls(thePackages)
}
#
UpgradeR <- function(){
#http://cran.r-project.org/bin/windows/base/rw-FAQ.html "What's the best way to upgrade?"
#TODO copy (say) R\win-library\3.0 to R\win-library\3.1
update.packages(checkBuilt=TRUE, ask=FALSE)
}
#
# whereLibs <- function(){
# #to figure out where your packages are located
# .libPaths()
# }
#
#
#
# resetAllVariables <- function(){
# is.connection <- function(x) inherits(x, "connection")
#
# get_connections <- function(envir = parent.frame())
# {
# Filter(is.connection, mget(ls(envir = envir), envir = envir))
# }
#
# close_all_connections <- function()
# {
# lapply(c(sys.frames(), globalenv(), baseenv()),
# function(e) lapply(get_connections(e), close))
# }
#
# close_all_connections()
#
# reset_options <- function()
# {
# is_win <- .Platform$OS.type == "windows"
# options(
# add.smooth = TRUE,
# browserNLdisabled = FALSE,
# CBoundsCheck = FALSE,
# check.bounds = FALSE,
# continue = "+ ",
# contrasts = c(
# unordered = "contr.treatment",
# ordered = "contr.poly"
# ),
# defaultPackages = c(
# "datasets",
# "utils",
# "grDevices",
# "graphics",
# "stats",
# "methods"
# ),
# demo.ask = "default",
# device = if(is_win) windows else x11,
# device.ask.default = FALSE,
# digits = 7,
# echo = TRUE,
# editor = "internal",
# encoding = "native.enc",
# example.ask = "default",
# expressions = 5000,
# help.search.types = c("vignette", "demo", "help"),
# help.try.all.packages = FALSE,
# help_type = "text",
# HTTPUserAgent = with(
# R.version,
# paste0(
# "R (",
# paste(major, minor, sep = "."),
# " ",
# platform,
# " ",
# arch,
# " ",
# os,
# ")"
# )
# ),
# internet.info = 2,
# keep.source = TRUE,
# keep.source.pkgs = FALSE,
# locatorBell = TRUE,
# mailer = "mailto",
# max.print = 99999,
# menu.graphics = TRUE,
# na.action = "na.omit",
# nwarnings = 50,
# OutDec = ".",
# pager = "internal",
# papersize = "a4",
# pdfviewer = file.path(R.home("bin"), "open.exe"),
# pkgType = if(is_win) "win.binary" else "source",
# prompt = "> ",
# repos = c(
# CRAN = "@CRAN@",
# CRANextra = "http://www.stats.ox.ac.uk/pub/RWin"
# ),
# scipen = 0,
# show.coef.Pvalues = TRUE,
# show.error.messages = TRUE,
# show.signif.stars = TRUE,
# str = list(
# strict.width = "no",
# digits.d = 3,
# vec.len = 4
# ),
# str.dendrogram.last = "`",
# stringsAsFactors = TRUE,
# timeout = 60,
# ts.eps = 1e-05,
# ts.S.compat = FALSE,
# unzip = "internal",
# useFancyQuotes = TRUE,
# verbose = FALSE,
# warn = 0,
# warning.length = 1000,
# width = 80,
# windowsTimeouts = c(100, 500)
# )
# )
#
#
#
# source(file.path(R.home("etc"), "Rprofile.site"))
#
# cleanHistory()
# }
# reset_options()
# }
#
#
Raw
permutation.R
# #Example
# results <- resampling(x=1:totalNumberOfGenes,
# sampleSize=NgenesIp,
# nboot=Nrp, theta)
resampling <- function (x, sampleSize, nboot, theta)
{
#with no replacement replace = FALSE
resampled <- matrix(t(replicate(nboot, sample(x,
size = sampleSize,
replace = FALSE))),
nrow = nboot)
#TESTS
#expect TRUE:
#all(apply(bootsam, 1, function(x){ length(unique(x)) == dim(bootsam)[2]}))
#expect FALSE:
#all(apply(bootsam, 2, function(x){ length(unique(x)) == dim(bootsam)[1]}))
thetastar <- apply(resampled, 1, theta)
return(list(thetastar = thetastar, theRandomData = resampled))
}
Raw
PlotPathways.R
#version 2 (1/1/2018)
#loadls("igraph")
###########
#This method plots a graph from a graphnel with the nodes on entrez IDs (like before send to ROntoTools)
# layout.auto(graph, dim=2, ...)
# layout.random(graph, params, dim=2)
# layout.circle(graph, params)
# layout.sphere(graph, params)
# layout.fruchterman.reingold(graph, ..., dim=2, params)
# layout.kamada.kawai(graph, ..., dim=2, params)
# layout.spring(graph, ..., params)
# layout.reingold.tilford(graph, ..., params)
# layout.fruchterman.reingold.grid(graph, ..., params)
# layout.lgl(graph, ..., params)
# layout.graphopt(graph, ..., params=list())
# layout.svd(graph, d=shortest.paths(graph), ...)
# layout.norm(layout, xmin = NULL, xmax = NULL, ymin = NULL, ymax = NULL,
# zmin = NULL, zmax = NULL)
# layoutP = layout.reingold.tilford
# layoutP = layout.fruchterman.reingold
# layoutP = layout.circle
# layoutP = layout.reingold.tilford (circular =T)
# layoutP = layout.reingold.tilford (pathway.igraph, circular =T)
# layoutP = layout.reingold.tilford (pathway.igraph, niter=10000)
plotPathwayIgraph <- function(pathway.i, name = "A pathway", entrezOrganism = "org.Hs.eg",
layoutP = layout.auto(pathway.i),
shape="rectangle", vertex.color="lightblue", edge.color="black"){
loadlib(paste(entrezOrganism,".db",sep=''))
entrez2Sym <- get(paste(entrezOrganism,"SYMBOL",sep=''))
pathway.igraph <- igraph.from.graphNEL(pathway.i)
entrezID <- V(pathway.igraph)$name
entrezID <- sub("hsa:","", entrezID)
symbols <-lapply(entrezID, function(x){get(x, entrez2Sym)} )
symbols <- unlist(symbols)
V(pathway.igraph)$name <- symbols
E(pathway.igraph)$color = edge.color
# layoutP = layout.reingold.tilford (pathway.igraph, 1,)
#
# layoutP = layout.norm(layoutP, xmin = 1, xmax = 600, ymin = 1, ymax = 600)
#
# layoutP = adjustLayout(layoutP, mindist = 80)
#
# layoutP = layout.norm(layoutP, xmin = 1, xmax = 600, ymin = 1, ymax = 600)
#
plot( pathway.igraph,
vertex.shape=shape,layout=layoutP,
vertex.color=vertex.color,
vertex.size =(nchar(symbols) * 20),
vertex.size2 = (nchar(symbols) * 5 ),
vertex.label.dist=0.1,
edge.arrow.size=0.5, main=name)
}
#######
# Example: plotPathway(kegg_pathways$"path:hsa04122", name = " ", vertex.color="cadetblue2")
# select the color from http://www.stat.columbia.edu/~tzheng/files/Rcolor.pdf
plotPathway <- function(pathway.i, name = "A pathway", entrezOrganism = "org.Hs.eg",
shape="box", vertex.color="lightblue", edge.color="black", fontsize = 10, showEdgeData = TRUE){
loadlib(paste(entrezOrganism,".db",sep=''), bioc = FALSE)
loadlib("ROntoTools", bioc = FALSE)
par(mgp=c(1.5,0.5,0), #axis title, axis labels and axis line
mai=c(0.6,0.6,0.02,0.02) #c(bottom, left, top, right)
)
entrez2Sym <- get(paste(entrezOrganism,"SYMBOL",sep=''))
entrezID <- nodes(pathway.i)
entrezID <- sub("hsa:","", entrezID)
symbols <-lapply(entrezID, function(x){get(x, entrez2Sym)} )
symbols <- unlist(symbols)
nodes(pathway.i) <- symbols
edgeLabels <- unlist(edgeWeights(pathway.i))
edgesFromTo <- names(edgeLabels)
if(showEdgeData){
edgeLabels <- unlist(edgeData(pathway.i))
}
nA = makeNodeAttrs(pathway.i, fixedSize=FALSE, height = "1",
width = "1", fillcolor = vertex.color,
shape = shape, fontsize = fontsize)
eAttrs <- list()
edgesFromTo <- gsub("\\.", "~", edgesFromTo)
edgesFromTo -> names(edgeLabels)
eAttrs$label <- edgeLabels
att = list(graph = list(rankdir = "LR", rank = ""))
plot(pathway.i, main = name, attrs = att, nodeAttrs = nA, edgeAttrs = eAttrs)
}
#plotTheNelGraph:plot Using graphNEL
#Params
#gR: the graph
#factorsE: the factors
plotPathway2Colors <- function(pathway.i, subclass, name = "A pathway", entrezOrganism = "org.Hs.eg",
twocolors = c("lightgreen", "black"), twoshapes = c("box", "box"),
twofontcolor = c("black", "white") , fontsize = 10, numberofchar = 6, showEdgeData = TRUE){
loadlib(paste(entrezOrganism,".db",sep=''))
entrez2Sym <- get(paste(entrezOrganism,"SYMBOL",sep=''))
entrezID <- nodes(pathway.i)
entrezID <- sub("hsa:","", entrezID)
symbols <-lapply(entrezID, function(x){if(exists(x, entrez2Sym)) get(x, entrez2Sym) else x} )
symbols <- substr(symbols, 1, numberofchar)
symbols <- unlist(symbols)
nodes(pathway.i) <- symbols
edgeLabels <- unlist(edgeWeights(pathway.i))
edgesFromTo <- names(edgeLabels)
if(showEdgeData){
edgeLabels <- unlist(edgeData(pathway.i))
}
nodeType <- 1 + (nodes(pathway.i) %in% subclass)
nA = makeNodeAttrs(pathway.i,fixedSize=FALSE,
height = "1", width = "1",
fillcolor = twocolors[nodeType], shape = twoshapes[nodeType],
fontcolor = twofontcolor[nodeType], fontsize = fontsize )
eAttrs <- list()
edgesFromTo <- gsub("\\.", "~", edgesFromTo)
edgesFromTo -> names(edgeLabels)
eAttrs$label <- edgeLabels
att = list(graph = list(rankdir = "LR", rank = ""))
plot(pathway.i, main = name, attrs = att, nodeAttrs = nA, edgeAttrs = eAttrs)
legend("bottomright", legend = c("Original", "New"),
col = twocolors,
pch = c(19,19), title = "Box color",
#text.col = c("black","lightblue"), # terrain.colors(10)[1:2]
#lwd = 2,
#cex = 0.5,
text.font = 3
)
}
plotPathway2ColorsIGrpah <- function(pathway.i, subclass, name = "A pathway", entrezOrganism = "org.Hs.eg", layoutP = layout.auto, shape="rectangle",
vertex.color="lightblue", vertex.color.subClass="red", edge.color="black"){
loadlib(paste(entrezOrganism,".db",sep=''))
entrez2Sym <- get(paste(entrezOrganism,"SYMBOL",sep=''))
pathway.igraph <- igraph.from.graphNEL(pathway.i)
entrezID <- V(pathway.igraph)$name
entrezID <- sub("hsa:","", entrezID)
symbols <-lapply(entrezID, function(x){if(exists(x, entrez2Sym)) get(x, entrez2Sym) else x} )
symbols <- unlist(symbols)
V(pathway.igraph)$name <- symbols
vertexColors <- sapply(symbols %in% subclass, function(x) {if(x) vertex.color else vertex.color.subClass})
E(pathway.igraph)$color = edge.color
plot( pathway.igraph, vertex.shape=shape,layout=layoutP,
vertex.color=vertexColors,
vertex.size = (nchar(symbols) * 20),
vertex.size2 = (nchar(symbols) * 5 ),
vertex.label.dist=0.1, edge.arrow.size=0.5, main=name)
}
#tkplot
#http://stackoverflow.com/questions/17627052/change-background-color-of-tkplot-igraph-r
tkplotplotPathway <- function(pathway.igraph){
loadlib(paste(entrezOrganism,".db",sep=''))
entrez2Sym <- get(paste(entrezOrganism,"SYMBOL",sep=''))
pathway.igraph <- igraph.from.graphNEL(pathway.i)
entrezID <- V(pathway.igraph)$name
entrezID <- sub("hsa:","", entrezID)
symbols <-lapply(entrezID, function(x){get(x, entrez2Sym)} )
symbols <- unlist(symbols)
V(pathway.igraph)$name <- symbols
E(pathway.igraph)$color = edge.color
# layoutP = layout.reingold.tilford (pathway.igraph, 1,)
#
# layoutP = layout.norm(layoutP, xmin = 1, xmax = 600, ymin = 1, ymax = 600)
#
# layoutP = adjustLayout(layoutP, mindist = 80)
#
# layoutP = layout.norm(layoutP, xmin = 1, xmax = 600, ymin = 1, ymax = 600)
#
temporalT <- tkplot( pathway.igraph,
vertex.shape=shape,layout=layoutP,
vertex.color=vertex.color,
vertex.size =(nchar(symbols) * 20),
vertex.size2 = (nchar(symbols) * 5 ),
vertex.label.dist=0.1,
edge.arrow.size=0.5, main=name)
tkconfigure(igraph:::.tkplot.get(temporalT)$canvas, "bg"="white")
tkplot( pathway.igraph,
vertex.label = V(pathway.igraph)$name,
vertex.label.color = vertex.color,
vertex.frame.color = vertex.color,
#vertex.shape=shape,
vertex.size =(nchar(symbols) * 5),
#vertex.size2 = (nchar(symbols) * 5 ),
vertex.label.dist=0.1,
edge.arrow.size=0.5,
edge.width = 0.2,
main=name,
edge.curved = T,
layout=layout.reingold.tilford)
}
rglplot <- function(g){
rglplot(g)
}
###############
#GRAPHNEL
#More Examples http://www.bioconductor.org/packages/release/bioc/vignettes/graph/inst/doc/graphAttributes.R
simpleTestGRAPHNEL <- function(){
library(graph)
set.seed(123)
V <- LETTERS[1:4]
edL <- vector("list", length=4)
names(edL) <- V
for(i in 1:4)
edL[[i]] <- list(edges=5-i, weights=runif(1))
gR <- new("graphNEL", nodes=V, edgeL=edL)
edges(gR)
edgeWeights(gR)
head(edL)
edL[[1]] <- list(edges=2)
edL$'C' <- list(edges='D')
gR <- new("graphNEL", nodes=V, edgeL=edL, edgemode="directed")
plot(gR)
V <- c(V,"z")
V
plot(gR)
gX <- addEdge("A", "C", gR, 1)
plot(gX)
#count adjacent nodes
nodes(gR)
adj(gR, "A")
set.seed(123)
g1 = randomEGraph(LETTERS[1:15], edges=100)
g1
plot(g1)
}
#Rgraphviz
#other graphs
#http://students.washington.edu/mclarkso/documents/gplot%20Ver2.pdf
#labels
#http://mae.ucdavis.edu/dsouza/Lectures/Rgraphviz.pdf
simpleTestRgraphviz <- function(){
library("Rgraphviz")
defAttrs <- getDefaultAttrs()
set.seed(123)
V <- letters[1:10]
M <- 1:4
g1 <- randomGraph(V, M, 0.2)
plot(g1, attrs = list(node = list(label = "foo", fillcolor = "lightgreen"),
edge = list(color = "cyan"), graph = list(rankdir = "LR")))
nAttrs <- list()
eAttrs <- list()
z <- strsplit(packageDescription("Rgraphviz")$Description, " ")[[1]]
z <- z[1:numNodes(g1)]
names(z) = nodes(g1)
nAttrs$label <- z
eAttrs$label <- c("a~h" = "Label 1", "c~h" = "Label 2")
attrs <- list(node = list(shape = "ellipse", fixedsize = FALSE))
plot(g1, nodeAttrs = nAttrs, edgeAttrs = eAttrs, attrs = attrs)
eAttrs <- list()
ew <- as.character(unlist(edgeWeights(g1)))
ew <- ew[setdiff(seq(along = ew), removedEdges(g1))]
names(ew) <- edgeNames(g1)
eAttrs$label <- ew
plot(g1, nodeAttrs = nAttrs, edgeAttrs = eAttrs, attrs = attrs)
nAttrs$color <- c(a = "red", b = "red", g = "green", d = "blue")
eAttrs$color <- c("a~d" = "blue", "c~h" = "purple")
nAttrs$fillcolor <- c(j = "yellow")
nAttrs$fontcolor <- c(e = "green", f = "red")
eAttrs$fontcolor <- c("a~h" = "green", "c~h" = "brown")
nAttrs
plot(g1, nodeAttrs = nAttrs, attrs = attrs)
attrs$node$shape <- "ellipse"
nAttrs$shape <- c(g = "box", f = "circle", j = "box", a = "plaintext")
nodes <- buildNodeList(g1)
edges <- buildEdgeList(g1)
nodes[[1]]
nodes <- buildNodeList(g1, nodeAttrs = nAttrs, defAttrs = defAttrs$node)
edges <- buildEdgeList(g1, edgeAttrs = eAttrs, defAttrs = defAttrs$edge)
nodes[[1]]
for (j in c("a~e", "a~h")) edges[[j]]@attrs$arrowhead <- "open"
vv <- agopen(name = "foo", nodes = nodes, edges = edges, attrs = attrs,
edgeMode = "undirected")
plot(vv)
}
plotsExamples <-function(){
data(graphExamples)
z <- graphExamples[[8]]
nNodes <- length(nodes(z))
nA <- list()
nA$fixedSize <- rep(FALSE, nNodes)
nA$height <- nA$width <- rep("1", nNodes)
nA$label <- rep("Z", nNodes)
nA$color <- rep("green", nNodes)
nA$fillcolor <- rep("orange", nNodes)
nA$shape <- rep("box", nNodes)
nA$fontcolor <- rep("blue", nNodes)
nA$fontsize <- rep(40, nNodes)
nA <- lapply(nA, function(x) {
names(x) <- nodes(z)
x
})
plot(z, nodeAttrs = nA)
}
bipartitegraphs<-function(){
set.seed(123)
nodes1 <- paste(0:7)
nodes2 <- letters[1:10]
ft <- cbind(sample(nodes1, 24, replace = TRUE), sample(nodes2,24, replace = TRUE))
ft <- ft[!duplicated(apply(ft, 1, paste, collapse = "")), ]
g <- ftM2graphNEL(ft, edgemode = "directed")
g
twocolors <- c("#D9EF8B", "#E0F3F8")
nodeType <- 1 + (nodes(g) %in% nodes1)
nA = makeNodeAttrs(g, fillcolor = twocolors[nodeType])
sg1 = subGraph(nodes1, g)
sgL = list(list(graph = sg1, cluster = FALSE, attrs = c(rank = "sink")))
att = list(graph = list(rankdir = "LR", rank = ""))
plot(g, attrs = att, nodeAttrs = nA, subGList = sgL)
}
#http://svitsrv25.epfl.ch/R-doc/library/graph/html/inEdges.html
getedgesincoming<-function(){
V <- LETTERS[1:4]
edL3 <- vector("list", length=4)
for(i in 1:4)
edL3[[i]] <- list(edges=(i%%4)+1, weights=i)
names(edL3) <- V
gR3 <- new("graphNEL", nodes=V, edgeL=edL3, "directed")
inEdges(c("A", "B"), gR3)
plot(gR3)
inEdges("A", gR3)
}
#FOR PAINTING NICE GRAPHS USE IGRAPH
#http://igraph.sourceforge.net/doc/R/graphNEL.html
simpleTestGRAPHNEL <- function(){
library (igraph)
plotPathway <- function(pathway.i){
pathway.i <- igraph.from.graphNEL(pathway.i)
shapes1 <- rep("square", length(E(pathway.i)))
shapes2 <- rep("cicle", length(E(pathway.i)))
theColors <- c("red", "black", "green")
twoshapes <- c("circle", "square")
nodeType <- 1 + (nodes(gR) %in% factorsE[,"factor"])
edgeType <- 2 + (E(pathway.i)$weight )
shapes = sapply(nodeType,function(x) twoshapes[x])
colors = sapply(nodeType,function(x) theColors[x])
E(pathway.i)$color = sapply(edgeType,function(x) theColors[x])
E(pathway.i)$weight = sapply(E(pathway.i)$weight,function(x) {if(x==-1 | x==1) 2 else 1 })
if(tkplotV){
temporalT <- tkplot(pathway.i, vertex.shape=shapes,layout=layoutP,
vertex.color=colors, vertex.size=4, edge.width=E(pathway.i)$weight,
edge.arrow.size=0.5)
tkconfigure(igraph:::.tkplot.get(temporalT)$canvas, "bg"="white")
}
else{
plot( pathway.i, vertex.shape=shapes,layout=layoutP,
vertex.color=colors, vertex.size=4, edge.width=E(pathway.i)$weight,
vertex.label.dist=0.3, edge.arrow.size=0.5) #edge.label=E(pathway.i)$weight,
}
}
#SOME PAINTS
#http://igraph.sourceforge.net/doc/R/plot.common.html
g <- graph.ring(10)
# plotting a random graph, set the parameters in the command arguments
g <- barabasi.game(100)
plot(g, layout=layout.fruchterman.reingold, vertex.size=4,
vertex.label.dist=0.5, vertex.color="red", edge.arrow.size=0.5)
# plot a random graph, different color for each component
g <- erdos.renyi.game(100, 1/100)
comps <- clusters(g)$membership
colbar <- rainbow(max(comps)+1)
V(g)$color <- colbar[comps+1]
plot(g, layout=layout.fruchterman.reingold, vertex.size=5, vertex.label=NA)
# plot communities in a graph
g <- graph.full(5) %du% graph.full(5) %du% graph.full(5)
g <- add.edges(g, c(0,5, 0,10, 5,10))
com <- spinglass.community(g, spins=5)
V(g)$color <- com$membership+1
g <- set.graph.attribute(g, "layout", layout.kamada.kawai(g))
plot(g, vertex.label.dist=1.5)
# draw a bunch of trees, fix layout
igraph.options(plot.layout=layout.reingold.tilford)
plot(graph.tree(20, 2))
plot(graph.tree(50, 3), vertex.size=3, vertex.label=NA)
tkplot(graph.tree(50, 2, mode="undirected"), vertex.size=10,
vertex.color="green")
#draw the edges
#https://sites.google.com/site/daishizuka/toolkits/sna/weighted-edges
plot.igraph(net,vertex.label=V(net)$name,layout=layout.fruchterman.reingold, edge.color="black",edge.width=E(net)$weight)
#name of the weigth
#http://cran.r-project.org/web/packages/igraph/igraph.pdf
plot( g3, layout=layout.circle, edge.label=E(g3)$weight)
#bipartite
g <- graph.ring(10)
bipartite.mapping(g)
## A star is fine, too
g2 <- graph.star(10)
bipartite.mapping(g2)
## A graph containing a triangle is not fine
g3 <- graph.ring(10)
g3 <- add.edges(g3, c(1,3))
bipartite.mapping(g3)
plot(g)
##vertex shape
#http://igraph.sourceforge.net/doc/R/igraph.vertex.shapes.html
# all vertex shapes, minus "raster", that might not be available
#http://igraph.sourceforge.net/screenshots2.html
shapes <- setdiff(vertex.shapes(), "")
g <- graph.ring(length(shapes))
set.seed(42)
plot(g, vertex.shape=shapes, vertex.label=shapes, vertex.label.dist=1,
vertex.size=15, vertex.size2=15,
vertex.pie=lapply(shapes, function(x) if (x=="pie") 2:6 else 0),
vertex.pie.color=list(heat.colors(5)), edge.color="black",edge.width=2)
# add new vertex shape, plot nothing with no clipping
add.vertex.shape("nil")
plot(g, vertex.shape="nil")
#R get edge weigth value
#http://igraph.sourceforge.net/doc/R/attributes.html
g <- graph.ring(10)
g <- set.graph.attribute(g, "name", "RING")
# It is the same as
g$name <- "RING"
g$name
g <- set.vertex.attribute(g, "color", value=c("red", "green"))
get.vertex.attribute(g, "color")
g <- set.edge.attribute(g, "weight", value=runif(ecount(g)))
get.edge.attribute(g, "weight")
E(g)
list.edge.attributes(g)
get.edge.attribute(g, "weight", index=1)
# The following notation is more convenient
g <- graph.star(LETTERS[1:9])
V(g)$color <- c("red", "green")
V(g)$color
E(g)$weight <- runif(ecount(g))
E(g)$weight
print(g, g=TRUE, v=TRUE, e=TRUE)
plot(g2)
#neighbors incoming
neighbors(g2, "B", mode="in")
E(g2)
V(g2)
incident(g, 1)
theName <- V(g2)[1]
are.connected(g2, theName, "B")
are.connected(g, 2, 4)
get.edges(g, 1)
get.edge(g, 4)
get.vertex(2)
ei <- get.edge.ids(g, c(1,2, 4,5))
#
library(igraph)
camp <- graph.formula(Harry:Steve:Don:Bert - Harry:Steve:Don:Bert,
Pam:Brazey:Carol:Pat - Pam:Brazey:Carol:Pat,
Holly - Carol:Pat:Pam:Jennie:Bill,
Bill - Pauline:Michael:Lee:Holly,
Pauline - Bill:Jennie:Ann,
Jennie - Holly:Michael:Lee:Ann:Pauline,
Michael - Bill:Jennie:Ann:Lee:John,
Ann - Michael:Jennie:Pauline,
Lee - Michael:Bill:Jennie,
Gery - Pat:Steve:Russ:John,
Russ - Steve:Bert:Gery:John,
John - Gery:Russ:Michael)
V(camp)$label <- V(camp)$name
set.seed(42) ## to make this reproducable
co <- layout.auto(camp)
plot(0, type="n", ann=FALSE, axes=FALSE, xlim=extendrange(co[,1]),
ylim=extendrange(co[,2]))
plot(camp, layout=co, rescale=FALSE, add=TRUE,
vertex.shape="rectangle",
vertex.size=(strwidth(V(camp)$label) + strwidth("oo")) * 100,
vertex.size2=strheight("I") * 2 * 100)
#nice plots
#http://rulesofreason.wordpress.com/tag/tkplot/
#http://cran.r-project.org/web/packages/igraph/igraph.pdf
library(igraph)
L3 <- LETTERS[1:8]
d <- data.frame(start = sample(L3, 16, replace = T), end = sample(L3, 16, replace = T),
weight = c(20,40,20,30,50,60,20,30,20,40,20,30,50,60,20,30))
g <- graph.data.frame(d, directed = T)
V(g)$name
E(g)$weight
ideg <- degree(g, mode = "in", loops = F)
col=rainbow(12) # For edge colors
plot.igraph(g,
vertex.label = V(g)$name, vertex.label.color = "gray20",
vertex.size = ideg*25 + 40, vertex.size2 = 30,
vertex.color = "gray90", vertex.frame.color = "gray20",
vertex.shape = "rectangle",
edge.arrow.size=0.5, edge.color=col, edge.width = E(g)$weight / 10,
edge.curved = T,
layout = layout.reingold.tilford)
###############################
#http://stackoverflow.com/questions/12088473/spacing-vertices-evenly-in-igraph-in-r
set.seed(1410)
df<-data.frame(
"site.x"=c(rep("a",3),rep("b",3),rep("c",3),rep("d",3)),
"site.y"=c(rep(c("e","f","g"),4)),
"bond.strength"=sample(1:100,12, replace=TRUE))
library(igraph)
df<-graph.data.frame(df)
V(df)$names <- c("a","b","c","d","e","f","g")
layOUT<-data.frame(x=c(rep(1,4),rep(2,3)),y=c(4:1,3:1))
layOUT[1,2] <- 3
E(df)[ bond.strength < 101 ]$color <- "red"
E(df)[ bond.strength < 67 ]$color <- "yellow"
E(df)[ bond.strength < 34 ]$color <- "green"
V(df)$color <- "white"
l<-as.matrix(layOUT)
plot(df,layout=l,vertex.size=10,vertex.label=V(df)$names,
edge.arrow.size=0.01,vertex.label.color = "black")
################################
#Other
#http://stackoverflow.com/questions/7521381/draw-network-in-r-control-edge-thickness-plus-non-overlapping-edges
Edges <- data.frame(
from = rep(c("asdsf","iorejoriejo","asfasfas","asdfsdf","fasdfsawwwwww"),each=5),
to = rep(c("asdsf","iorejoriejo","asfasfas","asdfsdf","fasdfsawwwwww"),times=5),
thickness = abs(rnorm(25)))
Edges <- subset(Edges,from!=to)
qgraph(Edges,esize=5,gray=TRUE)
}
#this looks nice http://www.vesnam.com/Rblog/page/2/
#http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3702256/
Raw
Plots.R
#version 2 (1/1/2018)
myDensity <-function(fileNam, ptitle)
{
#myDensityDes(fileNam, ptitle)
myDensityHis(fileNam, ptitle)
}
myDensityHis <-function(fileNam,ptitle)
{
font_size_times = 2.5 #big plot
pdf(file = fileNam , onefile = TRUE, pagecentre = FALSE, compress = FALSE)
par(mgp=c(2.5,0.8,0),
mai=c(1.2,1.5,1.0,0.5) ) #c(bottom, left, top, right)
hist(hx, freq=FALSE, breaks = 50, main="", col = "darkseagreen1" ,
lty = 1, #line type solid
cex.lab=font_size_times, #font size
xlab = "Cox p-value",
cex.sub = font_size_times,
cex.axis = font_size_times
)
abline(v=FinalPCP,col="red", lwd=3 )
legend("topright",
paste("p-value=",FinalPCP),
#text.font=2, #2 bold
cex=font_size_times-0.5,
box.lty=0
)
title(ptitle, line = 0.5,
adj=0.4, #centers the titles (0 left - 1 right)
cex.main=(font_size_times+0.3))
dev.off()
}
myDensityDes <-function(fileNam, ptitle)
{
pdf(file = fileNam , onefile = TRUE)
denshx <- density(
hx, kernel = "gaussian", from = 0, to = 1, bw = "SJ" )
plot(denshx, main = ptitle)
abline(v = FinalPCP,col = 4)
abline(v = crcPthpval,col = 4)
legend("topright", paste( "p-val=",round(FinalPCP,5),
"\n Cox p-val=", round(crcPthpval,5) ))
dev.off()
}
########
#plot.survfit manual
#https://stat.ethz.ch/R-manual/R-devel/library/survival/html/plot.survfit.html
# @param groups is a factor
# @param mainTitle is a string
# @param survivalData is data.frame with headers "PatientID" "Survival" "Death"
#> head(survivalGBM)
# PatientID Survival Death
#1 TCGA-02-0001-01C-01T-0179-07 358 1
#2 TCGA-02-0003-01A-01T-0179-07 144 1
#3 TCGA-02-0007-01A-01T-0179-07 705 1
#
# Example 1:
# pdf( file=paste(resultsFolder, "SNF_all_kaplan.pdf", sep="") , onefile=TRUE)
# plot.SurvivalK3 (groups, mainTitle, survivalData, "bottomright")
# dev.off()
#
#
# Example 2:
# mainTitle <- "b) Survival curve, SNF, selected genes"
# pdf(
# file = paste(resultsFolder, "SNF_union_kaplan.pdf", sep = ""), onefile =
# TRUE, width = 9, height = 7
# )
# plot.SurvivalK3 (groups, mainTitle, survivalData,"bottomright",
# colorsP = c("red","blue","black"))
# dev.off()
#
#
# Example 3:
# labelsGroups <- c(paste("group 1 (",table(groups)[2], ")", sep="" ),
# paste("group 2 (",table(groups)[1], ")", sep="" ),
# paste("group 3 (",table(groups)[3], ")", sep="" ) )
# colorsLines <- c("black","red","blue")
# colorsLabels <- c("red","black","blue")
# pdf( file = pdfName, onefile = TRUE, width = 9, height = 7)
# plot.SurvivalK3 (groups, mainTitle, survivalData,"bottomright",
# labelClu = labelsGroups,
# colorsP = colorsLines,
# colorsL = colorsLabels,
# centerT = centerT )
# dev.off()
#
#
#TRY THIS: http://www.r-statistics.com/2013/07/creating-good-looking-survival-curves-the-ggsurv-function/
plot.SurvivalK3 <- function(groups, mainTitle, survData,llocation,
colorsP = c("red","black","blue"),
colorsL = c("red","black","blue"),
labelClu = c("group1","group2", "group3"),
centerT = 0.4){
font_size_times = 1.5
require("survival")
par( mgp=c(1.8,0.3,0), #axis title, axis labels and axis line
mai=c(2.0,2.0,1.0,1.0) ) #c(bottom, left, top, right)
coxFit <- coxph(Surv(time = Survival, event = Death) ~ groups,
data = survData, ties="exact")
resCox <- mySumm(coxFit)
mfit <- survfit(Surv(time = Survival, event = Death) ~ groups, data = survData)
pVal<-format.pval(resCox$sctest[3], digits = 3)
plot(mfit, col=colorsP,
lty = 1, #line type solid
lwd=5, #line tightness
#main = mainTitle,
cex.lab=font_size_times, #font size
#xaxt="n", #do not plot x axe
xlab = "Years",
conf.int=F,
#cex.main=(font_size_times+0.1),
cex.sub =font_size_times,
cex.axis =font_size_times,
mark=3, #type of sensored mark http://www.statmethods.net/advgraphs/parameters.html
cex = 2, #size of mark
xscale = 365.25, #values in years
ylab="Survival probability",
xaxs="i", #survival curve should touch the y-axis
tck=0.01 #length of tick marks
)
title(mainTitle, line = 0.5,
adj=centerT, #centers the titles (0 left - 1 right)
cex.main=(font_size_times+0.3))
# axis(1, at = seq(0, max(survData$Survival), length.out = 15),
# labels=round(seq(0, 4.1, length.out = 15), digits=1), las=1) #customize x axe
legend(llocation,
labelClu,
col=colorsL,
text.col=colorsL,
text.font=2, #bold
fill=colorsL, #boxes
cex=font_size_times,
title = paste(" Cox p-value = ", pVal),
title.col="black")
}
#or use briss
#delete margins on plots
plot.NoMargins <- function(){
par(mgp=c(1.5,0.5,0), #axis title, axis labels and axis line
mai=c(0.6,0.6,0.02,0.02) #c(bottom, left, top, right)
)
}
# Example:
# pdf( file=paste(resultsFolder, "SNF_all_kaplan.pdf", sep="") , onefile=TRUE)
# plot.SurvivalK3 (groups, mainTitle, survivalData, "bottomright")
# dev.off()
#TRY THIS: http://www.r-statistics.com/2013/07/creating-good-looking-survival-curves-the-ggsurv-function/
bigplot.SurvivalK3 <- function(groups, mainTitle, survData,llocation){
require("survival")
par(mgp=c(2.0,0.4,0), #axis title, axis labels and axis line
mai=c(1.0,1.0,0.02,0.02) #c(bottom, left, top, right)
)
colorsP <- c("red","black","blue")
labelClu <- c("group1","group2", "group3")
coxFit <- coxph(Surv(time = Survival, event = Death) ~ groups,
data = survData, ties="exact")
resCox <- mySumm(coxFit)
mfit <- survfit(Surv(time = Survival, event = Death) ~ groups, data = survData)
pVal<-format.pval(resCox$sctest[3], digits = 4)
plot(mfit, col=colorsP, lty = 1, #line type solid
lwd=15, #line tightness
main = mainTitle,
cex.lab=2.5, #font size
xaxt="n", #do not plot x axe
xlab = "Years",conf.int=F,
ylab="Survival probability"
)
axis(1, at = seq(0, max(survData$Survival), length.out = 15),
labels=round(seq(0, 4.1, length.out = 15), digits=1), las=1) #customize x axe
legend(llocation,
labelClu, col=colorsP,
text.col=colorsP,
text.font=2, #bold
fill=colorsP, #boxes
cex=2,
title = paste(" Cox p-value = ", pVal),
title.col="black")
}
############
#data2plot for example gene expression data
#colorscales Greens, Greys https://plot.ly/r/heatmaps/
#save with preview
potheatmap <- function(data2plot, colorscale="Hot"){
require("plotly")
f <- list(
family = "Calibri, monospace",
size = 30,
color= "#000000"
)
x <- list(
tickfont = f
)
plot_ly(z = data2plot, cl.lim = c(0,1),colorscale=colorscale,
type = "heatmap")
layout(xaxis = x, yaxis = x)
}
plotLineSimple <- function (line1, line2, line3){
loadls("Rcpp ggplot2")
minY = min(line1, line2, line3)
maxY = max(line1, line2, line3)
plot (df1, xlim=c(minX, maxX), ylim=c(minY, maxY))
plot(line1,type="l", col = "red", ylim=c(minY, maxY))
lines(line2, col = "blue")
lines(line3, col = "green")
}
plotLines <- function (line1, line2, line3,lab1, lab2, lab3,xlab, ylab ){
loadls("Rcpp ggplot2")
df <- data.frame(x = rep(c(1:length(line1)), 3), y = c(line1, line2, line3),
colors = rep(c(lab1, lab2, lab3), each=length(line1)) )
plines <- ggplot(data = df, aes(x=x, y=y, col=colors ) )+ geom_line() + xlab(xlab) + ylab(ylab)
plines + guides(fill=guide_legend(title=NULL))
plines + theme(legend.title=element_blank())
plines
}
# Example
# full documentation http://www.inside-r.org/packages/cran/VennDiagram/docs/venn.diagram
plot.venn <- function (){
require("VennDiagram")
# sample four-set Venn Diagram
A <- sample(1:1000, 400, replace = FALSE);
B <- sample(1:1000, 600, replace = FALSE);
C <- sample(1:1000, 350, replace = FALSE);
D <- sample(1:1000, 550, replace = FALSE);
E <- sample(1:1000, 375, replace = FALSE);
venn.plot <- venn.diagram(
x = list(
A = A,
D = D,
B = B,
C = C
),
filename = "Venn_4set_pretty.tiff",
col = "transparent",
fill = c("cornflowerblue", "green", "yellow", "darkorchid1"),
alpha = 0.50,
# label.col = c("orange", "white", "darkorchid4", "white",
# "white", "white", "white", "white", "darkblue", "white",
# "white", "white", "white", "darkgreen", "white"),
cex = 0,
#fontfamily = "serif",
#fontface = "bold",
cat.col = c("darkblue", "darkgreen", "orange", "darkorchid4"),
cat.cex = 0,
cat.pos = 0,
cat.dist = 0.07,
#cat.fontfamily = "serif",
# rotation.degree = 270,
margin = 0.2
);
plot(venn.plot)
# sample five-set Venn Diagram
venn.plot <- venn.diagram(
x = list(
A = A,
B = B,
C = C,
D = D,
E = E
),
filename = "Venn_5set_pretty.tiff",
col = "black",
fill = c("dodgerblue", "goldenrod1", "darkorange1", "seagreen3", "orchid3"),
alpha = 0.50,
cex = c(1.5, 1.5, 1.5, 1.5, 1.5, 1, 0.8, 1, 0.8, 1, 0.8, 1, 0.8,
1, 0.8, 1, 0.55, 1, 0.55, 1, 0.55, 1, 0.55, 1, 0.55, 1, 1, 1, 1, 1, 1.5),
cat.col = c("dodgerblue", "goldenrod1", "darkorange1", "seagreen3", "orchid3"),
cat.cex = 1.5,
cat.fontface = "bold",
margin = 0.05
);
plot(venn.plot)
}
#save plot pdf nice
#pdf( file=paste(resultsFolder, "plots2.pdf", sep="") , onefile=TRUE)
# layout(matrix(c(1,1,2,3),ncol=2,byrow=TRUE),widths=c(1,1),heights=c(3,2))
# par(mar=c(0,0,5,0))
# plot.new(); plot.window(xlim=c(0,1),ylim=c(0,1))
# title("Title",cex.main=1.5)
# text(0.4,0.5,adj=c(0,0),lab=
# "> print(myList)
# $a
# [1] 1
#
# $b
# [1] 2")
# par(mar=c(5,4,1,1))
# boxplot(1:10)
# hist(1:10)
# dev.off()
Raw
survival_analysis.R
geneSetFormat <- function(genesOfInterest){
genesOfInterest <- gsub("hsa:","",genesOfInterest)
genesOfInterest <- lapply(genesOfInterest, find_sym )
genesOfInterest <- unlist(genesOfInterest[!is.na(genesOfInterest)])
#Number of mRNA on pathway
genesOfInterest <- intersect(colnames(geneEData), as.character(genesOfInterest))
genesOfInterest
}
#get symbols
#x <- org.Hs.egSYMBOL
# Get the gene symbol that are mapped to an entrez gene identifiers
# mapped_genes <- mappedkeys(x)
#example
find_sym <- function(g){
require(org.Hs.eg.db)
SYMB <- org.Hs.egSYMBOL[[g]]
if(is.null(SYMB) ){
return(NA)
}else{
return(SYMB)
}
}
#function copied from survival:::summary.coxph
mySumm <- function (object, conf.int = 0.95, scale = 1, ...)
{
cox <- object
beta <- cox$coefficients
if (is.null(cox$coefficients)) {
return(object)
}
nabeta <- !(is.na(beta))
beta2 <- beta[nabeta]
if (is.null(beta) | is.null(cox$var))
stop("Input is not valid")
se <- sqrt(diag(cox$var))
if (!is.null(cox$naive.var))
nse <- sqrt(diag(cox$naive.var))
rval <- list(call = cox$call, fail = cox$fail, na.action = cox$na.action,
n = cox$n, loglik = cox$loglik)
if (!is.null(cox$nevent))
rval$nevent <- cox$nevent
if (is.null(cox$naive.var)) {
tmp <- cbind(beta, exp(beta), se, beta/se, 1 - pchisq((beta/se)^2,
1))
dimnames(tmp) <- list(names(beta), c("coef", "exp(coef)",
"se(coef)", "z", "Pr(>|z|)"))
}
else {
tmp <- cbind(beta, exp(beta), nse, se, beta/se, 1 - pchisq((beta/se)^2,
1))
dimnames(tmp) <- list(names(beta), c("coef", "exp(coef)",
"se(coef)", "robust se", "z", "Pr(>|z|)"))
}
rval$coefficients <- tmp
if (conf.int) {
z <- qnorm((1 + conf.int)/2, 0, 1)
beta <- beta * scale
se <- se * scale
tmp <- cbind(exp(beta), exp(-beta), exp(beta - z * se),
exp(beta + z * se))
dimnames(tmp) <- list(names(beta), c("exp(coef)", "exp(-coef)",
paste("lower .", round(100 * conf.int, 2), sep = ""),
paste("upper .", round(100 * conf.int, 2), sep = "")))
rval$conf.int <- tmp
}
df <- length(beta2)
logtest <- -2 * (cox$loglik[1] - cox$loglik[2])
rval$logtest <- c(test = logtest, df = df, pvalue = 1 - pchisq(logtest,
df))
rval$sctest <- c(test = cox$score, df = df, pvalue = 1 -
pchisq(cox$score, df))
rval$rsq <- c(rsq = 1 - exp(-logtest/cox$n), maxrsq = 1 -
exp(2 * cox$loglik[1]/cox$n))
rval$waldtest <- c(test = as.vector(round(cox$wald.test,
2)), df = df, pvalue = 1 - pchisq(as.vector(cox$wald.test),
df))
if (!is.null(cox$rscore))
rval$robscore <- c(test = cox$rscore, df = df, pvalue = 1 -
pchisq(cox$rscore, df))
rval$used.robust <- !is.null(cox$naive.var)
if (!is.null(cox$concordance)) {
if (is.matrix(cox$concordance))
temp <- colSums(cox$concordance)
else temp <- cox$concordance
rval$concordance <- c(concordance = (temp[1] + temp[3]/2)/sum(temp[1:3]),
se = temp[5]/(2 * sum(temp[1:3])))
}
rval
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment