R/FunctionsForBatchAnalysis.R
In supatt-lab/r3Cseq: Analysis of Chromosome Conformation Capture and Next-generation Sequencing (3C-seq)
# TODO: These following functions were implemented to get the interaction regions for replicates data.
# Author: Supat Thongjuea
# Contact : supat.thongjuea@imm.ox.ac.uk or supat.thongjuea@gmail.com
#####
#####The functions below are completely migrated to BioC 3.9 on R 3.6
#####
########################################
setGeneric(
name="getBatchRawReads",
def=function(object){
standardGeneric("getBatchRawReads")
}
)
setMethod("getBatchRawReads",
signature(object = "r3CseqInBatch"),
function (object){
if(!is(object,"r3CseqInBatch")){
stop("Need to initialize the r3CseqInBatch object")
}
######check files########################
bamDir<-bamFilesDirectory(object)
##########get labels#######################
expFiles<-BamExpFiles(object)
contrFiles<-BamContrFiles(object)
#########checking directory and files####
if (file.exists(bamDir)){
setwd(file.path(bamDir))
}else{
stop(paste("Couldn't find path","-->",bamDir))
}
objName <- deparse(substitute(object))
exp.bam.files <-BamExpFiles(object)
contr.bam.files <-BamContrFiles(object)
for(i.file in c(exp.bam.files)){
if(file.exists(i.file) ==FALSE){
stop(paste("Couldn't find file","-->",i.file))
}
}
for(i.file in c(contr.bam.files)){
if(file.exists(i.file) ==FALSE){
stop(paste("Couldn't find file","-->",i.file))
}
}
exp.lib.size<-c()
exp.read.length<-c()
for(i in 1:length(exp.bam.files)){
input.bam<-exp.bam.files[i]
print(paste("start reading in...",input.bam, "....file in the experiment"))
what <- c("rname", "strand", "pos", "qwidth","seq")
param <- ScanBamParam(what=what,flag = scanBamFlag(isUnmappedQuery = FALSE))
exp.bam <- scanBam(input.bam, param=param)
exp.length <-width(exp.bam[[1]]$seq[1])
exp.GRanges<-GRanges(seqnames=as.vector(exp.bam[[1]]$rname),IRanges(start=exp.bam[[1]]$pos,width=exp.length),strand=exp.bam[[1]]$strand)
fileName<-gsub(".bam","",expFiles[i])
lib.size <-length(exp.GRanges)
exp.lib.size<-append(exp.lib.size,lib.size)
exp.read.length<-append(exp.read.length,exp.length)
save(exp.GRanges,file=paste(fileName,".rData",sep=""))
}
contr.lib.size<-c()
contr.read.length<-c()
###########Make count table for each bam file in the controls########
for(i in 1:length(contr.bam.files)){
input.bam<-contr.bam.files[i]
print(paste("start reading in...",input.bam, "....file in the control"))
what <- c("rname", "strand", "pos", "qwidth","seq")
param <- ScanBamParam(what=what,flag = scanBamFlag(isUnmappedQuery = FALSE))
contr.bam <- scanBam(input.bam, param=param)
contr.length <-width(contr.bam[[1]]$seq[1])
contr.GRanges<-GRanges(seqnames=as.vector(contr.bam[[1]]$rname),IRanges(start=contr.bam[[1]]$pos,width=contr.length),strand=contr.bam[[1]]$strand)
fileName<-gsub(".bam","",contrFiles[i])
lib.size <-length(contr.GRanges)
contr.lib.size<-append(contr.lib.size,lib.size)
contr.read.length<-append(contr.read.length,contr.length)
save(contr.GRanges,file=paste(fileName,".rData",sep=""))
}
expBatchLibrarySize(object)<-exp.lib.size
contrBatchLibrarySize(object)<-contr.lib.size
expBatchReadLength(object)<-exp.read.length
contrBatchReadLength(object)<-contr.read.length
assign(objName,object,envir=parent.frame())
invisible(1)
print("The raw read files are created!!!.")
}
)
setGeneric(
name="getBatchReadCountPerRestrictionFragment",
def=function(object,getReadsMethod = c("wholeReads", "adjacentFragmentEndsReads"),
nFragmentExcludedReadsNearViewpoint=2){
standardGeneric("getBatchReadCountPerRestrictionFragment")
}
)
setMethod("getBatchReadCountPerRestrictionFragment",
signature(object = "r3CseqInBatch"),
function (object,getReadsMethod,nFragmentExcludedReadsNearViewpoint){
if(!is(object,"r3CseqInBatch")){
stop("Need to initialize the r3CseqInBatch object")
}
selected_methods <- match.arg(getReadsMethod)
if(selected_methods==""){
selected_methods<-"wholeReads"
}
######check the number of fragment to exclude reads around the viewpoint###
nExcludedFragments <- nFragmentExcludedReadsNearViewpoint
if(nExcludedFragments <0 | nExcludedFragments>5){
stop("The number of fragment that uses for removing reads around the viewpoint must be 0-5 restriction fragments.")
}
######check files########################
objName <- deparse(substitute(object))
bamDir<-bamFilesDirectory(object)
#########checking directory and files####
if (file.exists(bamDir)){
setwd(file.path(bamDir))
}else{
stop(paste("Couldn't find path","-->",bamDir))
}
exp.bam.files <-BamExpFiles(object)
contr.bam.files <-BamContrFiles(object)
for(i.file in c(exp.bam.files)){
i.file<-gsub(".bam","",i.file)
i.file<-paste(i.file,".rData",sep="")
if(file.exists(i.file) ==FALSE){
stop(paste("Couldn't find file","-->",i.file, ". Please run getBatchRawReads function"))
}
}
for(i.file in c(contr.bam.files)){
i.file<-gsub(".bam","",i.file)
i.file<-paste(i.file,".rData",sep="")
if(file.exists(i.file) ==FALSE){
stop(paste("Couldn't find file","-->",i.file, ". Please run getBatchRawReads function"))
}
}
####Get organism#######################
orgName <- organismName(object)
####Initialized Restriction Enzyme Object##
enzymeDb <-new("repbaseEnzyme")
resEnzyme <-restrictionEnzyme(object)
####Build GRanges for restriction site#####
print("Fragmenting genome..........")
wholeFragments<-getWholeGenomeRestrictionFragments(enzymeDb,resEnzyme,orgName)
##########get labels#######################
expFiles<-BamExpFiles(object)
contrFiles<-BamContrFiles(object)
###########################################
exp.read.length<-expBatchReadLength(object)
contr.read.length<-contrBatchReadLength(object)
if(selected_methods=="wholeReads"){
###########Making GRange Object############
wholeFragments.GRanges<-GRanges(seqnames = as.character(wholeFragments$chromosome),IRanges(start=wholeFragments$start,end=wholeFragments$end))
###########Make count table for each bam file in the experiments########
readCountTable<-data.frame(chromosome=seqnames(wholeFragments.GRanges),start=start(wholeFragments.GRanges),end=end(wholeFragments.GRanges))
for(i in 1:length(exp.bam.files)){
print(paste("start counting reads in...",exp.bam.files[i], "....file in the experiment"))
fileName<-gsub(".bam","",expFiles[i])
fileName<-paste(fileName,".rData",sep="")
load(fileName)
exp.GRanges<-excludeReadsNearViewpoint(object,exp.GRanges,nExcludedFragments)
readsPerFragment <- countOverlaps(wholeFragments.GRanges,subject=exp.GRanges)
reads<-data.frame(nReads=readsPerFragment)
Name<-gsub(".bam","",expFiles[i])
colnames(reads)[1]<-Name
readCountTable<-cbind(readCountTable,reads)
}
###########Make count table for each bam file in the controls########
for(i in 1:length(contr.bam.files)){
print(paste("start counting reads in...",contr.bam.files[i], "....file in the control"))
fileName<-gsub(".bam","",contrFiles[i])
fileName<-paste(fileName,".rData",sep="")
load(fileName)
contr.GRanges<-excludeReadsNearViewpoint(object,contr.GRanges,nExcludedFragments)
readsPerFragment <- countOverlaps(wholeFragments.GRanges,subject=contr.GRanges)
reads<-data.frame(nReads=readsPerFragment)
Name<-gsub(".bam","",contrFiles[i])
colnames(reads)[1]<-Name
readCountTable<-cbind(readCountTable,reads)
}
#########Filter out fragment with no read###########################
n<-ncol(readCountTable)
readCountTable.f<-readCountTable[rowSums(readCountTable[,4:n]) >0,]
readCountTable.GRanges<-GRanges(seqnames=readCountTable.f$chromosome,
IRanges(start=readCountTable.f$start,end=readCountTable.f$end))
values(readCountTable.GRanges)<-readCountTable.f[,4:n]
readCountTable(object)<-readCountTable.GRanges
assign(objName,object,envir=parent.frame())
invisible(1)
print("The countTable is created!!!.")
}
if(selected_methods=="adjacentFragmentEndsReads"){
###########Making GRange Object############
wholeFragments.GRanges<-GRanges(seqnames=as.character(wholeFragments$chromosome),
IRanges(start=wholeFragments$start,end=wholeFragments$end))
###########Make count table for each bam file in the experiments########
readCountTable<-data.frame(chromosome=seqnames(wholeFragments.GRanges),
start=start(wholeFragments.GRanges),end=end(wholeFragments.GRanges))
for(i in 1:length(exp.bam.files)){
print(paste("start counting reads in...",exp.bam.files[i], "....file in the experiment"))
fileName<-gsub(".bam","",expFiles[i])
fileName<-paste(fileName,".rData",sep="")
load(fileName)
exp.GRanges<-excludeReadsNearViewpoint(object,exp.GRanges,nExcludedFragments)
########Make RangedData for the 5' and 3' of RE#######
wholeFragments_5_prime.GRanges<-GRanges(seqnames=seqnames(wholeFragments.GRanges),
IRanges(start=start(wholeFragments.GRanges),
end=start(wholeFragments.GRanges)+exp.read.length[i]))
wholeFragments_3_prime.GRanges<-GRanges(seqnames=seqnames(wholeFragments.GRanges),
IRanges(start=end(wholeFragments.GRanges)-exp.read.length[i],
end=end(wholeFragments.GRanges)))
readsPerFragment_5 <- countOverlaps(wholeFragments_5_prime.GRanges,subject=exp.GRanges)
readsPerFragment_3 <- countOverlaps(wholeFragments_3_prime.GRanges,subject=exp.GRanges)
combined_reads<-readsPerFragment_5+readsPerFragment_3
reads<-data.frame(nReads=combined_reads)
fileName<-gsub(".bam","",expFiles[i])
colnames(reads)[1]<-fileName
readCountTable<-cbind(readCountTable,reads)
}
for(i in 1:length(contr.bam.files)){
print(paste("start counting reads in...",contr.bam.files[i], "....file in the control"))
fileName<-gsub(".bam","",contrFiles[i])
fileName<-paste(fileName,".rData",sep="")
load(fileName)
contr.GRanges<-excludeReadsNearViewpoint(object,contr.GRanges,nExcludedFragments)
########Make RangedData for the 5' and 3' of RE#######
wholeFragments_5_prime.GRanges<-GRanges(seqnames=seqnames(wholeFragments.GRanges),
IRanges(start=start(wholeFragments.GRanges),
end=start(wholeFragments.GRanges)+contr.read.length[i]))
wholeFragments_3_prime.GRanges<-GRanges(seqnames=seqnames(wholeFragments.GRanges),
IRanges(start=end(wholeFragments.GRanges)-contr.read.length[i],
end=end(wholeFragments.GRanges)))
readsPerFragment_5 <- countOverlaps(wholeFragments_5_prime.GRanges,subject=contr.GRanges)
readsPerFragment_3 <- countOverlaps(wholeFragments_3_prime.GRanges,subject=contr.GRanges)
combined_reads<-readsPerFragment_5+readsPerFragment_3
reads<-data.frame(nReads=combined_reads)
fileName<-gsub(".bam","",contrFiles[i])
colnames(reads)[1]<-fileName
readCountTable<-cbind(readCountTable,reads)
}
#########Filter out fragment with no read###########################
n<-ncol(readCountTable)
readCountTable.f<-readCountTable[rowSums(readCountTable[,4:n]) >0,]
readCountTable.GRanges<-GRanges(seqnames=readCountTable.f$chromosome,
IRanges(start=readCountTable.f$start,end=readCountTable.f$end))
values(readCountTable.GRanges)<-readCountTable.f[,4:n]
readCountTable(object)<-readCountTable.GRanges
assign(objName,object,envir=parent.frame())
invisible(1)
print("The countTable is created!!!.")
}
}
)
setGeneric(
name="getBatchReadCountPerWindow",
def=function(object,windowSize=5e3,nFragmentExcludedReadsNearViewpoint=2,mode=c("non-overlapping", "overlapping")){
standardGeneric("getBatchReadCountPerWindow")
}
)
setMethod("getBatchReadCountPerWindow",
signature(object = "r3CseqInBatch"),
function (object,windowSize,nFragmentExcludedReadsNearViewpoint,mode){
objName <- deparse(substitute(object))
if(!is(object,"r3CseqInBatch")){
stop("Need to initialize the r3CseqInBatch object")
}
######check the number of fragment to exclude reads around the viewpoint###
nExcludedFragments <- nFragmentExcludedReadsNearViewpoint
if(nExcludedFragments <0 | nExcludedFragments>5){
stop("The number of fragment that uses for removing reads around the viewpoint must be 0-5 restriction fragments.")
}
#####Select the read count methods#########
mode <- match.arg(mode)
if(mode==""){
mode<-"non-overlapping"
}
######check files########################
objName <- deparse(substitute(object))
bamDir<-bamFilesDirectory(object)
#########checking directory and files####
if (file.exists(bamDir)){
setwd(file.path(bamDir))
}else{
stop(paste("Couldn't find path","-->",bamDir))
}
exp.bam.files <-BamExpFiles(object)
contr.bam.files <-BamContrFiles(object)
for(i.file in c(exp.bam.files)){
i.file<-gsub(".bam","",i.file)
i.file<-paste(i.file,".rData",sep="")
if(file.exists(i.file) ==FALSE){
stop(paste("Couldn't find file","-->",i.file, ". Please run getBatchRawReads function"))
}
}
for(i.file in c(contr.bam.files)){
i.file<-gsub(".bam","",i.file)
i.file<-paste(i.file,".rData",sep="")
if(file.exists(i.file) ==FALSE){
stop(paste("Couldn't find file","-->",i.file, ". Please run getBatchRawReads function"))
}
}
##########get labels#######################
expFiles<-BamExpFiles(object)
contrFiles<-BamContrFiles(object)
###########################################
wholeFragments.GRanges<-getFragmentsPerWindow(object,windowSize,mode)
###########Make count table for each bam file in the experiments########
readCountTable<-data.frame(chromosome=seqnames(wholeFragments.GRanges),
start=start(wholeFragments.GRanges),
end=end(wholeFragments.GRanges))
for(i in 1:length(exp.bam.files)){
print(paste("start counting reads in...",exp.bam.files[i], "....file in the experiment"))
fileName<-gsub(".bam","",expFiles[i])
fileName<-paste(fileName,".rData",sep="")
load(fileName)
exp.GRanges<-excludeReadsNearViewpoint(object,exp.GRanges,nExcludedFragments)
readsPerFragment <- countOverlaps(wholeFragments.GRanges,subject=exp.GRanges)
reads<-data.frame(nReads=readsPerFragment)
Name<-gsub(".bam","",expFiles[i])
colnames(reads)[1]<-Name
readCountTable<-cbind(readCountTable,reads)
}
###########Make count table for each bam file in the controls########
for(i in 1:length(contr.bam.files)){
print(paste("start counting reads in...",contr.bam.files[i], "....file in the control"))
fileName<-gsub(".bam","",contrFiles[i])
fileName<-paste(fileName,".rData",sep="")
load(fileName)
contr.GRanges<-excludeReadsNearViewpoint(object,contr.GRanges,nExcludedFragments)
readsPerFragment <- countOverlaps(wholeFragments.GRanges,subject=contr.GRanges)
reads<-data.frame(nReads=readsPerFragment)
Name<-gsub(".bam","",contrFiles[i])
colnames(reads)[1]<-Name
readCountTable<-cbind(readCountTable,reads)
}
#########Filter out fragment with no read###########################
n<-ncol(readCountTable)
readCountTable.f<-readCountTable[rowSums(readCountTable[,4:n]) >0,]
readCountTable.GRanges<-GRanges(seqnames=readCountTable.f$chromosome,
IRanges(start=readCountTable.f$start,
end=readCountTable.f$end))
values(readCountTable.GRanges)<-readCountTable.f[,4:n]
readCountTable(object)<-readCountTable.GRanges
assign(objName,object,envir=parent.frame())
invisible(1)
print("The countTable is created!!!.")
}
)
setGeneric(
name="calculateBatchRPM",
def=function(object,normalized_method=c("powerlawFittedRPM","normalRPM")){
standardGeneric("calculateBatchRPM")
}
)
setMethod("calculateBatchRPM",
signature(object = "r3CseqInBatch"),
function (object,normalized_method){
if(!is(object,"r3CseqInBatch")){
stop("Need the r3CseqInBatch object")
}
#####Select the method#########
selected_methods <- match.arg(normalized_method)
if(selected_methods==""){
selected_methods<-"powerlawFittedRPM"
}
##########get labels#######################
expFiles<-BamExpFiles(object)
contrFiles<-BamContrFiles(object)
##########Get library size#################
expLibs<-expBatchLibrarySize(object)
contrLibs<-contrBatchLibrarySize(object)
###########################################
objName <- deparse(substitute(object))
count.data<-readCountTable(object)
if(length(count.data)==0){
stop("No reads found in the countTable object, please run getBatchReadCountPerRestrictionFragment or getBatchReadCountPerWindow function!!")
}
if(selected_methods=="powerlawFittedRPM"){
RPMsTable<-data.frame(chromosome=seqnames(count.data),start=start(count.data),end=end(count.data))
for(i in 1:length(expFiles)){
print(paste("Calculating RPMs reads in...",expFiles[i], "....file in the experiment"))
fileName<-gsub(".bam","",expFiles[i])
count.data.f<-mcols(count.data)
count.i<-count.data.f[,colnames(count.data.f)==fileName]
values<-count.i
coef.i<-getPowerLawFittedCoeficient(values)
RPMs<-data.frame(RPMs=powerLawFittedRPM(values,coef.i))
colnames(RPMs)[1]<-fileName
RPMsTable<-cbind(RPMsTable,RPMs)
}
for(i in 1:length(contrFiles)){
print(paste("Calculating RPMs reads in...",contrFiles[i], "....file in the control"))
fileName<-gsub(".bam","",contrFiles[i])
count.data.f<-mcols(count.data)
count.i<-count.data.f[,colnames(count.data.f)==fileName]
values<-count.i
coef.i<-getPowerLawFittedCoeficient(values)
RPMs<-data.frame(RPMs=powerLawFittedRPM(values,coef.i))
colnames(RPMs)[1]<-fileName
RPMsTable<-cbind(RPMsTable,RPMs)
}
####################################
n<-ncol(RPMsTable)
RPMsTable.GRanges<-GRanges(seqnames=RPMsTable$chromosome,
IRanges(start=RPMsTable$start,end=RPMsTable$end))
values(RPMsTable.GRanges)<-RPMsTable[,4:n]
RPMsTable(object)<-RPMsTable.GRanges
assign(objName,object,envir=parent.frame())
invisible(1)
print("The RPMs table is created!!!.")
}
if(selected_methods=="normalRPM"){
RPMsTable<-data.frame(chromosome=seqnames(count.data),start=start(count.data),end=end(count.data))
for(i in 1:length(expFiles)){
print(paste("Calculating RPMs reads in...",expFiles[i], "....file in the experiment"))
fileName<-gsub(".bam","",expFiles[i])
count.data.f<-mcols(count.data)
count.i<-count.data.f[,colnames(count.data.f)==fileName]
values<-count.i
RPMs<-data.frame(RPMs=normalcalRPM(values,expLibs[i]))
colnames(RPMs)[1]<-fileName
RPMsTable<-cbind(RPMsTable,RPMs)
}
for(i in 1:length(contrFiles)){
print(paste("Calculating RPMs reads in...",contrFiles[i], "....file in the control"))
fileName<-gsub(".bam","",contrFiles[i])
count.data.f<-mcols(count.data)
count.i<-count.data.f[,colnames(count.data.f)==fileName]
values<-count.i
RPMs<-data.frame(RPMs=normalcalRPM(values,contrLibs[i]))
colnames(RPMs)[1]<-fileName
RPMsTable<-cbind(RPMsTable,RPMs)
}
####################################
n<-ncol(RPMsTable)
RPMsTable.GRanges<-GRanges(seqnames=RPMsTable$chromosome,IRanges(start=RPMsTable$start,end=RPMsTable$end))
values(RPMsTable.GRanges)<-RPMsTable[,4:n]
RPMsTable(object)<-RPMsTable.GRanges
assign(objName,object,envir=parent.frame())
invisible(1)
print("The RPMs table is created!!!.")
}
}
)
setGeneric(
name="getBatchInteractions",
def=function(object,method=c("union","intersection"),smoothing.parameter=0.1,fdr=0.05){
standardGeneric("getBatchInteractions")
}
)
setMethod("getBatchInteractions",
signature(object = "r3CseqInBatch"),
function (object,method,smoothing.parameter,fdr){
if(!is(object,"r3CseqInBatch")){
stop("Need the r3CseqInBatch object")
}
#####Select the method#########
selected_methods <- match.arg(method)
if(selected_methods==""){
selected_methods<-"union"
}
if(smoothing.parameter <0.06 | smoothing.parameter>0.4){
stop("The smoothing parameter must be 0.06-0.4 (default=0.1).")
}
##########################################
objName <- deparse(substitute(object))
##########get labels#######################
expFiles<-BamExpFiles(object)
contrFiles<-BamContrFiles(object)
##########get read count##################
count.data<-readCountTable(object)
if(length(count.data)==0){
stop("No reads found in the countTable object, please run getBatchReadCountPerRestrictionFragment or getBatchReadCountPerWindow function!!")
}
##########get RPM##################
rpm.data<-RPMsTable(object)
if(length(rpm.data)==0){
stop("No RPMs found, please run calculateBatchRPM function!!")
}
#############get interaction in experiments###########
expInteraction<-data.frame(chromosome=seqnames(rpm.data),start=start(rpm.data),end=end(rpm.data))
for(i in 1:length(expFiles)){
print(paste("Analyzing interaction regions ...",expFiles[i], ".... in the experiment"))
fileName<-gsub(".bam","",expFiles[i])
count.data.f<-mcols(count.data)
count.i<-count.data.f[,colnames(count.data.f)==fileName]
rpm.data.f<-mcols(rpm.data)
rpm.i<-rpm.data.f[,colnames(rpm.data.f)==fileName]
count.i<-GRanges(seqnames = seqnames(count.data),IRanges(start=start(count.data),
end=end(count.data)),nReads=count.i,RPMs=rpm.i)
count.i<-count.i[count.i$RPMs>1,]
expInt.i<-assign3CseqSigContact(object,count.i,smoothing.parameter,fdr)
expInt.i.frame<-data.frame(chromosome=seqnames(expInt.i),
start=start(expInt.i),
end=end(expInt.i),
nReads=expInt.i$nReads,
RPMs=expInt.i$RPMs,
p.value=expInt.i$p.value,
q.value=expInt.i$q.value)
expInteraction<-merge(expInteraction,expInt.i.frame,by.x=c("chromosome","start","end"),by.y=c("chromosome","start","end"),all=T)
expInt.i.frame<-data.frame(chromosome=seqnames(expInt.i),
start=start(expInt.i),
end=end(expInt.i),
nReads=expInt.i$nReads,
RPMs=expInt.i$RPMs,
p.value=expInt.i$p.value,
q.value=expInt.i$q.value)
write.table(expInt.i.frame,file=paste(fileName,".interaction.txt",sep=""),append=FALSE, sep="\t", quote=FALSE,row.names=FALSE, col.names=TRUE)
}
n.col<-ncol(expInteraction)
expInteraction.filtered<-expInteraction[rowSums(is.na(expInteraction[,4:n.col])) !=n.col-3,]
#############get interaction in controls###############
contrInteraction<-data.frame(chromosome=seqnames(rpm.data),start=start(rpm.data),end=end(rpm.data))
for(i in 1:length(contrFiles)){
print(paste("Analyzing interaction regions ...",contrFiles[i], ".... in the control"))
fileName<-gsub(".bam","",contrFiles[i])
count.data.f<-mcols(count.data)
count.i<-count.data.f[,colnames(count.data.f)==fileName]
rpm.data.f<-mcols(rpm.data)
rpm.i<-rpm.data.f[,colnames(rpm.data.f)==fileName]
count.i<-GRanges(seqnames = seqnames(count.data),IRanges(start=start(count.data),
end=end(count.data)),nReads=count.i,RPMs=rpm.i)
count.i<-count.i[count.i$RPMs>1,]
contrInt.i<-assign3CseqSigContact(object,count.i,smoothing.parameter,fdr)
contrInt.i.frame<-data.frame(chromosome=seqnames(contrInt.i),
start=start(contrInt.i),
end=end(contrInt.i),
nReads=contrInt.i$nReads,
RPMs=contrInt.i$RPMs,
p.value=contrInt.i$p.value,
q.value=contrInt.i$q.value)
contrInteraction<-merge(contrInteraction,contrInt.i.frame,by.x=c("chromosome","start","end"),by.y=c("chromosome","start","end"),all=T)
contrInt.i.frame<-data.frame(chromosome=seqnames(contrInt.i),
start=start(contrInt.i),
end=end(contrInt.i),
nReads=contrInt.i$nReads,
RPMs=contrInt.i$RPMs,
p.value=contrInt.i$p.value,
q.value=contrInt.i$q.value)
write.table(contrInt.i.frame,file=paste(fileName,".interaction.txt",sep=""),append=FALSE, sep="\t", quote=FALSE,row.names=FALSE, col.names=TRUE)
}
m.col<-ncol(contrInteraction)
contrInteraction.filtered<-contrInteraction[rowSums(is.na(contrInteraction[,4:m.col])) !=m.col-3,]
if(selected_methods=="intersection"){
##############################################
exp.intersec<-na.omit(expInteraction.filtered)
exp.readCounts<-exp.intersec[,c(seq(from=4,to=n.col,by=4))]
exp.RPMs<-exp.intersec[,c(seq(from=5,to=n.col,by=4))]
exp.pvalue<-exp.intersec[,c(seq(from=6,to=n.col,by=4))]
exp.nReads.mean<-round(rowMeans(exp.readCounts))
exp.RPMs.mean<-rowMeans(exp.RPMs)
exp.pvalue[exp.pvalue == 0] <- 1e-20
exp.fisher.pvalue<-apply(exp.pvalue,1,fishersMethod)
exp.fisher.qvalue<-qvalue(exp.fisher.pvalue, fdr.level=fdr, pi0.method="bootstrap")$qvalues
exp.intersec.interaction<-GRanges(seqnames=exp.intersec$chromosome,
IRanges(start=exp.intersec$start,end=exp.intersec$end),
nReads=exp.nReads.mean,RPMs=exp.RPMs.mean,p.value=exp.fisher.pvalue,q.value=exp.fisher.qvalue)
##############################################
contr.intersec<-na.omit(contrInteraction.filtered)
contr.readCounts<-contr.intersec[,c(seq(from=4,to=n.col,by=4))]
contr.RPMs<-contr.intersec[,c(seq(from=5,to=n.col,by=4))]
contr.pvalue<-contr.intersec[,c(seq(from=6,to=n.col,by=4))]
contr.nReads.mean<-round(rowMeans(contr.readCounts))
contr.RPMs.mean<-rowMeans(contr.RPMs)
contr.pvalue[contr.pvalue == 0] <- 1e-20
contr.fisher.pvalue<-apply(contr.pvalue,1,fishersMethod)
contr.fisher.qvalue<-qvalue(contr.fisher.pvalue, fdr.level=fdr, pi0.method="bootstrap")$qvalues
contr.intersec.interaction<-GRanges(seqnames=contr.intersec$chromosome,
IRanges(start=contr.intersec$start,end=contr.intersec$end),
nReads=contr.nReads.mean,RPMs=contr.RPMs.mean,
p.value=contr.fisher.pvalue,q.value=contr.fisher.qvalue)
##############################################
expInteractionRegions(object)<-exp.intersec.interaction
contrInteractionRegions(object)<-contr.intersec.interaction
assign(objName,object,envir=parent.frame())
invisible(1)
print("Interactions from the intersection method are created!!!.")
}
if(selected_methods=="union"){
##############################################
exp.union<-expInteraction.filtered
exp.readCounts<-exp.union[,c(seq(from=4,to=n.col,by=4))]
exp.readCounts[is.na(exp.readCounts)==TRUE]<-0
exp.RPMs<-exp.union[,c(seq(from=5,to=n.col,by=4))]
exp.RPMs[is.na(exp.RPMs)==TRUE]<-0
exp.pvalue<-exp.union[,c(seq(from=6,to=n.col,by=4))]
exp.pvalue[is.na(exp.pvalue)==TRUE]<-1
exp.nReads.mean<-round(rowMeans(exp.readCounts))
exp.RPMs.mean<-rowMeans(exp.RPMs)
exp.pvalue[exp.pvalue == 0] <- 1e-20
exp.fisher.pvalue<-apply(exp.pvalue,1,fishersMethod)
exp.fisher.qvalue<-qvalue(exp.fisher.pvalue, fdr.level=fdr, pi0.method="bootstrap")$qvalues
exp.union.interaction<-GRanges(seqnames=exp.union$chromosome,
IRanges(start=exp.union$start,end=exp.union$end),
nReads=exp.nReads.mean,RPMs=exp.RPMs.mean,
p.value=exp.fisher.pvalue,q.value=exp.fisher.qvalue)
##############################################
contr.union<-contrInteraction.filtered
contr.readCounts<-contr.union[,c(seq(from=4,to=n.col,by=4))]
contr.readCounts[is.na(contr.readCounts)==TRUE]<-0
contr.RPMs<-contr.union[,c(seq(from=5,to=n.col,by=4))]
contr.RPMs[is.na(contr.RPMs)==TRUE]<-0
contr.pvalue<-contr.union[,c(seq(from=6,to=n.col,by=4))]
contr.pvalue[is.na(contr.pvalue)==TRUE]<-1
contr.nReads.mean<-round(rowMeans(contr.readCounts))
contr.RPMs.mean<-rowMeans(contr.RPMs)
contr.pvalue[contr.pvalue == 0] <- 1e-20
contr.fisher.pvalue<-apply(contr.pvalue,1,fishersMethod)
contr.fisher.qvalue<-qvalue(contr.fisher.pvalue, fdr.level=fdr, pi0.method="bootstrap")$qvalues
contr.union.interaction<-GRanges(seqnames=contr.union$chromosome,
IRanges(start=contr.union$start,end=contr.union$end),
nReads=contr.nReads.mean,RPMs=contr.RPMs.mean,
p.value=contr.fisher.pvalue,q.value=contr.fisher.qvalue)
##############################################
expInteractionRegions(object)<-exp.union.interaction
contrInteractionRegions(object)<-contr.union.interaction
assign(objName,object,envir=parent.frame())
invisible(1)
print("Interactions from the union method are created!!!.")
}
}
)
supatt-lab/r3Cseq documentation built on May 6, 2019, 4:34 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# TODO: These following functions were implemented to get the interaction regions for replicates data.
# Author: Supat Thongjuea
# Contact : supat.thongjuea@imm.ox.ac.uk or supat.thongjuea@gmail.com
#####
#####The functions below are completely migrated to BioC 3.9 on R 3.6
#####
########################################
setGeneric(
name="getBatchRawReads",
def=function(object){
standardGeneric("getBatchRawReads")
}
)
setMethod("getBatchRawReads",
signature(object = "r3CseqInBatch"),
function (object){
if(!is(object,"r3CseqInBatch")){
stop("Need to initialize the r3CseqInBatch object")
}
######check files########################
bamDir<-bamFilesDirectory(object)
##########get labels#######################
expFiles<-BamExpFiles(object)
contrFiles<-BamContrFiles(object)
#########checking directory and files####
if (file.exists(bamDir)){
setwd(file.path(bamDir))
}else{
stop(paste("Couldn't find path","-->",bamDir))
}
objName <- deparse(substitute(object))
exp.bam.files <-BamExpFiles(object)
contr.bam.files <-BamContrFiles(object)
for(i.file in c(exp.bam.files)){
if(file.exists(i.file) ==FALSE){
stop(paste("Couldn't find file","-->",i.file))
}
}
for(i.file in c(contr.bam.files)){
if(file.exists(i.file) ==FALSE){
stop(paste("Couldn't find file","-->",i.file))
}
}
exp.lib.size<-c()
exp.read.length<-c()
for(i in 1:length(exp.bam.files)){
input.bam<-exp.bam.files[i]
print(paste("start reading in...",input.bam, "....file in the experiment"))
what <- c("rname", "strand", "pos", "qwidth","seq")
param <- ScanBamParam(what=what,flag = scanBamFlag(isUnmappedQuery = FALSE))
exp.bam <- scanBam(input.bam, param=param)
exp.length <-width(exp.bam[[1]]$seq[1])
exp.GRanges<-GRanges(seqnames=as.vector(exp.bam[[1]]$rname),IRanges(start=exp.bam[[1]]$pos,width=exp.length),strand=exp.bam[[1]]$strand)
fileName<-gsub(".bam","",expFiles[i])
lib.size <-length(exp.GRanges)
exp.lib.size<-append(exp.lib.size,lib.size)
exp.read.length<-append(exp.read.length,exp.length)
save(exp.GRanges,file=paste(fileName,".rData",sep=""))
}
contr.lib.size<-c()
contr.read.length<-c()
###########Make count table for each bam file in the controls########
for(i in 1:length(contr.bam.files)){
input.bam<-contr.bam.files[i]
print(paste("start reading in...",input.bam, "....file in the control"))
what <- c("rname", "strand", "pos", "qwidth","seq")
param <- ScanBamParam(what=what,flag = scanBamFlag(isUnmappedQuery = FALSE))
contr.bam <- scanBam(input.bam, param=param)
contr.length <-width(contr.bam[[1]]$seq[1])
contr.GRanges<-GRanges(seqnames=as.vector(contr.bam[[1]]$rname),IRanges(start=contr.bam[[1]]$pos,width=contr.length),strand=contr.bam[[1]]$strand)
fileName<-gsub(".bam","",contrFiles[i])
lib.size <-length(contr.GRanges)
contr.lib.size<-append(contr.lib.size,lib.size)
contr.read.length<-append(contr.read.length,contr.length)
save(contr.GRanges,file=paste(fileName,".rData",sep=""))
}
expBatchLibrarySize(object)<-exp.lib.size
contrBatchLibrarySize(object)<-contr.lib.size
expBatchReadLength(object)<-exp.read.length
contrBatchReadLength(object)<-contr.read.length
assign(objName,object,envir=parent.frame())
invisible(1)
print("The raw read files are created!!!.")
}
)
setGeneric(
name="getBatchReadCountPerRestrictionFragment",
def=function(object,getReadsMethod = c("wholeReads", "adjacentFragmentEndsReads"),
nFragmentExcludedReadsNearViewpoint=2){
standardGeneric("getBatchReadCountPerRestrictionFragment")
}
)
setMethod("getBatchReadCountPerRestrictionFragment",
signature(object = "r3CseqInBatch"),
function (object,getReadsMethod,nFragmentExcludedReadsNearViewpoint){
if(!is(object,"r3CseqInBatch")){
stop("Need to initialize the r3CseqInBatch object")
}
selected_methods <- match.arg(getReadsMethod)
if(selected_methods==""){
selected_methods<-"wholeReads"
}
######check the number of fragment to exclude reads around the viewpoint###
nExcludedFragments <- nFragmentExcludedReadsNearViewpoint
if(nExcludedFragments <0 | nExcludedFragments>5){
stop("The number of fragment that uses for removing reads around the viewpoint must be 0-5 restriction fragments.")
}
######check files########################
objName <- deparse(substitute(object))
bamDir<-bamFilesDirectory(object)
#########checking directory and files####
if (file.exists(bamDir)){
setwd(file.path(bamDir))
}else{
stop(paste("Couldn't find path","-->",bamDir))
}
exp.bam.files <-BamExpFiles(object)
contr.bam.files <-BamContrFiles(object)
for(i.file in c(exp.bam.files)){
i.file<-gsub(".bam","",i.file)
i.file<-paste(i.file,".rData",sep="")
if(file.exists(i.file) ==FALSE){
stop(paste("Couldn't find file","-->",i.file, ". Please run getBatchRawReads function"))
}
}
for(i.file in c(contr.bam.files)){
i.file<-gsub(".bam","",i.file)
i.file<-paste(i.file,".rData",sep="")
if(file.exists(i.file) ==FALSE){
stop(paste("Couldn't find file","-->",i.file, ". Please run getBatchRawReads function"))
}
}
####Get organism#######################
orgName <- organismName(object)
####Initialized Restriction Enzyme Object##
enzymeDb <-new("repbaseEnzyme")
resEnzyme <-restrictionEnzyme(object)
####Build GRanges for restriction site#####
print("Fragmenting genome..........")
wholeFragments<-getWholeGenomeRestrictionFragments(enzymeDb,resEnzyme,orgName)
##########get labels#######################
expFiles<-BamExpFiles(object)
contrFiles<-BamContrFiles(object)
###########################################
exp.read.length<-expBatchReadLength(object)
contr.read.length<-contrBatchReadLength(object)
if(selected_methods=="wholeReads"){
###########Making GRange Object############
wholeFragments.GRanges<-GRanges(seqnames = as.character(wholeFragments$chromosome),IRanges(start=wholeFragments$start,end=wholeFragments$end))
###########Make count table for each bam file in the experiments########
readCountTable<-data.frame(chromosome=seqnames(wholeFragments.GRanges),start=start(wholeFragments.GRanges),end=end(wholeFragments.GRanges))
for(i in 1:length(exp.bam.files)){
print(paste("start counting reads in...",exp.bam.files[i], "....file in the experiment"))
fileName<-gsub(".bam","",expFiles[i])
fileName<-paste(fileName,".rData",sep="")
load(fileName)
exp.GRanges<-excludeReadsNearViewpoint(object,exp.GRanges,nExcludedFragments)
readsPerFragment <- countOverlaps(wholeFragments.GRanges,subject=exp.GRanges)
reads<-data.frame(nReads=readsPerFragment)
Name<-gsub(".bam","",expFiles[i])
colnames(reads)[1]<-Name
readCountTable<-cbind(readCountTable,reads)
}
###########Make count table for each bam file in the controls########
for(i in 1:length(contr.bam.files)){
print(paste("start counting reads in...",contr.bam.files[i], "....file in the control"))
fileName<-gsub(".bam","",contrFiles[i])
fileName<-paste(fileName,".rData",sep="")
load(fileName)
contr.GRanges<-excludeReadsNearViewpoint(object,contr.GRanges,nExcludedFragments)
readsPerFragment <- countOverlaps(wholeFragments.GRanges,subject=contr.GRanges)
reads<-data.frame(nReads=readsPerFragment)
Name<-gsub(".bam","",contrFiles[i])
colnames(reads)[1]<-Name
readCountTable<-cbind(readCountTable,reads)
}
#########Filter out fragment with no read###########################
n<-ncol(readCountTable)
readCountTable.f<-readCountTable[rowSums(readCountTable[,4:n]) >0,]
readCountTable.GRanges<-GRanges(seqnames=readCountTable.f$chromosome,
IRanges(start=readCountTable.f$start,end=readCountTable.f$end))
values(readCountTable.GRanges)<-readCountTable.f[,4:n]
readCountTable(object)<-readCountTable.GRanges
assign(objName,object,envir=parent.frame())
invisible(1)
print("The countTable is created!!!.")
}
if(selected_methods=="adjacentFragmentEndsReads"){
###########Making GRange Object############
wholeFragments.GRanges<-GRanges(seqnames=as.character(wholeFragments$chromosome),
IRanges(start=wholeFragments$start,end=wholeFragments$end))
###########Make count table for each bam file in the experiments########
readCountTable<-data.frame(chromosome=seqnames(wholeFragments.GRanges),
start=start(wholeFragments.GRanges),end=end(wholeFragments.GRanges))
for(i in 1:length(exp.bam.files)){
print(paste("start counting reads in...",exp.bam.files[i], "....file in the experiment"))
fileName<-gsub(".bam","",expFiles[i])
fileName<-paste(fileName,".rData",sep="")
load(fileName)
exp.GRanges<-excludeReadsNearViewpoint(object,exp.GRanges,nExcludedFragments)
########Make RangedData for the 5' and 3' of RE#######
wholeFragments_5_prime.GRanges<-GRanges(seqnames=seqnames(wholeFragments.GRanges),
IRanges(start=start(wholeFragments.GRanges),
end=start(wholeFragments.GRanges)+exp.read.length[i]))
wholeFragments_3_prime.GRanges<-GRanges(seqnames=seqnames(wholeFragments.GRanges),
IRanges(start=end(wholeFragments.GRanges)-exp.read.length[i],
end=end(wholeFragments.GRanges)))
readsPerFragment_5 <- countOverlaps(wholeFragments_5_prime.GRanges,subject=exp.GRanges)
readsPerFragment_3 <- countOverlaps(wholeFragments_3_prime.GRanges,subject=exp.GRanges)
combined_reads<-readsPerFragment_5+readsPerFragment_3
reads<-data.frame(nReads=combined_reads)
fileName<-gsub(".bam","",expFiles[i])
colnames(reads)[1]<-fileName
readCountTable<-cbind(readCountTable,reads)
}
for(i in 1:length(contr.bam.files)){
print(paste("start counting reads in...",contr.bam.files[i], "....file in the control"))
fileName<-gsub(".bam","",contrFiles[i])
fileName<-paste(fileName,".rData",sep="")
load(fileName)
contr.GRanges<-excludeReadsNearViewpoint(object,contr.GRanges,nExcludedFragments)
########Make RangedData for the 5' and 3' of RE#######
wholeFragments_5_prime.GRanges<-GRanges(seqnames=seqnames(wholeFragments.GRanges),
IRanges(start=start(wholeFragments.GRanges),
end=start(wholeFragments.GRanges)+contr.read.length[i]))
wholeFragments_3_prime.GRanges<-GRanges(seqnames=seqnames(wholeFragments.GRanges),
IRanges(start=end(wholeFragments.GRanges)-contr.read.length[i],
end=end(wholeFragments.GRanges)))
readsPerFragment_5 <- countOverlaps(wholeFragments_5_prime.GRanges,subject=contr.GRanges)
readsPerFragment_3 <- countOverlaps(wholeFragments_3_prime.GRanges,subject=contr.GRanges)
combined_reads<-readsPerFragment_5+readsPerFragment_3
reads<-data.frame(nReads=combined_reads)
fileName<-gsub(".bam","",contrFiles[i])
colnames(reads)[1]<-fileName
readCountTable<-cbind(readCountTable,reads)
}
#########Filter out fragment with no read###########################
n<-ncol(readCountTable)
readCountTable.f<-readCountTable[rowSums(readCountTable[,4:n]) >0,]
readCountTable.GRanges<-GRanges(seqnames=readCountTable.f$chromosome,
IRanges(start=readCountTable.f$start,end=readCountTable.f$end))
values(readCountTable.GRanges)<-readCountTable.f[,4:n]
readCountTable(object)<-readCountTable.GRanges
assign(objName,object,envir=parent.frame())
invisible(1)
print("The countTable is created!!!.")
}
}
)
setGeneric(
name="getBatchReadCountPerWindow",
def=function(object,windowSize=5e3,nFragmentExcludedReadsNearViewpoint=2,mode=c("non-overlapping", "overlapping")){
standardGeneric("getBatchReadCountPerWindow")
}
)
setMethod("getBatchReadCountPerWindow",
signature(object = "r3CseqInBatch"),
function (object,windowSize,nFragmentExcludedReadsNearViewpoint,mode){
objName <- deparse(substitute(object))
if(!is(object,"r3CseqInBatch")){
stop("Need to initialize the r3CseqInBatch object")
}
######check the number of fragment to exclude reads around the viewpoint###
nExcludedFragments <- nFragmentExcludedReadsNearViewpoint
if(nExcludedFragments <0 | nExcludedFragments>5){
stop("The number of fragment that uses for removing reads around the viewpoint must be 0-5 restriction fragments.")
}
#####Select the read count methods#########
mode <- match.arg(mode)
if(mode==""){
mode<-"non-overlapping"
}
######check files########################
objName <- deparse(substitute(object))
bamDir<-bamFilesDirectory(object)
#########checking directory and files####
if (file.exists(bamDir)){
setwd(file.path(bamDir))
}else{
stop(paste("Couldn't find path","-->",bamDir))
}
exp.bam.files <-BamExpFiles(object)
contr.bam.files <-BamContrFiles(object)
for(i.file in c(exp.bam.files)){
i.file<-gsub(".bam","",i.file)
i.file<-paste(i.file,".rData",sep="")
if(file.exists(i.file) ==FALSE){
stop(paste("Couldn't find file","-->",i.file, ". Please run getBatchRawReads function"))
}
}
for(i.file in c(contr.bam.files)){
i.file<-gsub(".bam","",i.file)
i.file<-paste(i.file,".rData",sep="")
if(file.exists(i.file) ==FALSE){
stop(paste("Couldn't find file","-->",i.file, ". Please run getBatchRawReads function"))
}
}
##########get labels#######################
expFiles<-BamExpFiles(object)
contrFiles<-BamContrFiles(object)
###########################################
wholeFragments.GRanges<-getFragmentsPerWindow(object,windowSize,mode)
###########Make count table for each bam file in the experiments########
readCountTable<-data.frame(chromosome=seqnames(wholeFragments.GRanges),
start=start(wholeFragments.GRanges),
end=end(wholeFragments.GRanges))
for(i in 1:length(exp.bam.files)){
print(paste("start counting reads in...",exp.bam.files[i], "....file in the experiment"))
fileName<-gsub(".bam","",expFiles[i])
fileName<-paste(fileName,".rData",sep="")
load(fileName)
exp.GRanges<-excludeReadsNearViewpoint(object,exp.GRanges,nExcludedFragments)
readsPerFragment <- countOverlaps(wholeFragments.GRanges,subject=exp.GRanges)
reads<-data.frame(nReads=readsPerFragment)
Name<-gsub(".bam","",expFiles[i])
colnames(reads)[1]<-Name
readCountTable<-cbind(readCountTable,reads)
}
###########Make count table for each bam file in the controls########
for(i in 1:length(contr.bam.files)){
print(paste("start counting reads in...",contr.bam.files[i], "....file in the control"))
fileName<-gsub(".bam","",contrFiles[i])
fileName<-paste(fileName,".rData",sep="")
load(fileName)
contr.GRanges<-excludeReadsNearViewpoint(object,contr.GRanges,nExcludedFragments)
readsPerFragment <- countOverlaps(wholeFragments.GRanges,subject=contr.GRanges)
reads<-data.frame(nReads=readsPerFragment)
Name<-gsub(".bam","",contrFiles[i])
colnames(reads)[1]<-Name
readCountTable<-cbind(readCountTable,reads)
}
#########Filter out fragment with no read###########################
n<-ncol(readCountTable)
readCountTable.f<-readCountTable[rowSums(readCountTable[,4:n]) >0,]
readCountTable.GRanges<-GRanges(seqnames=readCountTable.f$chromosome,
IRanges(start=readCountTable.f$start,
end=readCountTable.f$end))
values(readCountTable.GRanges)<-readCountTable.f[,4:n]
readCountTable(object)<-readCountTable.GRanges
assign(objName,object,envir=parent.frame())
invisible(1)
print("The countTable is created!!!.")
}
)
setGeneric(
name="calculateBatchRPM",
def=function(object,normalized_method=c("powerlawFittedRPM","normalRPM")){
standardGeneric("calculateBatchRPM")
}
)
setMethod("calculateBatchRPM",
signature(object = "r3CseqInBatch"),
function (object,normalized_method){
if(!is(object,"r3CseqInBatch")){
stop("Need the r3CseqInBatch object")
}
#####Select the method#########
selected_methods <- match.arg(normalized_method)
if(selected_methods==""){
selected_methods<-"powerlawFittedRPM"
}
##########get labels#######################
expFiles<-BamExpFiles(object)
contrFiles<-BamContrFiles(object)
##########Get library size#################
expLibs<-expBatchLibrarySize(object)
contrLibs<-contrBatchLibrarySize(object)
###########################################
objName <- deparse(substitute(object))
count.data<-readCountTable(object)
if(length(count.data)==0){
stop("No reads found in the countTable object, please run getBatchReadCountPerRestrictionFragment or getBatchReadCountPerWindow function!!")
}
if(selected_methods=="powerlawFittedRPM"){
RPMsTable<-data.frame(chromosome=seqnames(count.data),start=start(count.data),end=end(count.data))
for(i in 1:length(expFiles)){
print(paste("Calculating RPMs reads in...",expFiles[i], "....file in the experiment"))
fileName<-gsub(".bam","",expFiles[i])
count.data.f<-mcols(count.data)
count.i<-count.data.f[,colnames(count.data.f)==fileName]
values<-count.i
coef.i<-getPowerLawFittedCoeficient(values)
RPMs<-data.frame(RPMs=powerLawFittedRPM(values,coef.i))
colnames(RPMs)[1]<-fileName
RPMsTable<-cbind(RPMsTable,RPMs)
}
for(i in 1:length(contrFiles)){
print(paste("Calculating RPMs reads in...",contrFiles[i], "....file in the control"))
fileName<-gsub(".bam","",contrFiles[i])
count.data.f<-mcols(count.data)
count.i<-count.data.f[,colnames(count.data.f)==fileName]
values<-count.i
coef.i<-getPowerLawFittedCoeficient(values)
RPMs<-data.frame(RPMs=powerLawFittedRPM(values,coef.i))
colnames(RPMs)[1]<-fileName
RPMsTable<-cbind(RPMsTable,RPMs)
}
####################################
n<-ncol(RPMsTable)
RPMsTable.GRanges<-GRanges(seqnames=RPMsTable$chromosome,
IRanges(start=RPMsTable$start,end=RPMsTable$end))
values(RPMsTable.GRanges)<-RPMsTable[,4:n]
RPMsTable(object)<-RPMsTable.GRanges
assign(objName,object,envir=parent.frame())
invisible(1)
print("The RPMs table is created!!!.")
}
if(selected_methods=="normalRPM"){
RPMsTable<-data.frame(chromosome=seqnames(count.data),start=start(count.data),end=end(count.data))
for(i in 1:length(expFiles)){
print(paste("Calculating RPMs reads in...",expFiles[i], "....file in the experiment"))
fileName<-gsub(".bam","",expFiles[i])
count.data.f<-mcols(count.data)
count.i<-count.data.f[,colnames(count.data.f)==fileName]
values<-count.i
RPMs<-data.frame(RPMs=normalcalRPM(values,expLibs[i]))
colnames(RPMs)[1]<-fileName
RPMsTable<-cbind(RPMsTable,RPMs)
}
for(i in 1:length(contrFiles)){
print(paste("Calculating RPMs reads in...",contrFiles[i], "....file in the control"))
fileName<-gsub(".bam","",contrFiles[i])
count.data.f<-mcols(count.data)
count.i<-count.data.f[,colnames(count.data.f)==fileName]
values<-count.i
RPMs<-data.frame(RPMs=normalcalRPM(values,contrLibs[i]))
colnames(RPMs)[1]<-fileName
RPMsTable<-cbind(RPMsTable,RPMs)
}
####################################
n<-ncol(RPMsTable)
RPMsTable.GRanges<-GRanges(seqnames=RPMsTable$chromosome,IRanges(start=RPMsTable$start,end=RPMsTable$end))
values(RPMsTable.GRanges)<-RPMsTable[,4:n]
RPMsTable(object)<-RPMsTable.GRanges
assign(objName,object,envir=parent.frame())
invisible(1)
print("The RPMs table is created!!!.")
}
}
)
setGeneric(
name="getBatchInteractions",
def=function(object,method=c("union","intersection"),smoothing.parameter=0.1,fdr=0.05){
standardGeneric("getBatchInteractions")
}
)
setMethod("getBatchInteractions",
signature(object = "r3CseqInBatch"),
function (object,method,smoothing.parameter,fdr){
if(!is(object,"r3CseqInBatch")){
stop("Need the r3CseqInBatch object")
}
#####Select the method#########
selected_methods <- match.arg(method)
if(selected_methods==""){
selected_methods<-"union"
}
if(smoothing.parameter <0.06 | smoothing.parameter>0.4){
stop("The smoothing parameter must be 0.06-0.4 (default=0.1).")
}
##########################################
objName <- deparse(substitute(object))
##########get labels#######################
expFiles<-BamExpFiles(object)
contrFiles<-BamContrFiles(object)
##########get read count##################
count.data<-readCountTable(object)
if(length(count.data)==0){
stop("No reads found in the countTable object, please run getBatchReadCountPerRestrictionFragment or getBatchReadCountPerWindow function!!")
}
##########get RPM##################
rpm.data<-RPMsTable(object)
if(length(rpm.data)==0){
stop("No RPMs found, please run calculateBatchRPM function!!")
}
#############get interaction in experiments###########
expInteraction<-data.frame(chromosome=seqnames(rpm.data),start=start(rpm.data),end=end(rpm.data))
for(i in 1:length(expFiles)){
print(paste("Analyzing interaction regions ...",expFiles[i], ".... in the experiment"))
fileName<-gsub(".bam","",expFiles[i])
count.data.f<-mcols(count.data)
count.i<-count.data.f[,colnames(count.data.f)==fileName]
rpm.data.f<-mcols(rpm.data)
rpm.i<-rpm.data.f[,colnames(rpm.data.f)==fileName]
count.i<-GRanges(seqnames = seqnames(count.data),IRanges(start=start(count.data),
end=end(count.data)),nReads=count.i,RPMs=rpm.i)
count.i<-count.i[count.i$RPMs>1,]
expInt.i<-assign3CseqSigContact(object,count.i,smoothing.parameter,fdr)
expInt.i.frame<-data.frame(chromosome=seqnames(expInt.i),
start=start(expInt.i),
end=end(expInt.i),
nReads=expInt.i$nReads,
RPMs=expInt.i$RPMs,
p.value=expInt.i$p.value,
q.value=expInt.i$q.value)
expInteraction<-merge(expInteraction,expInt.i.frame,by.x=c("chromosome","start","end"),by.y=c("chromosome","start","end"),all=T)
expInt.i.frame<-data.frame(chromosome=seqnames(expInt.i),
start=start(expInt.i),
end=end(expInt.i),
nReads=expInt.i$nReads,
RPMs=expInt.i$RPMs,
p.value=expInt.i$p.value,
q.value=expInt.i$q.value)
write.table(expInt.i.frame,file=paste(fileName,".interaction.txt",sep=""),append=FALSE, sep="\t", quote=FALSE,row.names=FALSE, col.names=TRUE)
}
n.col<-ncol(expInteraction)
expInteraction.filtered<-expInteraction[rowSums(is.na(expInteraction[,4:n.col])) !=n.col-3,]
#############get interaction in controls###############
contrInteraction<-data.frame(chromosome=seqnames(rpm.data),start=start(rpm.data),end=end(rpm.data))
for(i in 1:length(contrFiles)){
print(paste("Analyzing interaction regions ...",contrFiles[i], ".... in the control"))
fileName<-gsub(".bam","",contrFiles[i])
count.data.f<-mcols(count.data)
count.i<-count.data.f[,colnames(count.data.f)==fileName]
rpm.data.f<-mcols(rpm.data)
rpm.i<-rpm.data.f[,colnames(rpm.data.f)==fileName]
count.i<-GRanges(seqnames = seqnames(count.data),IRanges(start=start(count.data),
end=end(count.data)),nReads=count.i,RPMs=rpm.i)
count.i<-count.i[count.i$RPMs>1,]
contrInt.i<-assign3CseqSigContact(object,count.i,smoothing.parameter,fdr)
contrInt.i.frame<-data.frame(chromosome=seqnames(contrInt.i),
start=start(contrInt.i),
end=end(contrInt.i),
nReads=contrInt.i$nReads,
RPMs=contrInt.i$RPMs,
p.value=contrInt.i$p.value,
q.value=contrInt.i$q.value)
contrInteraction<-merge(contrInteraction,contrInt.i.frame,by.x=c("chromosome","start","end"),by.y=c("chromosome","start","end"),all=T)
contrInt.i.frame<-data.frame(chromosome=seqnames(contrInt.i),
start=start(contrInt.i),
end=end(contrInt.i),
nReads=contrInt.i$nReads,
RPMs=contrInt.i$RPMs,
p.value=contrInt.i$p.value,
q.value=contrInt.i$q.value)
write.table(contrInt.i.frame,file=paste(fileName,".interaction.txt",sep=""),append=FALSE, sep="\t", quote=FALSE,row.names=FALSE, col.names=TRUE)
}
m.col<-ncol(contrInteraction)
contrInteraction.filtered<-contrInteraction[rowSums(is.na(contrInteraction[,4:m.col])) !=m.col-3,]
if(selected_methods=="intersection"){
##############################################
exp.intersec<-na.omit(expInteraction.filtered)
exp.readCounts<-exp.intersec[,c(seq(from=4,to=n.col,by=4))]
exp.RPMs<-exp.intersec[,c(seq(from=5,to=n.col,by=4))]
exp.pvalue<-exp.intersec[,c(seq(from=6,to=n.col,by=4))]
exp.nReads.mean<-round(rowMeans(exp.readCounts))
exp.RPMs.mean<-rowMeans(exp.RPMs)
exp.pvalue[exp.pvalue == 0] <- 1e-20
exp.fisher.pvalue<-apply(exp.pvalue,1,fishersMethod)
exp.fisher.qvalue<-qvalue(exp.fisher.pvalue, fdr.level=fdr, pi0.method="bootstrap")$qvalues
exp.intersec.interaction<-GRanges(seqnames=exp.intersec$chromosome,
IRanges(start=exp.intersec$start,end=exp.intersec$end),
nReads=exp.nReads.mean,RPMs=exp.RPMs.mean,p.value=exp.fisher.pvalue,q.value=exp.fisher.qvalue)
##############################################
contr.intersec<-na.omit(contrInteraction.filtered)
contr.readCounts<-contr.intersec[,c(seq(from=4,to=n.col,by=4))]
contr.RPMs<-contr.intersec[,c(seq(from=5,to=n.col,by=4))]
contr.pvalue<-contr.intersec[,c(seq(from=6,to=n.col,by=4))]
contr.nReads.mean<-round(rowMeans(contr.readCounts))
contr.RPMs.mean<-rowMeans(contr.RPMs)
contr.pvalue[contr.pvalue == 0] <- 1e-20
contr.fisher.pvalue<-apply(contr.pvalue,1,fishersMethod)
contr.fisher.qvalue<-qvalue(contr.fisher.pvalue, fdr.level=fdr, pi0.method="bootstrap")$qvalues
contr.intersec.interaction<-GRanges(seqnames=contr.intersec$chromosome,
IRanges(start=contr.intersec$start,end=contr.intersec$end),
nReads=contr.nReads.mean,RPMs=contr.RPMs.mean,
p.value=contr.fisher.pvalue,q.value=contr.fisher.qvalue)
##############################################
expInteractionRegions(object)<-exp.intersec.interaction
contrInteractionRegions(object)<-contr.intersec.interaction
assign(objName,object,envir=parent.frame())
invisible(1)
print("Interactions from the intersection method are created!!!.")
}
if(selected_methods=="union"){
##############################################
exp.union<-expInteraction.filtered
exp.readCounts<-exp.union[,c(seq(from=4,to=n.col,by=4))]
exp.readCounts[is.na(exp.readCounts)==TRUE]<-0
exp.RPMs<-exp.union[,c(seq(from=5,to=n.col,by=4))]
exp.RPMs[is.na(exp.RPMs)==TRUE]<-0
exp.pvalue<-exp.union[,c(seq(from=6,to=n.col,by=4))]
exp.pvalue[is.na(exp.pvalue)==TRUE]<-1
exp.nReads.mean<-round(rowMeans(exp.readCounts))
exp.RPMs.mean<-rowMeans(exp.RPMs)
exp.pvalue[exp.pvalue == 0] <- 1e-20
exp.fisher.pvalue<-apply(exp.pvalue,1,fishersMethod)
exp.fisher.qvalue<-qvalue(exp.fisher.pvalue, fdr.level=fdr, pi0.method="bootstrap")$qvalues
exp.union.interaction<-GRanges(seqnames=exp.union$chromosome,
IRanges(start=exp.union$start,end=exp.union$end),
nReads=exp.nReads.mean,RPMs=exp.RPMs.mean,
p.value=exp.fisher.pvalue,q.value=exp.fisher.qvalue)
##############################################
contr.union<-contrInteraction.filtered
contr.readCounts<-contr.union[,c(seq(from=4,to=n.col,by=4))]
contr.readCounts[is.na(contr.readCounts)==TRUE]<-0
contr.RPMs<-contr.union[,c(seq(from=5,to=n.col,by=4))]
contr.RPMs[is.na(contr.RPMs)==TRUE]<-0
contr.pvalue<-contr.union[,c(seq(from=6,to=n.col,by=4))]
contr.pvalue[is.na(contr.pvalue)==TRUE]<-1
contr.nReads.mean<-round(rowMeans(contr.readCounts))
contr.RPMs.mean<-rowMeans(contr.RPMs)
contr.pvalue[contr.pvalue == 0] <- 1e-20
contr.fisher.pvalue<-apply(contr.pvalue,1,fishersMethod)
contr.fisher.qvalue<-qvalue(contr.fisher.pvalue, fdr.level=fdr, pi0.method="bootstrap")$qvalues
contr.union.interaction<-GRanges(seqnames=contr.union$chromosome,
IRanges(start=contr.union$start,end=contr.union$end),
nReads=contr.nReads.mean,RPMs=contr.RPMs.mean,
p.value=contr.fisher.pvalue,q.value=contr.fisher.qvalue)
##############################################
expInteractionRegions(object)<-exp.union.interaction
contrInteractionRegions(object)<-contr.union.interaction
assign(objName,object,envir=parent.frame())
invisible(1)
print("Interactions from the union method are created!!!.")
}
}
)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.