.Rhistory

load("/home/roger/.cache/.fr-AZdJ3l/resPT.RData")
View(resPT)
load("/home/roger/.cache/.fr-2K7h87/breastTCGA.RData")
breastTCGA
paste0("GO:",seq(1:10000))
GO.output <- paste0("GO:",seq(1:10000))
sapply(strsplit(GO.output),head,1)
sapply(strsplit(GO.output,":"),head,1)
time(sapply(strsplit(GO.output,":"),head,1))
system.time(sapply(strsplit(GO.output,":"),head,1))
system.time(sapply(strsplit(GO.output,":"),'[',1))
system.time(sapply(strsplit(GO.output,":"),'[',1))
system.time(sapply(strsplit(GO.output,":"),'[',1))
system.time(sapply(strsplit(GO.output,":"),'[',1))
system.time(sapply(strsplit(GO.output,":"),head,1))
system.time(sapply(strsplit(GO.output,":"),head,1))
system.time(sapply(strsplit(GO.output,":"),head,1))
system.time(vapply(strsplit(GO.output,":",fixed=TRUE),'[','',1))
system.time(vapply(strsplit(GO.output,":",fixed=TRUE),'[','',1))
system.time(vapply(strsplit(GO.output,":",fixed=TRUE),'[','',1))
system.time(vapply(strsplit(GO.output,":",fixed=TRUE),'[','',1))
system.time(vapply(strsplit(GO.output,":"),'[','',1))
system.time(vapply(strsplit(GO.output,":"),'[','',1))
system.time(vapply(strsplit(GO.output,":"),'[','',1))
system.time(vapply(strsplit(GO.output,":"),'[','',1))
system.time(vapply(strsplit(GO.output,":",fixed=TRUE),'[','',1))
system.time(vapply(strsplit(GO.output,":",fixed=TRUE),'[','',1))
system.time(vapply(strsplit(GO.output,":",fixed=TRUE),'[',1,FUN.VALUE=character(1)))
system.time(vapply(strsplit(GO.output,":",fixed=TRUE),'[',1,FUN.VALUE=character(1)))
system.time(vapply(strsplit(GO.output,":",fixed=TRUE),'[',1,FUN.VALUE=character(1)))
system.time(vapply(strsplit(GO.output,":",fixed=TRUE),'[',1,FUN.VALUE=character(1)))
system.time(sapply(strsplit(GO.output,":",fixed=TRUE),'[',1))
system.time(sapply(strsplit(GO.output,":",fixed=TRUE),'[',1))\
system.time(sapply(strsplit(GO.output,":",fixed=TRUE),'[',1))
sapply(strsplit(GO.output,":",fixed=TRUE),'[',1)
sapply(strsplit(GO.output,":")[,1]
sapply(strsplit(GO.output,":")[1,]
strsplit(GO.output,":")[1,]
GO.id.html <- paste0(sapply(strsplit(GO.output,":",fixed=TRUE),'[',1),'&#8203;',sapply(strsplit(GO.output,":",fixed=TRUE),'[',2))
GO.id.html
sapply(strsplit(GO.output,":",fixed=TRUE),'[',2)
library(ChIPQC)
?plotPCA
attributes <- c("Condition","ID")
replace(attributes,"Condition","DBA_CONDITION")
?replace
paste0("DBA_",toupper(attributes))
?makeHyperlinkString
?class
class(df)
library(openxlsx)
?addWorksheet
?writeDataTable
?addStyle
?createStyle
?vector
data.frame(a=c(1,2,3,4),b=c(4,5,6,7))
df <- data.frame(a=c(1,2,3,4),b=c(4,5,6,7))
names(df$a) <- c("pepit","pap","asd","ert")
df
names(df$a)
names(df$)
names(df)
names(df$a)
?writeData
?readLines
matrix(1:30,c(1,10))
matrix(1:30,c(10,1))
head(matrix(1:30,c(10,1)),20)
#setwd("/data/atac-seq/macs2/"); opts <- list(); opts$samplesheet <- "DBsheet_TimAtacSeq1.csv"
setwd("/data/atac-seq/macs2/"); opts <- list(); opts$samplesheet <- "DBsheet_TimAtacSeq1.csv"
?dba.contrast
?show
#setwd("/data/atac-seq/macs2/"); opts <- list(); opts$samplesheet <- "DBsheet_TimAtacSeq1.csv"
diffobject <- dba(sampleSheet=opts$samplesheet)
diffobject$samples
?message
?dba
list()
comparisons <- list()
comparisons
?append
/add
?add
comparisons <- list(comparisons,c("ASDASD","SASDASD"))
comparisons
list()
?list
comparisons <- NA
comparisons <- list(comparisons,c("ASDASD","SASDASD"))
comparisons
comparisons <- NULL
comparisons <- list(comparisons,c("ASDASD","SASDASD"))
comparisons
list(0)
trim()
comparisons
text
comparisons[1] <- c("ASDASD","SASDASD")
comparisons <- list()
comparisons[1] <- c("ASDASD","SASDASD")
comparisons[1]
comparisons[[1] <- c("ASDASD","SASDASD")
comparisons[[1]] <- c("ASDASD","SASDASD")
input <- "Roger vs Goliat"
comparisons[[1]] <- strsplit(gsub(" ","",input),"vs")[[1]]
comparisons
?message
comparisons
names(comparisons)
a <- vector()
a
a <- character(0)
a
c(a,"ASD")
c(a,"ASD","ASD")
diffobject
setwd("/home/roger/Documents/Genomics/MyScripts")
ini=45886885; end=45920909; chrom <- "22"; filename <- "chr22.vcf.gz"
GENOME.class <- readVCF(filename,numcols=10000,tid=chrom,from=ini,to=end,include.unknown=TRUE)
library(PopGenome)
opts$attributes <- "Condition"
attributes <- strsplit(opts$attributes,",")[[1]]
attributes
GENOME.class <- readVCF(filename,numcols=10000,tid=chrom,from=ini,to=end,include.unknown=TRUE)
GENOME.class
individuals <- get.individuals(GENOME.class)[[1]]
get.individuals(GENOME.class)
GENOME.class <- tryCatch({readVCF(filename,numcols=10000,tid=chrom,from=ini,to=end,include.unknown=TRUE)},error = function(e) {
message(e); write(sprintf("%s:%d-%d -- %s",chrom,ini,end,e),sprintf("error_chr%s_%s.log",chrom,pop),append=TRUE)
return(NULL)})
# Verify that the region object is not null (failure to load), contains variants and has been loaded onto R. If it is, then we assume
# that there are no variants, and therefore S,D and all related metrics are 0 (except for alpha)
if (is.null(region)||is.logical(region)||region@n.biallelic.sites==0) { # If readVCF fails, region=FALSE. If no variants, region=NULL
print("This region does not contain any variants: S, Pi and D set to 0")
nwin=floor((end-ini)/wsize)
windows <- cbind(start=seq(ini,end-wsize,by=wsize),end=seq(ini+wsize,end,by=10000))
newrows <- matrix(rep(c(0,0,0,0,0,0,0,NA,NA,0,0,0,0,NA,NA),nwin),nrow=nwin,byrow=TRUE)
colnames(newrows) <- c("S","Pi","DAF","Divsites","D","K","Unknown","Alpha","Fisher","Psel","Pneu","Dsel","Dneu","msel","mneu")
export <- data.frame(population=pop,chr=rep(paste(c("chr",chrom),collapse=""),nwin),windows,newrows)
write.table(export,file=paste(db,".tab",sep=""),quote=FALSE,sep="\t",row.names=F,append=TRUE,
col.names=!file.exists(paste(db,".tab",sep=""))) # Column names written if file does not exist  quit()
quit()
}
region <- tryCatch({readVCF(filename,numcols=10000,tid=chrom,from=ini,to=end,include.unknown=TRUE)},error = function(e) {
message(e); write(sprintf("%s:%d-%d -- %s",chrom,ini,end,e),sprintf("error_chr%s_%s.log",chrom,pop),append=TRUE)
return(NULL)})
# Verify that the region object is not null (failure to load), contains variants and has been loaded onto R. If it is, then we assume
# that there are no variants, and therefore S,D and all related metrics are 0 (except for alpha)
if (is.null(region)||is.logical(region)||region@n.biallelic.sites==0) { # If readVCF fails, region=FALSE. If no variants, region=NULL
print("This region does not contain any variants: S, Pi and D set to 0")
nwin=floor((end-ini)/wsize)
windows <- cbind(start=seq(ini,end-wsize,by=wsize),end=seq(ini+wsize,end,by=10000))
newrows <- matrix(rep(c(0,0,0,0,0,0,0,NA,NA,0,0,0,0,NA,NA),nwin),nrow=nwin,byrow=TRUE)
colnames(newrows) <- c("S","Pi","DAF","Divsites","D","K","Unknown","Alpha","Fisher","Psel","Pneu","Dsel","Dneu","msel","mneu")
export <- data.frame(population=pop,chr=rep(paste(c("chr",chrom),collapse=""),nwin),windows,newrows)
write.table(export,file=paste(db,".tab",sep=""),quote=FALSE,sep="\t",row.names=F,append=TRUE,
col.names=!file.exists(paste(db,".tab",sep=""))) # Column names written if file does not exist  quit()
quit()
}
ini=50278879; end=50290610; chrom <- "22"; filename <- "chr22.vcf.gz"
message(e); write(sprintf("%s:%d-%d -- %s",chrom,ini,end,e),sprintf("error_chr%s_%s.log",chrom,pop),append=TRUE)
return(NULL)})
# Verify that the region object is not null (failure to load), contains variants and has been loaded onto R. If it is, then we assume
# that there are no variants, and therefore S,D and all related metrics are 0 (except for alpha)
if (is.null(region)||is.logical(region)||region@n.biallelic.sites==0) { # If readVCF fails, region=FALSE. If no variants, region=NULL
print("This region does not contain any variants: S, Pi and D set to 0")
nwin=floor((end-ini)/wsize)
windows <- cbind(start=seq(ini,end-wsize,by=wsize),end=seq(ini+wsize,end,by=10000))
newrows <- matrix(rep(c(0,0,0,0,0,0,0,NA,NA,0,0,0,0,NA,NA),nwin),nrow=nwin,byrow=TRUE)
colnames(newrows) <- c("S","Pi","DAF","Divsites","D","K","Unknown","Alpha","Fisher","Psel","Pneu","Dsel","Dneu","msel","mneu")
export <- data.frame(population=pop,chr=rep(paste(c("chr",chrom),collapse=""),nwin),windows,newrows)
write.table(export,file=paste(db,".tab",sep=""),quote=FALSE,sep="\t",row.names=F,append=TRUE,
col.names=!file.exists(paste(db,".tab",sep=""))) # Column names written if file does not exist  quit()
quit()
}
region <- tryCatch({readVCF(filename,numcols=10000,tid=chrom,from=ini,to=end,include.unknown=TRUE)},error = function(e) {
message(e); write(sprintf("%s:%d-%d -- %s",chrom,ini,end,e),sprintf("error_chr%s_%s.log",chrom,pop),append=TRUE)
return(NULL)})
# Verify that the region object is not null (failure to load), contains variants and has been loaded onto R. If it is, then we assume
# that there are no variants, and therefore S,D and all related metrics are 0 (except for alpha)
if (is.null(region)||is.logical(region)||region@n.biallelic.sites==0) { # If readVCF fails, region=FALSE. If no variants, region=NULL
cat("This region does not contain any variants: DAF set to NA\n\n")
export <- data.frame(chr=chrom,start=ini,end=end,DAF=NA)
write.table(export,file=paste(output.name,".tab",sep=""),quote=FALSE,sep="\t",row.names=F,append=TRUE,
col.names=!file.exists(paste(output.name,".tab",sep="")))  # Column names written if file does not exist
quit()
}
region
individuals <- get.individuals(region)[[1]]
humans <- individuals[which(individuals != c("Chimp","Chimp.2"))]
n <- length(humans)
bial <- get.biallelic.matrix(GENOME.class,1) # Biallelic matrix
bialhuman <- bial[1:n,,drop=F] # Remove outgroup (drop = F to keep 1 dimension)
polym <- apply(bialhuman,2,sum)>0# Sites polymorphic in humans: non REF (0) alleles
bialhuman
bial\
bial
get.biallelic.matrix(GENOME.class,1)
GENOME.class <- tryCatch({readVCF(filename,numcols=10000,tid=chrom,from=ini,to=end,include.unknown=TRUE)},error = function(e) {
message(e); write(sprintf("%s:%d-%d -- %s",chrom,ini,end,e),sprintf("error_chr%s_%s.log",chrom,pop),append=TRUE)
return(NULL)})
# Verify that the region object is not null (failure to load), contains variants and has been loaded onto R. If it is, then we assume
# that there are no variants, and therefore S,D and all related metrics are 0 (except for alpha)
if (is.null(GENOME.class)||is.logical(GENOME.class)||GENOME.class@n.biallelic.sites==0) { # If readVCF fails, region=FALSE. If no variants, region=NULL
cat("This region does not contain any variants: DAF set to NA\n\n")
export <- data.frame(chr=chrom,start=ini,end=end,DAF=NA)
write.table(export,file=paste(output.name,".tab",sep=""),quote=FALSE,sep="\t",row.names=F,append=TRUE,
col.names=!file.exists(paste(output.name,".tab",sep="")))  # Column names written if file does not exist
quit()
}
individuals <- get.individuals(GENOME.class)[[1]]
humans <- individuals[which(individuals != c("Chimp","Chimp.2"))]
n <- length(humans)
bial <- get.biallelic.matrix(GENOME.class,1) # Biallelic matrix
bialhuman <- bial[1:n,,drop=F] # Remove outgroup (drop = F to keep 1 dimension)
bialhuman
colSums(is.na(bialhuman))
misshuman <- colSums(is.na(bialhuman))>0
misshuman
sum(misshuman)
polyal <- GENOME.class@region.data@polyallelic.sites[[1]]
polyal
comparison <- c("MOLM1","MUTZ3")
diffobject$masks$comparison[]1
diffobject$masks$comparison[1]
diffobject$masks$eval(comparison[1])
diffobject$masks$as.symbol(comparison[1])
diffobject$masks[as.symbol(comparison[1])]
diffobject$masks[comparison[1]]
diffobject$masks[[comparison[1]]]
diffobject$masks$MOLM1
GENOME.class@n.sites-sum(misshuman,na.rm=T)-polysites
GENOME.class@n.sites
50290610-50278879
?PopGenome
get.diversity(GENOME.class)
system.time(diversity.stats(GENOME.class,pi=TRUE))
system.time(F_ST.stats(GENOME.class)))
system.time(F_ST.stats(GENOME.class))
system.time(F_ST.stats(GENOME.class,mode="nucleotide"))
GENOME.class <- F_ST.stats(GENOME.class,mode="nucleotide")
GENOME.class@Pi
GENOME.class <- diversity.stats(GENOME.class,pi=TRUE)
GENOME.class@Pi
GENOME.class@Pi/GENOME.class@n.sites
GENOME.class@Pi/GENOME.class@n.sites2
GENOME.class@n.sites2
GENOME.class@n.sites
m <- GENOME.class@n.sites-sum(misshuman,na.rm=T)-polysites
freqs <- apply(bialhuman,2,table)
k <- sum(freqs[1,]*freqs[2,])
Pi <- function(k,m,n) {
comb <- choose(n,2) # Binomial coefficient = combination without repetition
pi <- k/(comb*m)
return(round(pi,7))
}
Pi(k,m,n)
polysites <- sum(polyal) # N polyallelic sites in that region
m <- GENOME.class@n.sites-sum(misshuman,na.rm=T)-polysites
freqs <- apply(bialhuman,2,table)
k <- sum(freqs[1,]*freqs[2,])
Pi <- function(k,m,n) {
comb <- choose(n,2) # Binomial coefficient = combination without repetition
pi <- k/(comb*m)
return(round(pi,7))
}
Pi(k,m,n)
k <- sum(freqs[1,]*freqs[2,])
freqs <- apply(bialhuman,2,table)
k <- sum(freqs[1,]*freqs[2,])
freqs
bialhuman <- bial[1:n,,drop=F] # Remove outgroup (drop = F to keep 1 dimension)
misshuman <- colSums(is.na(bialhuman))>0 # Sites missing in humans (e.g. structural variants)
polym <- apply(bialhuman,2,sum)>0 & !misshuman# Sites polymorphic in humans: non REF (0) alleles
bialhuman <- bialhuman[,polym,drop=F] # Remove sites monomorphic in humans
polyal <- GENOME.class@region.data@polyallelic.sites[[1]] # Positions of all polyalleles
polysites <- sum(polyal) # N polyallelic sites in that region
bialhuman
freqs <- apply(bialhuman,2,table)
k <- sum(freqs[1,]*freqs[2,])
k
Pi(k,m,n)
Pi <- function(k,m,n) {
comb <- choose(n,2) # Binomial coefficient = combination without repetition
pi <- k/(comb*m)
return(pi,7))
}
Pi <- function(k,m,n) {
comb <- choose(n,2) # Binomial coefficient = combination without repetition
pi <- k/(comb*m)
return(pi,7)
}
Pi(k,m,n)
Pi <- function(k,m,n) {
comb <- choose(n,2) # Binomial coefficient = combination without repetition
pi <- k/(comb*m)
return(pi)
}
Pi(k,m,n)
GENOME.class@Pi/GENOME.class@n.sites
m
GENOME.class@n.sites
sum(misshuman,na.rm=T)
polysites <- length(polyal) # N polyallelic sites in that region
m <- GENOME.class@n.sites-sum(misshuman,na.rm=T)-polysites
freqs <- apply(bialhuman,2,table)
k <- sum(freqs[1,]*freqs[2,])
Pi <- function(k,m,n) {
comb <- choose(n,2) # Binomial coefficient = combination without repetition
pi <- k/(comb*m)
return(pi)
}
Pi(k,m,n)
GENOME.class@Pi/GENOME.class@n.sites
Pi(k,m+1,n)
Pi(k,m+2,n)
Pi(k,m+10,n)
Pi(k,m-1,n)
Pi(k,m-20,n)
Pi(k,m-100,n)
Pi(k,m-90,n)
Pi(k,m-95,n)
Pi(k,m-96,n)
Pi(k,m-97,n)
?round
ceiling(1/2)
diffobject
diffobject$masks
diffobject$masks$consensus
diffobject$masks$All
?dba.plotVenn
comparisons
comparisons <- c("Roger vs Goliat","Andre vs Levy","Candy vs Cloud")
strsplit(comparisons,"vs")
unlist(strsplit(comparisons,"vs"))
comp.unique <- unique(unlist(strsplit(gsub(" ","",comparison),"vs")))
comp.unique
comparisons
comp.unique <- unique(unlist(strsplit(gsub(" ","",comparisons),"vs")))
comp.unique
diffobject$masks
diffobject$masks[[comp.unique]]
comp.unique
diffobject$masks[[c("MOLM1","MUTZ1")]]
diffobject$masks[["MOLM1"]]
diffobject$masks$comparisons|diffobject$masks$CD34negCD15neg
diffobject$masks$MOLM1|diffobject$masks$CD34negCD15neg
diffobject$masks$None
diffobject$masks$None
names(diffobject$masks)
names(diffobject$masks$narrow
names(diffobject$masks$narrow)
names(diffobject$masks$``)
names(diffobject$masks$All)
diffobject$masks$All
diffobject$masks[[""]]
dba.plotVenn(diffobject,c(TRUE,TRUE,TRUE)
)
dba.plotVenn(diffobject,c(TRUE,TRUE,TRUE,FALSE,FALSE))
dba.plotVenn(diffobject,c(TRUE,TRUE,TRUE,FALSE,FALSE,TRUE))
diffobject$samples
dba.plotVenn(diffobject,diffobject$masks$CD34posCD15neg|diffobject$masks$CD34negCD15neg|diffobject$masks$Bulk,main="Binding site overlaps - All samples")
dba.plotVenn(diffobject,c(TRUE,TRUE,TRUE,FALSE,FALSE,FALSE))
dba.plotVenn(diffobject,c(TRUE,TRUE,TRUE,FALSE,FALSE,TRUE))
diffobject$samples %in% c("MOLM1","Bulk")
diffobject$samples$SampleID %in% c("MOLM1","Bulk")
diffobject$samples[["Condition"]] %in% c("MOLM1","Bulk")
comp.unique <- c("MOLM1","MUTZ3")
dba.plotVenn(diffobject,diffobject$samples$SampleID %in% comp.unique,main="Binding site overlaps - All samples")
diffobject$masks$MUTZ3
diffobject$masks$``
dba.mask(diffobject)
dba.mask(diffobject,attribute="DBA_CONDITION")
dba.mask(diffobject,attribute=DBA_CONDITION)
dba.mask(diffobject)
dba.mask(diffobject,DBA_CONDITION)
dba.mask(diffobject,DBA_CONDITION,"MUTZ3")
comparison
dba.mask(diffobject,attribute,comparison[1])
attributes
attributes <- paste0("DBA_",toupper(attributes))
dba.mask(diffobject,attribute,comparison[1])
attribute = attributes[1]
attribute
dba.mask(diffobject,attribute,comparison[1]))
dba.mask(diffobject,attribute,comparison[1])
attribute
dba.mask(diffobject,attribute=attribute,comparison[1])
dba.mask(diffobject,attribute=attribute,value=comparison[1])
dba.mask(diffobject,as.symbol(attribute),comparison[1])
dba.mask(diffobject,attribute=as.symbol(attribute),comparison[1])
as.symbol(attribute)
dba.mask(diffobject,attribute=eval(as.symbol(attribute)),comparison[1])
DBA_CONDITION
?dba.contrast
?regex
?regexpr
comparison
input
input2 <- "1 vs 2"
split <- strsplit(gsub(" ","",input2),"vs")[[1]]
split
grep(split,"[0-9]",perl=TRUE)
grep("[0-9]",split,perl=TRUE)
split[0] = "r"
grep("[0-9]",split,perl=TRUE)
split
grep("[0-9]",split,perl=TRUE)
grep("a",split,perl=TRUE)
grep("1",split,perl=TRUE)
grep("2",split,perl=TRUE)
grep([0-9],split,perl=TRUE)
grep("\\d",split,perl=TRUE)
split
split[0]
split[1]
split[1] = "r"
grep("\\d",split,perl=TRUE)
colnames(diffobject$samples)
diffobject$samples$SampleID
split
split <- c("ASD","ASD2")
split
grep("[0-9]",split)
grep("[0-9]",split,fixed=TRUE)
grep("\\w",split,fixed=TRUE)
split <- c("1","2")
grep("[0-9]",split,fixed=TRUE)
grep("[0-9]",split,fixed=TRUE,perl=TRUE)
grep("[0-9]",split,perl=TRUE)
split
split <- c("ASD","ASD2")
grep("[0-9]",split,perl=TRUE)
grep("^[0-9]$",split,perl=TRUE)
?dba.plotPCA
dba.plotPCA(diffobject,attributes = attributes) # CHANGE???
attributes
dba.plotPCA(diffobject,attributes = DBA_CONDITION) # CHANGE???
?dba.plotPCA
attributes <- c("DBA_CONDITION","DBA_ID")
attributes
eval(as.symbol(attributes))
DBA_ID
?eval
sapply(attributes,function(x)eval(as.symbol(x)))
atts <- sapply(attributes,function(x)eval(as.symbol(x)))
atts
?dba.plotPCA
?plotPrincomp
?CHIPQC
?ChIPQCreport
source("http://bioconductor.org/biocLite.R")
biocLite("Basic4Cseq")
library(base64enc)
library(Basic4Cseq)
biocLite("BSgenome.Hsapiens.UCSC.hg19")
library(BSgenome.Hsapiens.UCSC.hg19)
fragmentData = createVirtualFragmentLibrary(chosenGenome = Hsapiens, firstCutter = targets[1],
secondCutter = targets[2], readLength = 2,  chromosomeName = "test", libraryName = "")
library(Basic4Cseq)
fragmentData = createVirtualFragmentLibrary(chosenGenome = Hsapiens, firstCutter = targets[1],
secondCutter = targets[2], readLength = 2,  chromosomeName = "test", libraryName = "")