extract_3C_profile.R

#!/bin/env Rscript

# Initialize the variables
library_path <- "RLIB"
paired_bed_file <- "test_data/SRR292238X.reassigned.array.bed"
fragment <- "chr14:54607949-54608839"
outfile <- "test_data/chr4_54607949_54608839.bedgraph"
verbose <- FALSE

# Load the libraries
require(getopt, quietly = TRUE)

# set the commandline option specifications
specification <- matrix(c(
	"matrix", "m", 1, "character",
	"fragment", "f", 1, "character",
	"outfile", "o", 1, "character",
	"library", "l", 1, "character",
	"verbose", "v", 0, "logical",
	"help", "h", 0, "logical"
), byrow = TRUE, ncol = 4)

# parse the commandline options
commandline_options <- getopt(specification)

# print the usage if help was asked
if(!is.null(commandline_options[["help"]])) {
	cat(getopt(specification, usage = TRUE))
	q(status = 1)
}
if(!is.null(commandline_options[["verbose"]])) {
	verbose <- TRUE	
}
if(!is.null(commandline_options[["matrix"]])){
	paired_bed_file <- commandline_options[["matrix"]]
}
if(!is.null(commandline_options[["fragment"]])){
	fragment <- commandline_options[["fragment"]]
}
if(!is.null(commandline_options[["outfile"]])){
	outfile <- commandline_options[["outfile"]]
}
if(!is.null(commandline_options[["library"]])){
	library_path <- commandline_options[["library"]]
}
if(verbose){
	message(paste("Matrix bed file:", paired_bed_file))
	message(paste("Fragment:", fragment))
	message(paste("Output file:", outfile))
	message(paste("Library path:", library_path))
}

# Load the libraries
#
#
source(file.path(library_path, "bed_readers.R"))

# Load the input data
# 
paired_bed_data <- read_paired_bed(paired_bed_file)

key_pairs <- with(paired_bed_data, cbind(
	paste(chr_a, ":", start_a, "-", end_a, sep=""),
	paste(chr_b, ":", start_b, "-", end_b, sep="")
))

# Get the rows with fragment equal to specified fragment
# 
tf_a <- key_pairs[,1] == fragment
tf_b <- key_pairs[,2] == fragment

# stop if no fragments
if(sum(tf_a) + sum(tf_b) == 0){
	stop(paste0("No interactions found for fragment ", fragment))
}

# Extract the targets
#
x <- paired_bed_data[tf_a, c( 4, 5, 6, 8)]
colnames(x) <- c("chr", "start", "end", "signal")
y <- paired_bed_data[tf_b, c( 1, 2, 3, 8)]
colnames(y) <- c("chr", "start", "end", "signal")

# Reorder the extracted data
# 
d <- rbind(x, y)
i <- order(d$chr, d$start, d$end)
d <- d[i,]

# Write the output
# 
write.table(d, file = outfile, sep = "\t", row.names = FALSE, col.names = FALSE, quote = FALSE)