#' @export
#comparing ChIPseq data with the number p-scores of intersected contacts
ChIPSeqCorContacts <- function(Experiments_4C,rearranged_rawData,ChIPseqVScontacts_correlation_plots)
{
#choosing a genome size file
ls_genomes <- system("ls ~/Analyze4C/genomes/",intern=TRUE)
genome_file.name <- readline(prompt=cat("\nChoose the file of the organism you are dealing with:\n",ls_genomes,"",sep="\n"))
genome_sizes <- read.table(paste("~/Analyze4C/genomes/",genome_file.name,sep=""))
numOFchroms <- nrow(genome_sizes)
chr_names <- rep(0,numOFchroms)
for(t in 1:numOFchroms)
{
chr_names[t] <- paste("chr",t,sep="")
}
translocated_flag <- 0
files_flag <- 0 #this marks if the user does or doesn't want to choose more files to use (0 choose more, 1 don't choose more)
z <- 1 #counts how many files are chosen
secs_rmv <- "no" #this says if there were any sections removed
CO_ps_type <- c()
CO_ps <- c()
CO_ps_chroms <- c()
CO_peaks_type <- c()
CO_peaks <- c()
CO_peaks_chroms <- c()
wind_flag <- 0 #this marks if the user wants to change the window sizes
by_winds <- "no"
winds <- "NA"
psORreads <- "pScores" #are p-score or reads correlated with the ChIPseq data
save_flag <- 0 #marks if there was any plot saved, if yes then the details are entered into 'ChIPseqVScontacts_correlation_plots'
#the user decides if to correlate individual RE sites or by windows
choice4 <- as.integer(readline(prompt=cat("\nchoose a method of comparison:\n\n1) intersect ChIPseq peaks with contact bands\n2) create windows, add the values of ChIPseq peaks (p-scores or tags) and contacts, and correlate them\n\n")))
if(choice4 == 2)
{
by_winds <- "yes"
winds <- c()
}
#asking if to compare multiple ChIPseq files against one p-score file or one ChIPseq file against multiple p-score files
choice1 <- as.integer(readline(prompt=cat("\nchoose how you would like to compare the files:\n1) compare one ChIPseq file against multiple p-score files\n2) compare one p-score file against multiple ChIPseq files\n\n")))
if(choice1 == 1)
{
#choosing one peaks file
#choosing multiple p-scores
#getting the peaks file
ls_files_peaks <- system("ls ~/Analyze4C/ChIPseq/peaks | grep 'tagsANDpscores.bed$'",intern=TRUE)
repeat
{
file.name_peaks <- readline(prompt=cat("\nThese are the ChIPseq data (tags or p-scores) files available\nwhich would you like to use?\n",ls_files_peaks,"",sep="\n"))
ind_files_peaks <- pmatch(file.name_peaks,ls_files_peaks,nomatch=0)
if(ind_files_peaks == 0)
{
cat("no such file exists.\nplease try again.\n\n")
}
else
{
#importing the data from the file
peaks <- read.table(paste("~/Analyze4C/ChIPseq/peaks/",file.name_peaks,sep=""))
break
}
#here i might need to edit and sort the data in peaks
#peaks <- peaks[peaks[,1]!="Pt" & peaks[,1]!="Mt",]
#peaks <- peaks[order(peaks[,1],peaks[,2]),]
#names(peaks) <- c("chr","first","last","peaks")
}
peaks_Experiments <- file.name_peaks
num_of_peaks_files <- 1
ps_datName_all <- c() #this will contain all the names the user gives the p-score data for the plot
}
else
{
#choosing one p-score file
repeat
{
choice3 <- as.integer(readline(prompt=cat("\nwhat folder of p-score files would you like to choose from? (enter the number)\n1) no bp windows\n2) with bp windows\n\n")))
if(choice3 == 1) #no bp windows
{
ps_folder <- "no_bp"
}
else if(choice3 == 2) #with bp windows
{
ps_folder <- "with_bp"
}
ps_files <- system(paste("ls ~/Analyze4C/pScores/",ps_folder,sep=""),intern=TRUE)
ps.file <- readline(prompt=cat("These are the p-score files available\nwhich would you like to use?\n",ps_files,"",sep="\n"))
ind_files2 <- pmatch(ps.file,ps_files,nomatch=0)
if(ind_files2 == 0)
{
cat("no such file exists.\nplease try again.\n\n")
}
else
{
#here we check to see if the p-score file chosen is of a translocation switched type
if(identical(grep("translocatedSwitched",ps.file),integer(0)) == "FALSE") #the p-score file is of translocation switched and needs to be removed
{
#the translocation area needs to be removed both from the p-score and from the peaks file
translocated_flag <- 1
}
break
}
}
ps_Experiments <- ps.file
num_of_pScore_files <- 1
#importing the p-score data from the file
ps <- read.table(paste("~/Analyze4C/pScores/",ps_folder,"/",ps.file,sep=""))
#finding the specific p-score files information in Experiments_4C
sp1 <- unlist(strsplit(ps.file,"[.]"))
sp2 <- unlist(strsplit(sp1[1],"_"))
#getting the size of the window per RE sites
sp3 <- sp2[as.integer(grep("RE$",sp2))]
RE_wind <- as.numeric(strsplit(sp3,"RE"))
if(choice3 == 2)
{
#getting the size of the window per bp sites
sp4 <- sp2[as.integer(grep("bp$",sp2))]
bp_wind <- as.numeric(strsplit(sp4,"bp"))
}
#finding the specific p-scores files information in Experiments_4C
out <- findIn.Experiments_4C(sp2[1],sp2[2],sp2[3],Experiments_4C)
#getting the cis chromosome number
cis <- as.numeric(out[2])
peaks_datName_all <- c() #this will contain all the names the user gives the peaks data for the plot
}
intersected_pairs_all <- c() #this will contain all the intersections data from all the files
intersected_reads_pairs_all <- c() #this will contain all the intersections data from all the files, for the reads
while(files_flag == 0)
{
if(choice1 == 1)
{
cat("\np-score file number ",z,":\n",sep="")
#choosing a p-score file:
repeat
{
choice3 <- as.integer(readline(prompt=cat("\nwhat folder of p-score files would you like to choose from? (enter the number)\n1) no bp windows\n2) with bp windows\n\n")))
if(choice3 == 1) #no bp windows
{
ps_folder <- "no_bp"
}
else if(choice3 == 2) #with bp windows
{
ps_folder <- "with_bp"
}
ps_files <- system(paste("ls ~/Analyze4C/pScores/",ps_folder,sep=""),intern=TRUE)
ps.file <- readline(prompt=cat("These are the p-score files available\nwhich would you like to use?\n",ps_files,"",sep="\n"))
ind_files2 <- pmatch(ps.file,ps_files,nomatch=0)
if(ind_files2 == 0)
{
cat("no such file exists.\nplease try again.\n\n")
}
else
{
#here we check to see if the p-score file chosen is of a translocation switched type
if(identical(grep("translocatedSwitched",ps.file),integer(0)) == "FALSE") #the p-score file is of translocation switched and needs to be removed
{
#the translocation area needs to be removed both from the p-score and from the peaks file
translocated_flag <- 1
}
break
}
}
#gets the names of the experiments to later be entered into the 'ChIPseqVScontacts_correlation_plots' file
if(z==1)
{
ps_Experiments <- ps.file
}
else
{
ps_Experiments <- paste(ps_Experiments," ; ",ps.file,sep="")
}
#importing the p-score data from the file
ps <- read.table(paste("~/Analyze4C/pScores/",ps_folder,"/",ps.file,sep=""))
#finding the specific p-scores files information in Experiments_4C
sp1 <- unlist(strsplit(ps.file,"[.]"))
sp2 <- unlist(strsplit(sp1[1],"_"))
#getting the size of the window per RE sites
sp3 <- sp2[as.integer(grep("RE$",sp2))]
RE_wind <- as.numeric(strsplit(sp3,"RE"))
if(choice3 == 2)
{
#getting the size of the window per bp sites
sp4 <- sp2[as.integer(grep("bp$",sp2))]
bp_wind <- as.numeric(strsplit(sp4,"bp"))
}
#finding the specific p-scores files information in Experiments_4C
out <- findIn.Experiments_4C(sp2[1],sp2[2],sp2[3],Experiments_4C)
#getting the cis chromosome number
cis <- as.numeric(out[2])
}
else
{
#choosing an peaks file
cat("\npeaks file number ",z,":\n",sep="")
ls_files_peaks <- system("ls ~/Analyze4C/ChIPseq/peaks | grep 'tagsANDpscores.bed$'",intern=TRUE)
repeat
{
file.name_peaks <- readline(prompt=cat("\nThese are the ChIPseq data (tags or p-scores) files available\nwhich would you like to use?\n",ls_files_peaks,"",sep="\n"))
ind_files_peaks <- pmatch(file.name_peaks,ls_files_peaks,nomatch=0)
if(ind_files_peaks == 0)
{
cat("no such file exists.\nplease try again.\n\n")
}
else
{
#importing the data from the file
peaks <- read.table(paste("~/Analyze4C/ChIPseq/peaks/",file.name_peaks,sep=""))
break
}
#here i might need to edit and sort the data in peaks
#peaks <- peaks[peaks[,1]!="Pt" & peaks[,1]!="Mt",]
#peaks <- peaks[order(peaks[,1],peaks[,2]),]
#names(peaks) <- c("chr","first","last","peaks")
}
#gets the names of the experiments to later be entered into the 'ChIPseqVScontacts_correlation_plots' file
if(z==1)
{
peaks_Experiments <- file.name_peaks
}
else
{
peaks_Experiments <- paste(peaks_Experiments," ; ",file.name_peaks,sep="")
}
}
if(choice4 == 2 & wind_flag == 0)
{
genome_sizes_temp <- as.data.frame(cbind(gsub("^chr","",genome_sizes[,1]),genome_sizes[,2])) #removing 'chr' from the chromosome numbers, if there are any with 'chr'
write.table(genome_sizes_temp,"~/Analyze4C/temp/genome_sizes_temp.txt",row.names=FALSE,col.names=FALSE,quote=FALSE,sep="\t") #creating a temp file of the genome sizes
winSize <- as.integer(readline(prompt=cat("\nplease enter a window size (bp) that the whole genome will divided into:\n\n")))
system(paste("bedtools makewindows -g ~/Analyze4C/temp/genome_sizes_temp.txt -w",winSize,"> ~/Analyze4C/temp/windows.bed"))
}
#getting the cutoff for the p-score data
ans8 <- as.integer(readline(prompt=cat("\nenter the number of what type of cutoff method you would like to do for the p-score:\n\n1) all chromosomes with the same cutoff (in p-score)\n2) all chromosomes with the same cutoff percentage\n\n")))
if(ans8 == 1)
{
CO1 <- as.numeric(readline(prompt=cat("\nplease enter a p-score cutoff. Anything above the cutoff will be considered a positive RE site:\n\n")))
CO_ps_type <- paste(CO_ps_type,";","fixed")
CO_ps <- paste(CO_ps,";",CO1)
ans4 <- as.integer(readline(prompt=cat("\nshould the cutoff be applied to:\n1) the whole genome\n2) all trans\n3) only cis\n\n")))
}
else
{
ps_COper <- as.numeric(readline(prompt=cat("\nenter a cutoff percentage (between 0 and 1) that will be applied to the p-score\nonly the p-scores that are above or equal to the cutoff will be intersected with the peaks data.\nif you don't want any cutoff, enter 0:\n")))
CO_ps_type <- paste(CO_ps_type,";","top percent per chromosome")
CO_ps <- paste(CO_ps,";",ps_COper)
ans4 <- as.integer(readline(prompt=cat("\nshould the cutoff percentage be of and applied to:\n1) the whole genome\n2) all trans\n3) each chromosome separately\n\n")))
#getting the actual cutoffs (we do this here since we need to consider the translocated areas before removing them, if there are any, in order to calculate the cutoffs)
if(ans4 == 1) #whole genome
{
CO1 <- quantile(ps[,3],ps_COper)
CO_ps_chroms <- paste(CO_ps_chroms,";","all")
}
else if(ans4 == 2) #trans
{
CO1 <- quantile(ps[ps[,1]!=cis,3],ps_COper)
CO_ps_chroms <- paste(CO_ps_chroms,";","trans")
}
else if(ans4 == 3) #each chromosome separately
{
seperate_CO1 <- rep(0,numOFchroms)
for(k in 1:numOFchroms)
{
seperate_CO1[k] <- quantile(ps[ps[,1]==k,3],ps_COper)
}
CO_ps_chroms <- paste(CO_ps_chroms,";","each chromosome separately")
}
}
#choosing the value from which the cutoffs are created and the value which is going to be used for each peak above and equal to the cutoff
cat("\napplying cutoffs to the peaks data:\n")
ans11 <- as.integer(readline(prompt=cat("\nchoose an option:\n1) apply the cutoff on the p-scores, get the peak values from p-scores\n2) apply the cutoff on the p-scores, get the peak values from number of tags\n3) apply the cutoff on the number of tags, get the peak values from p-scores\n4) apply the cutoff on the number tags, get the peak values from number of tags\n\n")))
if(ans11==1)
{
CO_col <- 5
value_col <- 5
CO_source <- "p-scores"
value_source <- "p-scores"
}
else if(ans11==2)
{
CO_col <- 5
value_col <- 4
CO_source <- "p-scores"
value_source <- "tags"
}
else if(ans11==3)
{
CO_col <- 4
value_col <- 5
CO_source <- "tags"
value_source <- "p-scores"
}
else if(ans11==4)
{
CO_col <- 4
value_col <- 4
CO_source <- "tags"
value_source <- "tags"
}
#applying a cutoff to the peaks data
if(CO_col == 4)
{
ans9 <- as.integer(readline(prompt=cat("\nenter the number of what type of cutoff method you would like to do for the peaks:\n\n1) all chromosomes with the same cutoff (in number of tags)\n2) all chromosomes with the same cutoff percentage\n\n")))
}
else if(CO_col == 5)
{
ans9 <- as.integer(readline(prompt=cat("\nenter the number of what type of cutoff method you would like to do for the peaks:\n\n1) all chromosomes with the same cutoff (in p-scores)\n2) all chromosomes with the same cutoff percentage\n\n")))
}
if(ans9 == 1)
{
if(CO_col == 4)
{
CO2 <- as.numeric(readline(prompt=cat("\nplease enter a number of tags cutoff. Anything above the cutoff will be considered:\n\n")))
}
else if(CO_col == 5)
{
CO2 <- as.numeric(readline(prompt=cat("\nplease enter a p-score cutoff. Anything above the cutoff will be considered:\n\n")))
}
CO_peaks_type <- paste(CO_peaks_type,";","fixed")
CO_peaks <- paste(CO_peaks,";",CO2)
ans5 <- as.integer(readline(prompt=cat("\nshould the cutoff be applied to:\n1) the whole genome\n2) all trans\n3) cis\n\n")))
if(ans5 == 1) #whole genome
{
peaks_afterCO <- peaks[peaks[,CO_col]>=CO2,c(1:3,value_col)]
CO_peaks_chroms <- paste(CO_peaks_chroms,";","all")
}
else if(ans5 == 2) #trans
{
peaks_afterCO <- c()
for(k in 1:numOFchroms)
{
if(k!=cis)
{
peaks_temp <- peaks[peaks[,1]==k & peaks[,CO_col]>=CO2,c(1:3,value_col)]
peaks_afterCO <- rbind(peaks_afterCO,peaks_temp)
}
else
{
peaks_temp <- peaks[peaks[,1]==k,c(1:3,value_col)]
peaks_afterCO <- rbind(peaks_afterCO,peaks_temp)
}
}
CO_peaks_chroms <- paste(CO_peaks_chroms,";","trans")
}
else if(ans5 == 3) #cis
{
peaks_afterCO <- c()
for(k in 1:numOFchroms)
{
if(k==cis)
{
peaks_temp <- peaks[peaks[,1]==k & peaks[,CO_col]>=CO2,c(1:3,value_col)]
peaks_afterCO <- rbind(peaks_afterCO,peaks_temp)
}
else
{
peaks_temp <- peaks[peaks[,1]==k,c(1:3,value_col)]
peaks_afterCO <- rbind(peaks_afterCO,peaks_temp)
}
}
CO_peaks_chroms <- paste(CO_peaks_chroms,";","cis")
}
}
else
{
if(CO_col == 4)
{
peaks_COper <- as.numeric(readline(prompt=cat("\nenter a cutoff percentage (between 0 and 1) that will be applied to the number of tags of the peaks\nonly the peaks that are above or equal to the cutoff will be intersected to the reads.\nif you don't want any cutoff, enter 0:\n")))
}
else if(CO_col == 5)
{
peaks_COper <- as.numeric(readline(prompt=cat("\nenter a cutoff percentage (between 0 and 1) that will be applied to the p-scores of the peaks\nonly the peaks that are above or equal to the cutoff will be intersected to the reads.\nif you don't want any cutoff, enter 0:\n")))
}
CO_peaks_type <- paste(CO_peaks_type,";","top percent per chromosome")
CO_peaks <- paste(CO_peaks,";",peaks_COper)
ans5 <- as.integer(readline(prompt=cat("\nshould the cutoff percentage be of:\n1) the whole genome\n2) all trans\n3) each chromosome separately\n\n")))
if(ans5 == 1) #whole genome
{
CO2 <- quantile(peaks[,CO_col],peaks_COper)
peaks_afterCO <- peaks[peaks[,CO_col]>=CO2,c(1:3,value_col)]
CO_peaks_chroms <- paste(CO_peaks_chroms,";","all")
}
else if(ans5 == 2) #trans
{
peaks_afterCO <- c()
CO2 <- quantile(peaks[peaks[,1]!=cis,CO_col],peaks_COper)
for(k in 1:numOFchroms)
{
if(k!=cis)
{
peaks_temp <- peaks[peaks[,1]==k & peaks[,CO_col]>=CO2,c(1:3,value_col)]
peaks_afterCO <- rbind(peaks_afterCO,peaks_temp)
}
else
{
peaks_temp <- peaks[peaks[,1]==k,c(1:3,value_col)]
peaks_afterCO <- rbind(peaks_afterCO,peaks_temp)
}
}
CO_peaks_chroms <- paste(CO_peaks_chroms,";","trans")
}
else if(ans5 == 3) #each chromosome separately
{
peaks_afterCO <- c()
for(k in 1:numOFchroms)
{
CO2 <- quantile(peaks[peaks[,1]==k,CO_col],peaks_COper)
peaks_temp <- peaks[peaks[,1]==k & peaks[,CO_col]>=CO2,c(1:3,value_col)]
peaks_afterCO <- rbind(peaks_afterCO,peaks_temp)
}
CO_peaks_chroms <- paste(CO_peaks_chroms,";","each chromosome separately")
}
}
#adjusting the row numbers
rownames(peaks_afterCO) <- seq(length=nrow(peaks_afterCO))
#removing the translocated areas from 'rawData' and 'peaks_afterCO' if they were switched
if(translocated_flag == 1)
{
cat("\nremoving translocated areas:\nthe purpose of this is in order to be able and intersect the p-scores with peaks bands\n")
cat("the translocated areas indices mess up the order and will cause an error\n")
#getting the line numbers that need to be removed
for(d in 1:length(rearranged_rawData$Experiment))
{
sTemp <- unlist(strsplit(rearranged_rawData$Experiment[d],"_"))
dtTemp <- unlist(strsplit(rearranged_rawData$Date_and_Time[d]," "))
tTemp <- dtTemp[2]
dtTemp[2] <- paste(unlist(strsplit(tTemp,":")),collapse="")
if(all(sTemp %in% sp2) & !is.na(rearranged_rawData$Added_Sections_Lines[d]) & all(dtTemp %in% sp2))
{
lin_nums1 <- rearranged_rawData$Added_Sections_Lines[d]
lin_nums2 <- as.integer(unlist(strsplit(lin_nums1,", ")))
#getting the raw data file that the p-scores file was created from:
if(!identical(grep("covRemoved",sp1),integer(0)))
{
ls_raws <- system("ls ~/Analyze4C/rawData/coverage_removed",intern=TRUE)
}
else
{
ls_raws <- system("ls ~/Analyze4C/rawData/rearranged",intern=TRUE)
}
dAndtspl <- paste(dtTemp[1],"_",dtTemp[2],sep="")
raws_inds <- grep(dAndtspl,ls_raws)
for(o in raws_inds)
{
if(identical(grep("removed",ls_raws[o]),integer(0)))
{
if(!identical(grep("covRemoved",sp1),integer(0)))
{
#getting the raw data
raw_dat <- read.table(paste("~/Analyze4C/rawData/coverage_removed/",ls_raws[o],sep=""))
break
}
else
{
#getting the raw data
raw_dat <- read.table(paste("~/Analyze4C/rawData/rearranged/",ls_raws[o],sep=""))
break
}
}
}
secs_rmv <- "yes"
break
}
}
#the following actions will work under these specific conditions:
#1) the rearranged raw data file contains all the original chromosomes.
#2) the line numbers in 'rearranged_rawData' file are divided so that each pair are of one specific chromosome (meaning that 'first' and 'second' are not of two different chromosomes)
#removing the sections from the peaks and ps data
all_chroms <- seq(numOFchroms)
pointer <- 1
peaks_fin <- c()
#this is for removing sections from the raw data
rawData_removed <- cbind(raw_dat,rep(1,nrow(raw_dat)))
for(r in seq(1,length(lin_nums2),by=2))
{
first <- lin_nums2[r]
second <- lin_nums2[r+1]
ch <- raw_dat[first,1]
if(ch != raw_dat[second,1])
{
cat("\nerror:there is something wrong with the line numbers entered into the file 'rearranged_rawData'.\nthe first and last lines must be of the same chromosome\n\n")
return()
}
#this is for removing sections from the raw data
rawData_removed[first:second,ncol(rawData_removed)] <- 0
#taking all the chromosomes that aren't getting the removal treatment, and adding them to the peaks and ps data
while(all_chroms[pointer]!=ch)
{
peaks_t <- peaks_afterCO[peaks_afterCO[,1]==all_chroms[pointer],]
peaks_fin <- rbind(peaks_fin,peaks_t)
pointer <- pointer + 1
}
peaks_t <- c() #just in case no conditions are met, this way we won't get the previous peaks_t (this shouldn't happen, it's just a precaution)
#getting the indices to remove
if(first==1)#if first is the 1st line in the file
{
if(second!=nrow(raw_dat)) #if second isn't the last of the file
{
if(raw_dat[second,1]!=raw_dat[(second+1),1]) #if it covers all the chromosome (if second is the last in chromosome) (1)
{
peaks_t <- c()
}
else if(raw_dat[second,1]==raw_dat[(second+1),1]) #if second isn't the last of the chromosome (2)
{
#remove anything below the index at second
peaks_t <- peaks_afterCO[(peaks_afterCO[,1]==ch)&(peaks_afterCO[,2]>=raw_dat[(second+1),2]),]
}
}
else #the line numbers cover all the data, so none of it is used (9)
{
peaks_t <- c()
}
}
else if(raw_dat[first,1]!=raw_dat[(first-1),1]) #if it is the first line in chromosome
{
if(second!=nrow(raw_dat)) #if second isn't last line
{
if(raw_dat[second,1]!=raw_dat[(second+1),1]) #if second is the last in chromosome (4)
{
peaks_t <- c()
}
else if(raw_dat[second,1]==raw_dat[(second+1),1]) #if second isn't last in chromosome (3)
{
#remove anything below the index at second
peaks_t <- peaks_afterCO[(peaks_afterCO[,1]==ch)&(peaks_afterCO[,2]>=raw_dat[(second+1),2]),]
}
}
else #if second is last line of the file (7)
{
peaks_t <- c()
}
}
else if(raw_dat[first,1]==raw_dat[(first-1),1])#if first is in the middle of a chromosome
{
if(second!=nrow(raw_dat)) #if second isn't last line
{
if(raw_dat[second,1]!=raw_dat[(second+1),1]) #if second is the last in chromosome (6)
{
#remove anything above the index at first
peaks_t <- peaks_afterCO[(peaks_afterCO[,1]==ch)&(peaks_afterCO[,2]<=raw_dat[(first-1),2]),]
}
else if(raw_dat[second,1]==raw_dat[(second+1),1]) #if second isn't last in chromosome (5)
{
#remove anything above the index at first and below the index at second
peaks_t <- peaks_afterCO[(peaks_afterCO[,1]==ch)&((peaks_afterCO[,2]<=raw_dat[(first-1),2])|(peaks_afterCO[,2]>=raw_dat[(second+1),2])),]
}
}
else #if second is last line of the file (8)
{
#remove anything above the index at first
peaks_t <- peaks_afterCO[(peaks_afterCO[,1]==ch)&(peaks_afterCO[,2]<=raw_dat[(first-1),2]),]
}
}
peaks_fin <- rbind(peaks_fin,peaks_t)
pointer <- pointer + 1
}
#adding the last chromosomes data (if there are any) that didn't get the section removal treatment
while(pointer<=numOFchroms)
{
peaks_t <- peaks_afterCO[peaks_afterCO[,1]==all_chroms[pointer],]
peaks_fin <- rbind(peaks_fin,peaks_t)
pointer <- pointer + 1
}
peaks_afterCO <- peaks_fin
rownames(peaks_afterCO) <- seq(length=nrow(peaks_afterCO))
#this is for removing sections from the raw data
rawData_removed <- rawData_removed[rawData_removed[,ncol(rawData_removed)]==1,-ncol(rawData_removed)]
rownames(rawData_removed) <- seq(length=nrow(rawData_removed))
}
if(z == 1) #maybe i could remove this condition and have that every time you could enter which ever chrom you want, even different ones
{
choice2 <- as.integer(readline(prompt=cat("\nwhat chromosomes should the intersection be applied to? (enter the option number)\n1) all chromosomes\n2) trans\n3) cis\n4) specific chromosome\n\n")))
}
#getting the p-score and peaks data by chromosome according to the request of the user, and also the chromosomes chosen
if(choice2 == 1) #all chromosomes
{
ps_dat <- ps
peaks_dat <- peaks_afterCO
chroms <- c(1:numOFchroms)
int_chroms <- "all" #this will be entered later in 'ChIPseqVScontacts_correlation_plots.txt'
}
else if(choice2 == 2) #trans
{
ps_dat <- ps[ps[,1] != cis,]
peaks_dat <- peaks_afterCO[peaks_afterCO[,1] != cis,]
chroms <- Filter(function(x) x!=cis,1:numOFchroms) #this here excludes cis chromosome number in the sequence, another way would be - (1:numOFchroms)[-cis]
int_chroms <- "trans" #this will be entered later in 'ChIPseqVScontacts_correlation_plots.txt'
}
else if(choice2 == 3) #cis
{
ps_dat <- ps[ps[,1] == cis,]
peaks_dat <- peaks_afterCO[peaks_afterCO[,1] == cis,]
chroms <- cis
int_chroms <- "cis" #this will be entered later in 'ChIPseqVScontacts_correlation_plots.txt'
}
else if(choice2 == 4) #by chromosome
{
chr_chosen <- as.integer(readline(prompt=cat("\nchoose a chromosome you would like to use:\n\n")))
ps_dat <- ps[ps[,1] == chr_chosen,]
peaks_dat <- peaks_afterCO[peaks_afterCO[,1] == chr_chosen,]
chroms <- chr_chosen
int_chroms <- paste("chr",chr_chosen,sep="") #this will be entered later in 'ChIPseqVScontacts_correlation_plots.txt'
}
#add an index to each p-score file
ps_bed <- c()
for(u in chroms)
{
if(nrow(ps_dat[ps_dat[,1]==u,]) != 0)
{
ps_chrom_temp <- ps_dat[ps_dat[,1]==u,]
inds_temp <- ps_chrom_temp[,2]+1
if(tail(inds_temp,1)>genome_sizes[u,2]) #if the addition of the extra index is larger than the size of the chromosome, we will remove that p-score line
{
inds_temp <- head(inds_temp,-1)
ps_chrom_temp <- head(ps_chrom_temp,-1)
}
ps_chrom_bed <- cbind(ps_chrom_temp[,1],ps_chrom_temp[,2],inds_temp,ps_chrom_temp[,3])
ps_bed <- rbind(ps_bed,ps_chrom_bed)
}
}
ps_bed <- data.frame(ps_bed)
write.table(ps_bed,"~/Analyze4C/temp/ps_bed.bed",row.names=FALSE,col.names=FALSE,quote=FALSE,sep="\t") #creating a temp file of the p-score data chosen
write.table(peaks_dat,"~/Analyze4C/temp/peaks.bed",row.names=FALSE,col.names=FALSE,quote=FALSE,sep="\t") #creating a temp file of the peaks data chosen
##!
#getting the raw data that corresponds to the p-score file, in order to get the corresponding reads for each RE site
#cat("\nthis is the current p-score file that is being tested:\n",ps.file,"\n")
#repeat
#{
# choice <- as.integer(readline(prompt=cat("\nenter the number of the folder from which contains the raw data that the p-score file was created from (look at the name of the p-score file for clues):\n\n1) original\n2) rearranged\n3) coverage removed\n\n")))
# if(choice == 2)
# {
# ls_files <- system("ls ~/Analyze4C/rawData/rearranged",intern=TRUE)
# }
# else if(choice == 3)
# {
# ls_files <- system("ls ~/Analyze4C/rawData/coverage_removed",intern=TRUE)
# }
# else
# {
# ls_files <- system("ls ~/Analyze4C/rawData/original",intern=TRUE)
# }
#
# if(length(ls_files) != 0)
# {
# file.name <- readline(prompt=cat("\nchoose from these raw data files (make sure the name is similar to the p-score files):\n",ls_files,"\n",sep="\n"))
# ind_files <- pmatch(file.name,ls_files,nomatch=0)
# if(ind_files == 0)
# {
# cat("no such file exists.\nplease try again.\n\n")
# }
# else
# {
# break
# }
# }
# else
# {
# cat("\nthis folder is empty. please choose a different folder from which you would like to use a file\n\n")
# }
#}
#importing the data from the file
#if(choice == 2)
#{
# rawData <- read.table(paste("~/Analyze4C/rawData/rearranged/",file.name,sep=""))
#}
#else if(choice == 3)
#{
# rawData <- read.table(paste("~/Analyze4C/rawData/coverage_removed/",file.name,sep=""))
#}
#else
#{
# rawData <- read.table(paste("~/Analyze4C/rawData/original/",file.name,sep=""))
#}
##!
#add an index to the rawData_removed file
rawData_bed <- c()
for(u in chroms)
{
if(nrow(rawData_removed[rawData_removed[,1]==u,]) != 0)
{
rawData_chrom_temp <- rawData_removed[rawData_removed[,1]==u,]
rawData_inds_temp <- rawData_chrom_temp[,2]+1
if(tail(rawData_inds_temp,1)>genome_sizes[u,2]) #if the addition of the extra index is larger than the size of the chromosome, we will remove that p-score line
{
rawData_inds_temp <- head(rawData_inds_temp,-1)
rawData_chrom_temp <- head(rawData_chrom_temp,-1)
}
rawData_chrom_bed <- cbind(rawData_chrom_temp[,1],rawData_chrom_temp[,2],rawData_inds_temp,rawData_chrom_temp[,3])
rawData_bed <- rbind(rawData_bed,rawData_chrom_bed)
}
}
rawData_bed <- data.frame(rawData_bed)
write.table(rawData_bed,"~/Analyze4C/temp/rawData_bed.bed",row.names=FALSE,col.names=FALSE,quote=FALSE,sep="\t") #creating a temp file of the rawData data chosen
#removing the sections from the data:
#intersecting the bed files so that we get the indices that aren't removed (since the way we removed was with line numbers of the files, we can't just remove them from the p-score data since we removed some lines from the cis when calculating the p-scores)
system("bedtools intersect -a ~/Analyze4C/temp/ps_bed.bed -b ~/Analyze4C/temp/rawData_bed.bed > ~/Analyze4C/temp/out.bed")
system("bedtools intersect -a ~/Analyze4C/temp/ps_bed.bed -b ~/Analyze4C/temp/rawData_bed.bed -c > ~/Analyze4C/temp/out2.bed") #this output will have all the same lines as psBed.bed but with an indication if the indices intersected or not (0 not intersected, 1 intersected)
#there was a problem with intersected repeats, where if the same indices appeared in the translocated area that needs to be removed and an area that doesn'table
#this will cause repeats in the intersections and some indices to be out of order
#i couldn't figure out an automatic way to fix this
#for now the user will fix it manually
cat("\nview the 'out.bed' file for any repeats of intesections.\nthis happens if there is an index that appears in the translocated area and in areas that aren't.\n")
cat("bedtools doesn't know how to distinguish between 'good' intersections like these and 'bad' ones.\nthese repeats should be at the beginning of the chromosome (if the translocation occured on the left telomere of the chromosome.\n")
cat("you could confirm that these are the lines that need to be erased by looking at the 'out2.bed' which shows that there are some indices that appear early on and in translocated areas that have been intersected.\n")
cat("most chances are that you will see a '1' at the line of some indices with other indices that have a '0' surrounding them\n")
lns <- c()
fl <- 1
cat("\n\nentering the lines that need to be removed:\n")
while(fl == 1)
{
ans <- readline(prompt=cat("\nare there any (more) lines that need to be erased?\ny/n\n\n"))
if(ans == "y")
{
new_ln <- as.integer(readline(prompt=cat("\nenter the line number that needs to be removed, as it appears in 'out.bed':\n")))
lns <- c(lns,new_ln)
}
else
{
fl <- 0
}
}
ps_bed2 <- read.table("~/Analyze4C/temp/out.bed")
if(!is.null(lns))
{
ps_bed2 <- ps_bed2[-1*lns,] #removing the specific lines
}
#applying the cutoff to the p-scores
if(ans4 == 1) #whole genome
{
ps_bed_afterCO <- ps_bed2[ps_bed2[,4]>=CO1,]
}
else if(ans4 == 2) #trans
{
ps_bed_afterCO <- c()
for(k in 1:numOFchroms)
{
if(k!=cis)
{
ps_temp <- ps_bed2[ps_bed2[,1]==k & ps_bed2[,4]>=CO1,]
ps_bed_afterCO <- rbind(ps_bed_afterCO,ps_temp)
}
else
{
ps_temp <- ps_bed2[ps_bed2[,1]==k,]
ps_bed_afterCO <- rbind(ps_bed_afterCO,ps_temp)
}
}
}
else if(ans4 == 3) #cis or each chromosome separately
{
ps_bed_afterCO <- c()
for(k in 1:numOFchroms)
{
if(ans8 == 1) #only cis
{
if(k==cis)
{
ps_temp <- ps_bed2[ps_bed2[,1]==k & ps_bed2[,4]>=CO1,]
ps_bed_afterCO <- rbind(ps_bed_afterCO,ps_temp)
}
else
{
ps_temp <- ps_bed2[ps_bed2[,1]==k,]
ps_bed_afterCO <- rbind(ps_bed_afterCO,ps_temp)
}
}
else #each chromosome separately
{
ps_temp <- ps_bed2[ps_bed2[,1]==k & ps_bed2[,4]>=seperate_CO1[k],]
ps_bed_afterCO <- rbind(ps_bed_afterCO,ps_temp)
}
}
}
#adjusting the row numbers
rownames(ps_bed_afterCO) <- seq(length=nrow(ps_bed_afterCO))
write.table(ps_bed_afterCO,"~/Analyze4C/temp/ps_bed_afterCO.bed",row.names=FALSE,col.names=FALSE,quote=FALSE,sep="\t")
#ps_dat2 <- cbind(ps_bed2[,1],ps_bed2[,2],ps_bed2[,4])
system("bedtools intersect -a ~/Analyze4C/temp/ps_bed_afterCO.bed -b ~/Analyze4C/temp/rawData_bed.bed -loj > ~/Analyze4C/temp/reads.bed")
reads <- read.table("~/Analyze4C/temp/reads.bed",header=FALSE,stringsAsFactors=FALSE,sep="\t")
reads_bed <- as.data.frame(cbind(ps_bed_afterCO[,1],ps_bed_afterCO[,2],ps_bed_afterCO[,3],reads[,8]))
write.table(reads_bed,"~/Analyze4C/temp/reads_bed.bed",row.names=FALSE,col.names=FALSE,quote=FALSE,sep="\t")
#intersecting
if(choice4 == 1) #if by individual
{
system("bedtools intersect -a ~/Analyze4C/temp/ps_bed_afterCO.bed -b ~/Analyze4C/temp/peaks.bed -wo > ~/Analyze4C/temp/intersected.bed")
#if there were no intersections (i think this is the best way to deal with a case like this)
if(file.info("~/Analyze4C/temp/intersected.bed")$size == 0)
{
cat("\nthere were no intersections between the p-scores and peaks data\n\n")
return()
}
#if there were intersections:
intersected <- read.table("~/Analyze4C/temp/intersected.bed")
intersected_pairs <- as.data.frame(cbind(intersected[,4],intersected[,8]))
system("bedtools intersect -a ~/Analyze4C/temp/reads_bed.bed -b ~/Analyze4C/temp/peaks.bed -wo > ~/Analyze4C/temp/intersected_reads.bed")
#if there were no intersections (i think this is the best way to deal with a case like this)
if(file.info("~/Analyze4C/temp/intersected_reads.bed")$size == 0)
{
cat("\nthere were no intersections between the reads and peaks data\n\n")
return()
}
#if there were intersections:
intersected_reads <- read.table("~/Analyze4C/temp/intersected_reads.bed")
intersected_reads_pairs <- as.data.frame(cbind(intersected_reads[,4],intersected_reads[,8]))
}
else #if by window
{
system("bedtools map -a ~/Analyze4C/temp/windows.bed -b ~/Analyze4C/temp/ps_bed_afterCO.bed -c 4 > ~/Analyze4C/temp/ps_windows.bed")
ps_windows <- read.table("~/Analyze4C/temp/ps_windows.bed",header=FALSE,stringsAsFactors=FALSE,sep="\t")
ps_windows[ps_windows[,4]==".",4] <- 0
#write.table(ps_windows,"~/Analyze4C/temp/ps_windows.txt",row.names=FALSE,col.names=FALSE,quote=FALSE,sep="\t")
system("bedtools map -a ~/Analyze4C/temp/windows.bed -b ~/Analyze4C/temp/peaks.bed -c 4 > ~/Analyze4C/temp/peaks_windows.bed")
peaks_windows <- read.table("~/Analyze4C/temp/peaks_windows.bed",header=FALSE,stringsAsFactors=FALSE,sep="\t")
peaks_windows[peaks_windows[,4]==".",4] <- 0
#write.table(peaks_windows,"~/Analyze4C/temp/peaks_windows.txt",row.names=FALSE,col.names=FALSE,quote=FALSE,sep="\t")
intersected_pairs <- as.data.frame(cbind(as.numeric(ps_windows[,4]),as.numeric(peaks_windows[,4])))
system("bedtools map -a ~/Analyze4C/temp/windows.bed -b ~/Analyze4C/temp/reads_bed.bed -c 4 > ~/Analyze4C/temp/reads_windows.bed")
reads_windows <- read.table("~/Analyze4C/temp/reads_windows.bed",header=FALSE,stringsAsFactors=FALSE,sep="\t")
reads_windows[reads_windows[,4]==".",4] <- 0
#write.table(reads_windows,"~/Analyze4C/temp/reads_windows.txt",row.names=FALSE,col.names=FALSE,quote=FALSE,sep="\t")
intersected_reads_pairs <- as.data.frame(cbind(as.numeric(reads_windows[,4]),as.numeric(peaks_windows[,4])))
}
#adding the name given by user to the data, and creating a column for it
if(choice1 == 1)
{
ps_datName <- readline(prompt=cat("\nhow would you like this p-score data to be refered to in the plot? (no spaces)\n\n"))
ps_datName_all <- c(ps_datName_all,ps_datName)
intersected_pairs <- cbind(intersected_pairs,rep(ps_datName,nrow(intersected_pairs)))
intersected_reads_pairs <- cbind(intersected_reads_pairs,rep(ps_datName,nrow(intersected_reads_pairs)))
}
else
{
peaks_datName <- readline(prompt=cat("\nhow would you like this peaks data to be refered to in the plot? (no spaces)\n\n"))
peaks_datName_all <- c(peaks_datName_all,peaks_datName)
intersected_pairs <- cbind(intersected_pairs,rep(peaks_datName,nrow(intersected_pairs)))
intersected_reads_pairs <- cbind(intersected_reads_pairs,rep(peaks_datName,nrow(intersected_reads_pairs)))
}
colnames(intersected_pairs) <- c("pScore","peaks","names")
#intersected_pairs <- as.data.frame(intersected_pairs)
intersected_pairs_all <- rbind(intersected_pairs_all,intersected_pairs)
colnames(intersected_reads_pairs) <- c("reads","peaks","names")
#intersected_reads_pairs <- as.data.frame(intersected_reads_pairs)
intersected_reads_pairs_all <- rbind(intersected_reads_pairs_all,intersected_reads_pairs)
#spearman
cat("\nthe spearman correlation coefficient for the p-score vs. peaks:\n")
print(cor.test(intersected_pairs$pScore,intersected_pairs$peaks,method="spearman",exact=FALSE))
#print(cor.test(as.numeric(intersected_pairs$pScore),as.numeric(intersected_pairs$peaks),method="spearman",exact=FALSE))
#pearson
cat("\nthe pearson correlation coefficient for the p-score vs. peaks:\n")
print(cor.test(intersected_pairs$pScore,intersected_pairs$peaks))
#print(cor.test(as.numeric(intersected_pairs$pScore),as.numeric(intersected_pairs$peaks)))
#spearman
cat("\nthe spearman correlation coefficient for the reads vs. peaks:\n")
print(cor.test(intersected_reads_pairs_all$reads,intersected_reads_pairs_all$peaks,method="spearman",exact=FALSE))
#print(cor.test(as.numeric(intersected_reads_pairs_all$reads),as.numeric(intersected_reads_pairs_all$peaks),method="spearman",exact=FALSE))
#pearson
cat("\nthe pearson correlation coefficient for the reads vs. peaks:\n")
print(cor.test(intersected_reads_pairs_all$reads,intersected_reads_pairs_all$peaks))
#print(cor.test(as.numeric(intersected_reads_pairs_all$reads),as.numeric(intersected_reads_pairs_all$peaks)))
if(choice1 == 1)
{
#asking if to choose another p-score file to correlate
ans6 <- readline(prompt=cat("\nwould you like to choose another p-score file to intersect with the peaks data?\ny/n\n\n"))
if(ans6 == "y")
{
z <- z + 1
system("rm ~/Analyze4C/temp/ps_bed.bed")
system("rm ~/Analyze4C/temp/ps_bed_afterCO.bed")
system("rm ~/Analyze4C/temp/reads.bed")
system("rm ~/Analyze4C/temp/reads_bed.bed")
system("rm ~/Analyze4C/temp/rawData_bed.bed")
system("rm ~/Analyze4C/temp/out.bed")
system("rm ~/Analyze4C/temp/out2.bed")
if(choice4 == 1)
{
system("rm ~/Analyze4C/temp/intersected.bed")
system("rm ~/Analyze4C/temp/intersected_reads.bed")
}
else
{
winds <- c(winds," ; ",winSize)
system("rm ~/Analyze4C/temp/ps_windows.bed")
system("rm ~/Analyze4C/temp/reads_windows.bed")
system("rm ~/Analyze4C/temp/peaks_windows.bed")
choice5 <- readline(prompt=cat("\nwould you like to choose a different window size?\ny/n\n\n"))
if(choice5 == "y")
{
system("rm ~/Analyze4C/temp/windows.bed")
wind_flag <- 0
}
else
{
wind_flag <- 1
}
}
}
else
{
files_flag <- 1
num_of_pScore_files <- z
#erasing all the files from temp
system("rm ~/Analyze4C/temp/*")
}
}
else
{
#asking if to choose another peaks file to correlate
ans6 <- readline(prompt=cat("\nwould you like to choose another peaks file to intersect with the peaks data?\ny/n\n\n"))
if(ans6 == "y")
{
z <- z + 1
system("rm ~/Analyze4C/temp/peaks.bed")
system("rm ~/Analyze4C/temp/reads.bed")
system("rm ~/Analyze4C/temp/reads_bed.bed")
system("rm ~/Analyze4C/temp/rawData_bed.bed")
if(choice4 == 1)
{
system("rm ~/Analyze4C/temp/intersected.bed")
system("rm ~/Analyze4C/temp/intersected_reads.bed")
}
else
{
winds <- c(winds," ; ",winSize)
system("rm ~/Analyze4C/temp/peaks_windows.bed")
system("rm ~/Analyze4C/temp/ps_windows.bed")
system("rm ~/Analyze4C/temp/reads_windows.bed")
choice5 <- readline(prompt=cat("\nwould you like to choose a different window size?\ny/n\n\n"))
if(choice5 == "y")
{
system("rm ~/Analyze4C/temp/windows.bed")
wind_flag <- 0
}
else
{
wind_flag <- 1
}
}
}
else
{
files_flag <- 1
num_of_peaks_files <- z
#erasing all the files from temp
system("rm ~/Analyze4C/temp/*")
}
}
}
#getting the time and date
DandT1 <- toString(Sys.time())
DandT2 <- gsub(" ","_",DandT1)
DandT2 <- gsub(":","",DandT2)
#intersected_pairs_all <- as.data.frame(intersected_pairs_all)
cat("\ncreating the p-score VS. peaks plot\nplease wait...\n\n")
#print(ggplot2::ggplot(intersected_pairs_all, ggplot2::aes(x=pScore,y=peaks,color=names)) + ggplot2::scale_x_log10() + ggplot2::scale_y_log10() + ggplot2::geom_point() + ggplot2::geom_smooth(method=lm,se=FALSE,fullrange=TRUE))
#print(ggplot2::ggplot(intersected_pairs_all, ggplot2::aes(x=pScore,y=peaks,color=names)) + ggplot2::geom_point() + ggplot2::geom_smooth(method=lm,se=FALSE,fullrange=TRUE) + ggplot2::scale_y_continuous(limits = c(0,5000)) + ggplot2::ggtitle("contact VS. expression - correlation") + ggplot2::theme(plot.title = ggplot2::element_text(size=15,hjust = 0.5)))
print(ggplot2::ggplot(intersected_pairs_all, ggplot2::aes(x=pScore,y=peaks,color=names)) + ggplot2::geom_point() + ggplot2::geom_smooth(method=lm,se=FALSE,fullrange=TRUE) + ggplot2::ggtitle("contact VS. ChIPseq - correlation") + ggplot2::theme(plot.title = ggplot2::element_text(size=15,hjust = 0.5)))
#saving the plot if the user wants
ans7 <- readline(prompt=cat("\nwould you like to save this plot?\ny/n\n\n"))
if(ans7 == "y")
{
#save the plot
nm <- paste("ChIPseqVScontacts_correlation_plots_pScores_vs_peaks_",DandT2,".png",sep="")
wd <- as.numeric(readline(prompt=cat("\nenter the width size of the plot:\n\n")))
ht <- as.numeric(readline(prompt=cat("\nenter the height size of the plot:\n\n")))
cat("\nsaving the plot\nplease wait...\n\n")
ggplot2::ggsave(paste("~/Analyze4C/plots/",nm,sep=""),width=wd,height=ht)
save_flag <- 1
}
#intersected_reads_pairs_all <- as.data.frame(intersected_reads_pairs_all)
cat("\ncreating the reads VS. peaks plot\nplease wait...\n\n")
#print(ggplot2::ggplot(intersected_reads_pairs_all, ggplot2::aes(x=reads,y=peaks,color=names)) + ggplot2::scale_x_log10() + ggplot2::scale_y_log10() + ggplot2::geom_point() + ggplot2::geom_smooth(method=lm,se=FALSE,fullrange=TRUE))
#print(ggplot2::ggplot(intersected_reads_pairs_all, ggplot2::aes(x=reads,y=peaks,color=names)) + ggplot2::geom_point() + ggplot2::geom_smooth(method=lm,se=FALSE,fullrange=TRUE) + ggplot2::scale_y_continuous(limits = c(0,5000)) + ggplot2::ggtitle("contact VS. expression - reads correlation") + ggplot2::theme(plot.title = ggplot2::element_text(size=15,hjust = 0.5)))
print(ggplot2::ggplot(intersected_reads_pairs_all, ggplot2::aes(x=reads,y=peaks,color=names)) + ggplot2::geom_point() + ggplot2::geom_smooth(method=lm,se=FALSE,fullrange=TRUE) + ggplot2::ggtitle("contact VS. ChIPseq - reads correlation") + ggplot2::theme(plot.title = ggplot2::element_text(size=15,hjust = 0.5)))
#saving the plot if the user wants
ans10 <- readline(prompt=cat("\nwould you like to save this plot?\ny/n\n\n"))
if(ans10 == "y")
{
#save the plot
nm <- paste("ChIPseqVScontacts_correlation_plots_reads_vs_peaks_",DandT2,".png",sep="")
wd <- as.numeric(readline(prompt=cat("\nenter the width size of the plot:\n\n")))
ht <- as.numeric(readline(prompt=cat("\nenter the height size of the plot:\n\n")))
cat("\nsaving the plot\nplease wait...\n\n")
ggplot2::ggsave(paste("~/Analyze4C/plots/",nm,sep=""),width=wd,height=ht)
save_flag <- 1
}
if(save_flag == 1)
{
#replacing the " ; " at the start of the lines with nothing
CO_ps_type <- gsub("^ ; ","",CO_ps_type)
CO_ps <- gsub("^ ; ","",CO_ps)
CO_ps_chroms <- gsub("^ ; ","",CO_ps_chroms)
CO_peaks_type <- gsub("^ ; ","",CO_peaks_type)
CO_peaks <- gsub("^ ; ","",CO_peaks)
CO_peaks_chroms <- gsub("^ ; ","",CO_peaks_chroms)
if(choice4 == 2)
{
winds <- gsub("^ ; ","",winds)
}
#saving the parameters and details to ChIPseqVScontacts_correlation_plots.txt
ChIPseqVScontacts_correlation_plots[nrow(ChIPseqVScontacts_correlation_plots)+1,] <- c(psORreads,num_of_pScore_files,ps_Experiments,num_of_peaks_files,peaks_Experiments,CO_ps_type,CO_ps,CO_ps_chroms,CO_source,CO_peaks_type,CO_peaks,CO_peaks_chroms,int_chroms,secs_rmv,by_winds,winds,value_source,DandT1)
#sorting the list of experiments by bait alphabetically (and sorting the row indices)
ChIPseqVScontacts_correlation_plots <- ChIPseqVScontacts_correlation_plots[order(ChIPseqVScontacts_correlation_plots$Date_and_Time),]
rownames(ChIPseqVScontacts_correlation_plots) <- seq(length=nrow(ChIPseqVScontacts_correlation_plots))
#adding the new data to the file (by erasing the old one and creating a new one)
system("rm ChIPseqVScontacts_correlation_plots.txt")
write.table(ChIPseqVScontacts_correlation_plots,"ChIPseqVScontacts_correlation_plots.txt",sep="\t",row.names = FALSE,col.names = TRUE,quote=FALSE)
}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.