barplot-heatmap-visualization.R

# cmdline arguments to be given in the following order
# 1. inputFilename
# 2. width
# 3. height
# 4. outputFilename (csv format)
# 5. # of sequences to visualize 
# (This is useful because there could be hugely many sequences in the 
# collection and visualizing them all together may not be a good idea!)
#
# author: snikumbh@mpi-inf.mpg.de

args <- commandArgs(TRUE)
inputFilename <- args[1]
widthVal <- args[2]
heightVal <- args[3]
outputFilename <- args[4]
numOfSequences <- args[5]

# If a file with this name already exists, this overwrites it!
pdf('barplot-visualization-variable-len-seq.pdf', width=10.0, height=6.0)
csv_filename <- paste0(outputFilename, '.csv')
flds <- count.fields(csv_filename, sep=',') # get the number of segments in all sequences
seqs <- read.csv(csv_filename, fill=TRUE, nrows=length(flds), header=F, sep=",", col.names=paste0("Seg-", seq(1,max(flds),by=1)))
seqSegs <- seqs[1:numOfSequences,]
seqSegs_ones <- seqSegs
seqSegs_ones[!is.na(seqSegs)] <- 1 # sequences which had less than the maximum number of segments for any sequence in the collection, will have NA

ncolors <- max(flds)
color_palette <- heat.colors(ncolors)

# If a different color palette is required could use the following with experimenting:
# jet.colors <- colorRampPalette(c("#00007F", "blue", "#007FFF", "cyan",
#                     "#7FFF7F", "yellow", "#FF7F00"))#, "red", "#7F0000"))

ypos <- apply(seqSegs, 2, cumsum)
ypos <- ypos - seqSegs/2

xpos <- barplot(t(seqSegs_ones), col=c("white"), axes=F, horiz=T, axisnames=F)
cp <- 0
for (i in 1:ncol(seqSegs_ones)){
    temp <- color_palette[seqSegs[seq(1,nrow(seqSegs),by=1), i]]
    cp <- cbind(cp, temp)
}
new_cp <- cp[,-1]

for (i in 1:nrow(seqSegs_ones)){
    barplot(t(seqSegs_ones[i,]), col=new_cp[i,], horiz=T, add=T, space=(i-1)*1.2, axes=F, xlab="", axisnames=F)
}

axis(1,at=0.5+seq(0,ncolors-1,1), labels=seq(1,ncolors,1))
# Additional axis showing the nucleotide ranges, if required
# axis(1,at=seq(0,ncolors-1,1), labels=seq(1,ncolors,1))
legend("topright", legend = seq(1,ncolors), fill = color_palette, title="Segment ranks", ncol = 2, x.intersp = 0.1, bty="n", cex=0.8)

dev.off()