)
rm(list=ls(all=TRUE))#
graphics.off()#
library("methylumi")#
require("amap")#
require("Hmisc")#
require("MASS")#
require("limma")#
require("lattice")
install.packages("/Users/fmancuso/Downloads/Hmisc_3.7-0.tar.gz",#
                 lib="/Library/Frameworks/R.framework/Versions/2.11/Resources/library/",#
                 repos=NULL,#
                 dependencies=TRUE)
rm(list=ls(all=TRUE))#
graphics.off()#
library("methylumi")#
require("amap")#
require("Hmisc")#
require("MASS")#
require("limma")#
require("lattice")
install.packages("/Users/fmancuso/Downloads/Hmisc_3.7-0.tar.gz",#
                 lib="/Library/Frameworks/R.framework/Versions/2.11/Resources/library/",#
                 repos=NULL,#
                 dependencies=TRUE)
install.packages()
help("install.packages")
install.packages("/Users/fmancuso/Downloads/Hmisc_3.7-0.tar.gz",#
                 lib="/Library/Frameworks/R.framework/Versions/2.11/Resources/library/",#
                 repos="http://cran.r-project.org",#
                 dependencies=TRUE)
install.packages("Hmisc",#
                 lib="/Library/Frameworks/R.framework/Versions/2.11/Resources/library/",#
                 repos="http://cran.r-project.org",#
                 dependencies=TRUE)
install.packages("Hmisc",#
                 lib="/Library/Frameworks/R.framework/Versions/2.11/Resources/library/",#
                 repos="http://cran.r-project.org",#
                 dependencies=F)
install.packages("/Users/fmancuso/Downloads/Hmisc_3.7-0(2).tar.gz",#
                 lib="/Library/Frameworks/R.framework/Versions/2.11/Resources/library/",#
                 repos="NULL",#
                 dependencies=TRUE)
install.packages("/Users/fmancuso/Downloads/Hmisc_3.6-1.tar.gz",#
                 lib="/Library/Frameworks/R.framework/Versions/2.11/Resources/library/",#
                 repos="NULL",#
                 dependencies=TRUE)
install.packages("/Users/fmancuso/Downloads/Hmisc_3.6-1.tar.gz",#
                 lib="/Library/Frameworks/R.framework/Versions/2.11/Resources/library/",#
                 repos=NULL,#
                 dependencies=TRUE)
require(Hmisc)
install.packages("Hmisc", dep=TRUE, type="source")
sweave("/Users/fmancuso/Analysis/Rscripts/datasheet.Rnw")
Sweave("/Users/fmancuso/Analysis/Magda/Prova3.Rnw")
?RweaveLatex
Sweave("/Users/fmancuso/Analysis/Magda/Prova3.Rnw", driver = RweaveLatex(stylepath = TRUE))
Sweave("/Users/fmancuso/Analysis/Magda/Prova3.Rnw", driver = RweaveLatex())
testfile <- system.file("Sweave", "Sweave-test-1.Rnw", package = "utils")#
#
## enforce par(ask=FALSE)#
options(device.ask.default=FALSE)
tools::texi2dvi("Prova.tex", pdf=TRUE)
tools::texi2dvi("Prova3.tex", pdf=TRUE)
Sys.time
Sys.time()
as.POSIXlt(Sys.time(),"GMT")
source("http://bioconductor.org/biocLite.R")#
    biocLite("rtracklayer")
load("/Users/fmancuso/Downloads/R656/R656/R656tr.RData")
R656tr.RData
ls()
R656tr
Bioconductor()
###An example of how to use the function#
list1=letters[1:6];#
list2=letters[6:10];#
list3=letters[6:15];#
list4=letters[6:21];#
list5=letters[6:26];#
plot_venn5(list1, list2, list3, list4, list5);#
####
=================================================#
#
### here is code for the function#
#
plot_venn5=function(listA, listB, listC, listD, listE){#
#
all_ids=unique(c(listA, listB, listC, listD, listE));#
all_n=length(all_ids);#
#
#create counting matrix#
#
all_matrix=cbind(affy=rep(0,all_n),est=rep(0,all_n),glomchip=rep(0,all_n),sage=rep(0,all_n),stanford=rep(0,all_n));#
rownames(all_matrix)=all_ids;#
colnames(all_matrix)=c(substitute(listA), substitute(listB),#
substitute(listC), substitute(listD), substitute(listE)); #
all_matrix[all_ids %in% listA, 1]=1;#
all_matrix[all_ids %in% listB, 2]=1;#
all_matrix[all_ids %in% listC, 3]=1;#
all_matrix[all_ids %in% listD, 4]=1;#
all_matrix[all_ids %in% listE, 5]=1;#
#
#elements only in one data set.#
#
only1= apply(all_matrix,1,sum) == 1;#
nA=sum(all_matrix[only1, 1]);#
nB=sum(all_matrix[only1, 2]);#
nC=sum(all_matrix[only1, 3]);#
nD=sum(all_matrix[only1, 4]);#
nE=sum(all_matrix[only1, 5]);#
#
#elements overlapped in two data set.#
#
only2= apply(all_matrix,1,sum) == 2;#
nAB=sum(only2 & all_matrix[,1] & all_matrix[,2]);#
nAC=sum(only2 & all_matrix[,1] & all_matrix[,3]);#
nAD=sum(only2 & all_matrix[,1] & all_matrix[,4]);#
nAE=sum(only2 & all_matrix[,1] & all_matrix[,5]);#
nBC=sum(only2 & all_matrix[,2] & all_matrix[,3]);#
nBD=sum(only2 & all_matrix[,2] & all_matrix[,4]);#
nBE=sum(only2 & all_matrix[,2] & all_matrix[,5]);#
nCD=sum(only2 & all_matrix[,3] & all_matrix[,4]);#
nCE=sum(only2 & all_matrix[,3] & all_matrix[,5]);#
nDE=sum(only2 & all_matrix[,4] & all_matrix[,5]);#
#
#elements overlapped in three data set.#
#
only3= apply(all_matrix,1,sum) == 3;#
nABC=sum(only3 & all_matrix[,1] & all_matrix[,2] & all_matrix[,3]);#
nABD=sum(only3 & all_matrix[,1] & all_matrix[,2] & all_matrix[,4]);#
nABE=sum(only3 & all_matrix[,1] & all_matrix[,2] & all_matrix[,5]);#
nACD=sum(only3 & all_matrix[,1] & all_matrix[,3] & all_matrix[,4]);#
nACE=sum(only3 & all_matrix[,1] & all_matrix[,3] & all_matrix[,5]);#
nADE=sum(only3 & all_matrix[,1] & all_matrix[,4] & all_matrix[,5]);#
nBCD=sum(only3 & all_matrix[,2] & all_matrix[,3] & all_matrix[,4]);#
nBCE=sum(only3 & all_matrix[,2] & all_matrix[,3] & all_matrix[,5]);#
nBDE=sum(only3 & all_matrix[,2] & all_matrix[,4] & all_matrix[,5]);#
nCDE=sum(only3 & all_matrix[,3] & all_matrix[,4] & all_matrix[,5]);#
#
#elements overlapped in four data set.#
#
only4= apply(all_matrix,1,sum) == 4;#
nABCD=sum(only4 & all_matrix[,1] & all_matrix[,2] & all_matrix[,3] &#
all_matrix[,4]);#
nABCE=sum(only4 & all_matrix[,1] & all_matrix[,2] & all_matrix[,3] &#
all_matrix[,5]);#
nABDE=sum(only4 & all_matrix[,1] & all_matrix[,2] & all_matrix[,4] &#
all_matrix[,5]);#
nACDE=sum(only4 & all_matrix[,1] & all_matrix[,3] & all_matrix[,4] &#
all_matrix[,5]);#
nBCDE=sum(only4 & all_matrix[,2] & all_matrix[,3] & all_matrix[,4] &#
all_matrix[,5]);#
#
#elements overlapped in five data set.#
#
all5= apply(all_matrix,1,sum) == 5;#
nABCDE=sum(all5);#
#
#make the plot.#
#
elps=cbind(150*cos(seq(0,2*pi,len=1000)), 60*sin(seq(0,2*pi,len=1000)));#
#
relocate_elp=function(e, alpha, x, y){#
phi=(alpha/180)*pi;#
xr=e[,1]*cos(phi)+e[,2]*sin(phi);#
yr=-e[,1]*sin(phi)+e[,2]*cos(phi);#
xr=x+xr;#
yr=y+yr;#
return(cbind(xr, yr));#
}#
#
par(mar=c(1,1,1,1)); #
plot(c(0, 400), c(0, 400), type="n", axes=F, ylab="", xlab="");#
#
polygon(relocate_elp(elps, 90,200, 250));#
polygon(relocate_elp(elps, 162,250, 220));#
polygon(relocate_elp(elps, 234,250, 150));#
polygon(relocate_elp(elps, 306,180, 125));#
polygon(relocate_elp(elps, 378,145, 200));#
#
#label the data set name.#
#
text(50, 280, colnames(all_matrix)[1]);#
text(170,400, colnames(all_matrix)[2]);#
text(350,300, colnames(all_matrix)[3]);#
text(350,20, colnames(all_matrix)[4]);#
text(55,10, colnames(all_matrix)[5]);#
#
#label the numbers#
#
text(61, 228, nA);#
text(194, 329, nB);#
text(321, 245, nC);#
text(290, 81, nD);#
text(132, 69, nE);#
#
text(146, 250, nAB, cex=0.7); #
text(123, 188, nAC, cex=0.7); #
text(275, 152, nAD, cex=0.7); #
text(137, 146, nAE, cex=0.7); #
text(243, 268, nBC, cex=0.7); #
text(175, 267, nBD, cex=0.7); #
text(187, 117, nBE, cex=0.7); #
text(286, 188, nCD, cex=0.7); #
text(267, 235, nCE, cex=0.7); #
text(228, 105, nDE, cex=0.7); #
#
text(148, 210, nABC,cex=0.7);#
text(159, 253, nABD,cex=0.7); #
text(171, 141, nABE,cex=0.7); #
text(281, 175, nACD,cex=0.7); #
text(143, 163, nACE,cex=0.7); #
text(252, 145, nADE,cex=0.7); #
text(205, 255, nBCD,cex=0.7); #
text(254, 243, nBCE,cex=0.7); #
text(211, 118, nBDE,cex=0.7); #
text(267, 211, nCDE,cex=0.7); #
#
text(170, 231,nABCD,cex=0.7); #
text(158, 169,nABCE,cex=0.7); #
text(212, 139,nABDE,cex=0.7);#
text(263, 180,nACDE,cex=0.7); #
text(239, 232,nBCDE,cex=0.7);#
#
text(204,190,nABCDE); #
}
plot_venn5=function(listA, listB, listC, listD, listE){#
#
all_ids=unique(c(listA, listB, listC, listD, listE));#
all_n=length(all_ids);#
#
#create counting matrix#
#
all_matrix=cbind(affy=rep(0,all_n),est=rep(0,all_n),glomchip=rep(0,all_n),sage=rep(0,all_n),stanford=rep(0,all_n));#
rownames(all_matrix)=all_ids;#
colnames(all_matrix)=c(substitute(listA), substitute(listB),#
substitute(listC), substitute(listD), substitute(listE)); #
all_matrix[all_ids %in% listA, 1]=1;#
all_matrix[all_ids %in% listB, 2]=1;#
all_matrix[all_ids %in% listC, 3]=1;#
all_matrix[all_ids %in% listD, 4]=1;#
all_matrix[all_ids %in% listE, 5]=1;#
#
#elements only in one data set.#
#
only1= apply(all_matrix,1,sum) == 1;#
nA=sum(all_matrix[only1, 1]);#
nB=sum(all_matrix[only1, 2]);#
nC=sum(all_matrix[only1, 3]);#
nD=sum(all_matrix[only1, 4]);#
nE=sum(all_matrix[only1, 5]);#
#
#elements overlapped in two data set.#
#
only2= apply(all_matrix,1,sum) == 2;#
nAB=sum(only2 & all_matrix[,1] & all_matrix[,2]);#
nAC=sum(only2 & all_matrix[,1] & all_matrix[,3]);#
nAD=sum(only2 & all_matrix[,1] & all_matrix[,4]);#
nAE=sum(only2 & all_matrix[,1] & all_matrix[,5]);#
nBC=sum(only2 & all_matrix[,2] & all_matrix[,3]);#
nBD=sum(only2 & all_matrix[,2] & all_matrix[,4]);#
nBE=sum(only2 & all_matrix[,2] & all_matrix[,5]);#
nCD=sum(only2 & all_matrix[,3] & all_matrix[,4]);#
nCE=sum(only2 & all_matrix[,3] & all_matrix[,5]);#
nDE=sum(only2 & all_matrix[,4] & all_matrix[,5]);#
#
#elements overlapped in three data set.#
#
only3= apply(all_matrix,1,sum) == 3;#
nABC=sum(only3 & all_matrix[,1] & all_matrix[,2] & all_matrix[,3]);#
nABD=sum(only3 & all_matrix[,1] & all_matrix[,2] & all_matrix[,4]);#
nABE=sum(only3 & all_matrix[,1] & all_matrix[,2] & all_matrix[,5]);#
nACD=sum(only3 & all_matrix[,1] & all_matrix[,3] & all_matrix[,4]);#
nACE=sum(only3 & all_matrix[,1] & all_matrix[,3] & all_matrix[,5]);#
nADE=sum(only3 & all_matrix[,1] & all_matrix[,4] & all_matrix[,5]);#
nBCD=sum(only3 & all_matrix[,2] & all_matrix[,3] & all_matrix[,4]);#
nBCE=sum(only3 & all_matrix[,2] & all_matrix[,3] & all_matrix[,5]);#
nBDE=sum(only3 & all_matrix[,2] & all_matrix[,4] & all_matrix[,5]);#
nCDE=sum(only3 & all_matrix[,3] & all_matrix[,4] & all_matrix[,5]);#
#
#elements overlapped in four data set.#
#
only4= apply(all_matrix,1,sum) == 4;#
nABCD=sum(only4 & all_matrix[,1] & all_matrix[,2] & all_matrix[,3] &#
all_matrix[,4]);#
nABCE=sum(only4 & all_matrix[,1] & all_matrix[,2] & all_matrix[,3] &#
all_matrix[,5]);#
nABDE=sum(only4 & all_matrix[,1] & all_matrix[,2] & all_matrix[,4] &#
all_matrix[,5]);#
nACDE=sum(only4 & all_matrix[,1] & all_matrix[,3] & all_matrix[,4] &#
all_matrix[,5]);#
nBCDE=sum(only4 & all_matrix[,2] & all_matrix[,3] & all_matrix[,4] &#
all_matrix[,5]);#
#
#elements overlapped in five data set.#
#
all5= apply(all_matrix,1,sum) == 5;#
nABCDE=sum(all5);#
#
#make the plot.#
#
elps=cbind(150*cos(seq(0,2*pi,len=1000)), 60*sin(seq(0,2*pi,len=1000)));#
#
relocate_elp=function(e, alpha, x, y){#
phi=(alpha/180)*pi;#
xr=e[,1]*cos(phi)+e[,2]*sin(phi);#
yr=-e[,1]*sin(phi)+e[,2]*cos(phi);#
xr=x+xr;#
yr=y+yr;#
return(cbind(xr, yr));#
}#
#
par(mar=c(1,1,1,1)); #
plot(c(0, 400), c(0, 400), type="n", axes=F, ylab="", xlab="");#
#
polygon(relocate_elp(elps, 90,200, 250));#
polygon(relocate_elp(elps, 162,250, 220));#
polygon(relocate_elp(elps, 234,250, 150));#
polygon(relocate_elp(elps, 306,180, 125));#
polygon(relocate_elp(elps, 378,145, 200));#
#
#label the data set name.#
#
text(50, 280, colnames(all_matrix)[1]);#
text(170,400, colnames(all_matrix)[2]);#
text(350,300, colnames(all_matrix)[3]);#
text(350,20, colnames(all_matrix)[4]);#
text(55,10, colnames(all_matrix)[5]);#
#
#label the numbers#
#
text(61, 228, nA);#
text(194, 329, nB);#
text(321, 245, nC);#
text(290, 81, nD);#
text(132, 69, nE);#
#
text(146, 250, nAB, cex=0.7); #
text(123, 188, nAC, cex=0.7); #
text(275, 152, nAD, cex=0.7); #
text(137, 146, nAE, cex=0.7); #
text(243, 268, nBC, cex=0.7); #
text(175, 267, nBD, cex=0.7); #
text(187, 117, nBE, cex=0.7); #
text(286, 188, nCD, cex=0.7); #
text(267, 235, nCE, cex=0.7); #
text(228, 105, nDE, cex=0.7); #
#
text(148, 210, nABC,cex=0.7);#
text(159, 253, nABD,cex=0.7); #
text(171, 141, nABE,cex=0.7); #
text(281, 175, nACD,cex=0.7); #
text(143, 163, nACE,cex=0.7); #
text(252, 145, nADE,cex=0.7); #
text(205, 255, nBCD,cex=0.7); #
text(254, 243, nBCE,cex=0.7); #
text(211, 118, nBDE,cex=0.7); #
text(267, 211, nCDE,cex=0.7); #
#
text(170, 231,nABCD,cex=0.7); #
text(158, 169,nABCE,cex=0.7); #
text(212, 139,nABDE,cex=0.7);#
text(263, 180,nACDE,cex=0.7); #
text(239, 232,nBCDE,cex=0.7);#
#
text(204,190,nABCDE); #
}
list1=letters[1:6];#
list2=letters[6:10];#
list3=letters[6:15];#
list4=letters[6:21];#
list5=letters[6:26];#
plot_venn5(list1, list2, list3, list4, list5);
install.packages("RGG", repos="http://R-Forge.R-project.org")
install.packages("RGG", repos="http://R-Forge.R-project.org", type="source")
install.packages(c("rJava","JavaGD","iplots"))#
install.packages("JGR",type="mac.binary")
install.packages("/Users/fmancuso/Downloads/RGG_0.1.tar.gz",repos = NULL,lib="/Users/fmancuso/Library/R/2.11/library")
library(RGG)
library(JGR)
install.packages('JGR',,'http://www.rforge.net/')
install.packages('JGR',type="source2",'http://www.rforge.net/')
install.packages('JGR',type="source",'http://www.rforge.net/')
install.packages('JGR',,'http://www.rforge.net/',type="source")
install.packages('iplot',,'http://www.rforge.net/',type="source")
remove.packages {JGR}
remove.packages(JGR)
remove.packages(JGR,lib="/Users/fmancuso/Library/R/2.11/library")
require(JGR)
remove.packages(JGR,lib="/Users/fmancuso/Library/R/2.11/library")
install.packages("rattle", dependencies=c("Depends", "Suggests"))
library(rattle)
?require
require(amap,quietly=T)
require(Hmisc,quietly=T)
require(Hmisc,quietly=F)
library(Hmisc,quietly=T, verbose=F)
?require
suppressPackageStartupMessages("amap")
require(amap)
message("ABC", "DEF")#
suppressMessages(message("ABC"))#
#
testit <- function() {#
  message("testing package startup messages")#
  packageStartupMessage("initializing ...", appendLF = FALSE)#
  Sys.sleep(1)#
  packageStartupMessage(" done")#
}#
#
testit()#
suppressPackageStartupMessages(testit())#
suppressMessages(testit()
)
a <- 1:10
a
b <- 10:20
b
b <- 11:20
a*b
as.matrix(a*b)
a %*% b
as.matrix(a) %*% as.matrix(b)
as.matrix(a) %*% b
library(Rscript)
exit
q()
/Users/fmancuso/Analysis/Magda/PaperPipeline/HumMeth27.rda
"/Users/fmancuso/Analysis/Magda/PaperPipeline/HumMeth27.rda"
library(HumMeth27QCReport)
remove.packages(HumMeth27QCReport)
remove.packages("HumMeth27QCReport")
library(HumMeth27QCReport)
require("Biobase")
library(methylumi)
source("http://www.bioconductor.org/biocLite.R")#
biocLite("Biobase")
library(methylumi)
library(Biobase)
library(HumMeth27QCReport)
?hist
log2(1)
require(IlluminaHumanMethylation27k.db)
x <- IlluminaHumanMethylation27kACCNUM
x
str(x)
mapped_probes <- mappedkeys(x)
mapped_probes
??MethilumiM
??MethylumiM
a<-c(2.4,8)
stdv(a)
sd(a)
mean(a)
median(a)
a<-c(2e5,8e9)
median(a)
sd(a)
b<-log10(a)
sd(b)
?seq
require(HumMeth27QCReport)
library(snapCGH)
datadir <- "/Users/fmancuso/Analysis/Imma/AgilentCGH"
targets <- readTargets("targets.txt", path = datadir)
targets
RG1 <- read.maimages(targets$FileName, path = datadir, source = "genepix")
?read.image
?read.maimages
RG1 <- read.maimages(targets$FileName, path = datadir, source = "agilent")
RG1 <- read.maimages(targets$FileName, path = datadir, source = "agilent.median")
warnings()
RG1 <- read.maimages(targets$FileName, path = datadir, source = "agilent")
warnings()
RG1 <- read.maimages(targets$FileName, path = datadir, source = "agilent.median")
RG1 <- read.maimages(targets$FileName, path = datadir, source = "generic")
RG1 <- read.maimages(targets$FileName, path = datadir, source = "agilent")
RG1 <- read.maimages("US45102909_251469373282stp_S01_CGH_107_Sep09.txt", path = datadir, source = "agilent")
RG1 <- read.maimages("aUS45102909_251469373282stp_S01_CGH_107_Sep09.txt", path = datadir, source = "agilent")
RG1 <-read.maimages(targets$File_names,source="agilent",#
columns=list(R="rProcessedSignal",G="gProcessedSignal"))
RG1 <-read.maimages(targets$Filename,source="agilent",#
columns=list(R="rProcessedSignal",G="gProcessedSignal"))
RG1 <-read.maimages(targets$Filename,source="agilent",path = datadir,#
columns=list(R="rProcessedSignal",G="gProcessedSignal"))
RG1 <- read.maimages("US45102909_251469373282stp_S01_CGH_107_Sep09.txt", path = datadir, source = "agilent", columns=list(R="rProcessedSignal",G="gProcessedSignal"))
RG1 <- read.maimages(targets$FileName, path = datadir, source = "agilent", columns=list(R="rProcessedSignal",G="gProcessedSignal"))
RG1 <- read.maimages("US45102909_251469373282stp_S01_CGH_107_Sep09.txt", path = datadir, source = "agilent", columns=list(R="rProcessedSignal",G="gProcessedSignal"))
RG1 <- read.maimages(targets$FileName[1:5], path = datadir, source = "agilent", columns=list(R="rProcessedSignal",G="gProcessedSignal"))
targets$FileName[1:5]
targets <- readTargets("targets.txt", path = datadir)
targets$FileName[1:5]
RG1 <- read.maimages(targets$FileName, path = datadir, source = "agilent", columns=list(R="rProcessedSignal",G="gProcessedSignal"))
RG2 <- readPositionalInfo (RG1,source="agilent")
str(RG2)
RG2$design <- c(-1, -1)
?backgroundCorrect
RG3 <- backgroundCorrect(RG2, method="normexp", normexp.method="saddle", verbose=TRUE)
?normalizeWithinArrays
MA <- normalizeWithinArrays(RG3, method = "median")
?processCGH
MA2 <- processCGH(MA, method.of.averaging = mean, ID = "ID")
MA2 <- processCGH(MA, method.of.averaging = mean, ID = ID)
MA$genes$ProbeName
MA2 <- processCGH(MA, method.of.averaging = mean, ID = MA$genes$ProbeName)
MA$genes$Position
MA$genes$Chr
rm(list=ls(all=TRUE))#
graphics.off()#
#
library("snapCGH")#
datadir <- "/Users/fmancuso/Analysis/Imma/AgilentCGH"#
#
targets <- readTargets("targets.txt", path = datadir)#
RG1 <- read.maimages(targets$FileName, path = datadir, source = "agilent", columns=list(R="rProcessedSignal",G="gProcessedSignal"))#
#RG2 <- readPositionalInfo (RG1,source="agilent")#
RG1$design <- c(-1, -1)#
#
### background correction#
RG2 <- backgroundCorrect(RG1, method="normexp", normexp.method="saddle", verbose=TRUE)#
### normalization interArray#
MA <- normalizeWithinArrays(RG2, method = "median")#
### duplicates are removed following the averaging#
MA2 <- processCGH(MA, method.of.averaging = mean, ID = MA$genes$ProbeName)
 MA$genes$Position
head(MA$genes$ProbeName)
RG2 <- readPositionalInfo (RG1,source="agilent")
RG2
RG3 <- backgroundCorrect(RG2, method="normexp", normexp.method="saddle", verbose=TRUE)
MA <- normalizeWithinArrays(RG3, method = "median")
MA1 <- normalizeWithinArrays(RG3, method = "median")#
MA1$genes$Position <- MA1$genes$Start#
### duplicates are removed following the averaging#
MA2 <- processCGH(MA1, method.of.averaging = mean, ID = MA$genes$ProbeName)
Ma1
MA1
MA2 <- processCGH(MA1, method.of.averaging = mean)
MA$genes$Position
colnames(MA1$genes)[10] <- 'Position'
colnames(MA1$genes)[10]
colnames(MA1$genes)[8]
colnames(MA1$genes)[7]
MA1 <- normalizeWithinArrays(RG3, method = "median")#
colnames(MA1$genes)[7] <- 'Position'#
### duplicates are removed following the averaging#
MA2 <- processCGH(MA1, method.of.averaging = mean, ID = MA$genes$ProbeName)
colnames(MA1$genes)[3]
colnames(MA1$genes)[2]
colnames(MA1$genes)[4]
colnames(MA1$genes)[4] <- 'ID'
MA1 <- normalizeWithinArrays(RG3, method = "median")#
colnames(MA1$genes)[7] <- 'Position'#
### duplicates are removed following the averaging#
MA2 <- processCGH(MA1, method.of.averaging = mean)
MA2 <- processCGH(MA1, method.of.averaging = mean, ID = "ID")
str(MA1$genes)
colnames(MA1$genes)[4] <- 'ID'
str(MA1$genes)
MA2 <- processCGH(MA1, method.of.averaging = mean, ID = "ID")
require(HumMeth27QCReport)
library(HumMeth27QCReport)
require(HumMeth27QCReport)
Dir <-""/Users/fmancuso/Analysis/Magda/prova""
Dir <-""/Users/fmancuso/Analysis/Magda/prova"
require(HumMeth27QCReport)
library(codetools)
checkUsageEnv(.GlobalEnv)
NormCheck <- function(Dir, pval, ClustMethod, ChrX)   {#
	#
	if (missing(ClustMethod) || is.null(ClustMethod)) {#
		ClustMethod <- "euclidean"#
	}#
	if (missing(ChrX) || is.null(ChrX)) {#
		ChrX <- "FALSE"#
	}#
	#
	dataFiles <- ImportData(Dir)#
	AverageBeta <- dataFiles$AverageBeta#
	Ctrl <- dataFiles$Ctrl#
	samps_mod <- dataFiles$samps_mod#
	try(DiscarderII <- data.frame("ToDiscard"=dataFiles$DiscarderII), silent=F)#
#
	lumiMethy <- lumiMethyR(AverageBeta, lib="IlluminaHumanMethylation27k.db")#
#
	lumiMethy.c.adj <- lumiMethyC(lumiMethy)  #### color balance adjustment#
	#
	avgPval <- colMeans(detection(lumiMethy.c.adj),na.rm = T)#
	toKeep <- (avgPval < pval)#
	try(toKeep[as.character(DiscarderII$ToDiscard)] <- FALSE, silent=T)#
	#
	lumiMethy.norm <- lumiMethyN(lumiMethy.c.adj[,toKeep], method="quantile")   #### quantile normalization based on color balance adjusted data#
#
	if (ChrX == "TRUE" || ChrX == "T")   {#
		x <- IlluminaHumanMethylation27kCHR#
		fData(lumiMethy.norm)$chrom <- as.factor(unlist(mget(featureNames(lumiMethy), x)))#
		lumiMethy.norm.f <- data.frame(exprs(lumiMethy.norm), "Chr"= fData(lumiMethy.norm)$chrom)#
		lumiMethy.norm.f <- lumiMethy.norm.f[lumiMethy.norm.f$Chr != "X",]#
		lumiMethy.norm.f <- lumiMethy.norm.f[!is.na(lumiMethy.norm.f$Chr),]#
		lumiMethy.norm.f <- subset(lumiMethy.norm.f, select = -c(Chr))#
	}   #
	#
	if (ChrX == "FALSE" || ChrX == "F")   {#
		lumiMethy.norm <- lumiMethy.norm#
		lumiMethy.norm.f <- data.frame(exprs(lumiMethy.norm))#
	}#
	#
	pdf(file="ExplorativeAnalysis.pdf",width=15,height=9)#
	par(oma=c(1,1,1,1), mar=c(7.1,3.1,3.1,1.1))#
	pca.samples.plot(lumiMethy.norm.f)#
	cluster.samples.plot(lumiMethy.norm.f,main="Hierarchical Clustering",method=ClustMethod)#
	dev.off()#
	#
	NormBetasVal <- as.data.frame(exprs(lumiMethy.norm[toKeep,]))#
	WriteXLS("NormBetasVal", ExcelFileName = "QC_Analysis.xls", SheetNames = "Normalized Beta Values", #
			 verbose = FALSE, Encoding = c("UTF-8", "latin1"), perl = "perl",BoldHeaderRow=T)#
	#
	return(list(lumiMethy.norm=lumiMethy.norm, toKeep=toKeep))#
	ChrX <- NULL#
	rm(ChrX)#
}
checkUsage(NormCheck)
checkUsage(HumMeth27QCReport)
?subset
aaa <- barplot(c(3731,570.3,874.9,3440,363.8,358.7,3255,546.8,794.4))
aaa<-c(3731,570.3,874.9,3440,363.8,358.7,3255,546.8,794.4)
tool <- c("Decon","Distiller","DTASuperCharge","MaxQuant (v. 1.0.13.13)","MaxQuant (v. 1.1.1.14)","MSn_EXTRACT","PD","Raw2MSn","VEMS")
?barplot
barplot(aaa,names.arg =tool, las=2, ylim=c(0,max(aaa)+500), col=orange)
barplot(aaa,names.arg =tool, las=2, ylim=c(0,max(aaa)+500), col="orange")
par(oma=c(1,1,1,1), mar=c(10.1,4.1,3.1,1.1))
barplot(aaa,names.arg =tool, las=2, ylim=c(0,max(aaa)+500), col="orange")
barplot(aaa,names.arg =tool, las=2, ylim=c(0,max(aaa)+500), col="orange", ylab="Size [Mb]")
barplot(aaa,names.arg =tool, las=2, ylim=c(0,max(aaa)+500), col="orange", ylab="Size [Mb]", main="Size of peak list")
abline(h=seq(0,4000,500),col="grey")
bbb <- c(3359.2,195.8,285.8,423.1,182.9,232,198.1,235.1,337.4)
barplot(bbb,names.arg =tool, las=2, ylim=c(0,max(aaa)+500), col="orange", ylab="Size [Mb]", main="Size of peak list")
quartz()
par(oma=c(1,1,1,1), mar=c(10.1,4.1,3.1,1.1))
barplot(bbb,names.arg =tool, las=2, ylim=c(0,max(aaa)+500), col="orange", ylab="Size [Mb]", main="Size of peak list")
abline(h=seq(0,4000,500),col="grey")
tool <- c("DeconMSn","Distiller","DTASuperCharge","MaxQuant (v. 1.0.13.13)","MaxQuant (v. 1.1.1.14)","MSn_EXTRACT","PD","Raw2MSn","VEMS")
barplot(bbb,names.arg =tool, las=2, ylim=c(0,max(aaa)+500), col="orange", ylab="Size [Mb]", main="Size of peak list")
abline(h=seq(0,4000,500),col="grey")
barplot(aaa,names.arg =tool, las=2, ylim=c(0,max(aaa)+500), col="orange", ylab="Size [Mb]", main="Size of peak list")
abline(h=seq(0,4000,500),col="grey")
quartz()
par(oma=c(1,1,1,1), mar=c(10.1,4.1,3.1,1.1))
match005 <- c(48469,41225,10633,49195,49865,50675,46823,49096,47470)
barplot(match005,names.arg =tool, las=2, ylim=c(0,max(aaa)+500), col="lightblue", ylab="# of matches", main="Number of matches at p=0.05")
barplot(match005,names.arg =tool, las=2, ylim=c(0,max(match005)+500), col="lightblue", ylab="# of matches", main="Number of matches at p=0.05")
barplot(match005,names.arg =tool, las=2, ylim=c(0,max(match005)+5000), col="lightblue", ylab="# of matches", main="Number of matches at p=0.05")
barplot(match005,names.arg =tool, las=2, ylim=c(0,max(match005)+10000), col="lightblue", ylab="# of matches", main="Number of matches at p=0.05")
abline(h=seq(0,60000,5000),col="grey")
barplot(match005,names.arg =tool, las=2, ylim=c(0,max(match005)+10000), col="blue", ylab="# of matches", main="Number of matches at p=0.05")
abline(h=seq(0,60000,5000),col="grey")
library(HumMeth27QCReport)
install.packages(”HumMeth27QCReport", dependencies = TRUE)
install.packages("HumMeth27QCReport", dependencies = TRUE)
require(HumMeth27QCReport)
    source("http://www.bioconductor.org/biocLite.R")#
    biocLite("Biobase")
require(HumMeth27QCReport)
source("http://bioconductor.org/biocLite.R")#
pkgs <- rownames(installed.packages())#
biocLite(pkgs)
require(HumMeth27QCReport)
package.skeleton(HumMeth27QCReport_new, path="/Users/fmancuso/Analysis/Magda")
package.skeleton("HumMeth27QCReport_new", path="/Users/fmancuso/Analysis/Magda")
package.skeleton("HumMeth27QCReport_new", path="/Users/fmancuso/Analysis/Magda/", namespace=T)
?package.skeleton()
pwd
rwd
setwd("/Users/fmancuso/Analysis/Magda/HumMeth27QCReport.Rcheck/HumMeth27QCReport/doc")
Sweave(HumMeth27QCReport.Rnw)
Sweave("HumMeth27QCReport.Rnw")
ImportData <- function(Dir) {#
#
	setwd(Dir);#
#
	WorkFiles <- list.files(path = ".",pattern = ".*.txt")#
	AverageBeta <- WorkFiles[grep("AvgBeta", WorkFiles)]#
	Ctrls <- WorkFiles[grep("Control", WorkFiles)]#
	SamplesName <- WorkFiles[grep("Sample", WorkFiles)]#
	try(Discarder <- WorkFiles[grep("Discard", WorkFiles)], silent=T)#
#
	control <- read.table(Ctrls, header = TRUE, sep = "\t") #
	AvBeta <- read.table(AverageBeta, header = TRUE, sep = "\t", flush=T) #
	samps <- read.table(SamplesName, header = TRUE, sep = "\t")#
	DiscarderII <- NULL#
	try(DiscarderII <- read.table(Discarder, header = F, sep = "\t"), silent=T)#
	try(DiscarderII <- as.character(DiscarderII$V1), silent=T)#
#
	samps$SampleLabel <-  samps$Sample.ID#
	colnames(samps)[1] <- "Index"#
	colnames(samps)[2] <- "SampleID"#
#
	nsample <- length(samps$Index)#
	samps2 <- samps[with(samps, order(Index)), ]#
	#
	pdf(file="Sample.pdf",paper="a4r", fonts="Times")#
	SamplePDF <- data.frame(samps2$Index[1:(length(samps2$Index)/2)], samps2$SampleID[1:(length(samps2$SampleID)/2)],rep(" ",(length(samps2$SampleID)/2)), #
							samps2$Index[((length(samps2$Index)/2)+1):length(samps2$Index)], samps2$SampleID[((length(samps2$SampleID)/2)+1):length(samps2$SampleID)])#
	colnames(SamplePDF)<-c("Index","SampleID","","Index","SampleID")#
	textplot(SamplePDF, halign="center", valign="center", show.rownames = F)#
	title("Sample List")#
	dev.off()#
#
	return(list(ctrl=control,samples=samps, AverageBeta=AverageBeta, Ctrls=Ctrls, nsample=nsample, DiscarderII=DiscarderII))#
}
Dir <- "/Users/fmancuso/Analysis/Magda/prova"
HumMeth27QCReport(Dir,"Hum27",0.05,ChrX=F, "euclidean")
HumMeth27QCReport(Dir,"Hum27",0.05,ChrX=F, "euclidean")
ImportDataR <- ImportData(Dir)
str(ImportDataR)
setwd("/Users/fmancuso/Analysis/Magda/HumMeth27QCReport/inst/doc")
Sweave("HumMeth27QCReport.Rnw")
Sweave("HumMeth27QCReport.Rnw")
Sweave("HumMeth27QCReport.Rnw")
