Invest some time to learn the R-package ggplot2. I made an example for you below:

# create copy number data for 4 samples for 3 genes 
data<-data.frame(
    gene=sample(c("geneA","geneB","geneC"),10000,replace=TRUE),
    sampleName=sample(c("sample1","sample2","sample3","sample4"),10000,replace=TRUE),
    position=round(runif(n=10000,min=0,max=5000),digits=0),
    copy_number=rnorm(n=10000,mean=2,sd=0.1)
)

# create copy number gain in sample 4 gene A
data[data$gene=="geneA" & data$sampleName=="sample1",4]<-rnorm(
    mean=3,
    sd=0.1,
    n=length(
        data[data$gene=="geneA" & data$sampleName=="sample1",4]
    )
)

# and a focal loss at sample 3, gene B
data[data$gene=="geneB" & data$sampleName=="sample3" & data$position > 2000 & data$position < 4000,4]<-rnorm(
    mean=1,
    sd=0.1,
    n=length(
        data[data$gene=="geneB" & data$sampleName=="sample3" & data$position > 2000 & data$position < 4000,4]
    )
)

# plot copy numbers for all genes for all samples
library(ggplot2)
plot<-ggplot(data) + #add basic layer
geom_point(aes(x=position,y=copy_number)) + #add copy number points to plot
facet_wrap(gene ~ sampleName) + # seperate the CN numbers for each gene and for each sample
scale_y_continuous(limits=c(0,4)) # set y-limits

# and save plot to file
png("copy_numbers.png",1000,1000)
print(plot)
dev.off()

The fruit of my own searching, in parallel to asking this question is a tool called Dpeaks.

It is not exactly what i am after, but it produces plots (in various formats, drawn by R), based on user supplied variable-WIG file(s).

There is one Perl script to run the process. Another Perl script all the conversion of fixed-WIG to variable-WIG. By trial and error i found you can also start of with a BedGraph file and convert it, using BedToWig.sh to a variable-WIG file, see HERE for more information.

BTW MACS output the required variable-WIG format.

What we do here is to plot a heatmap where rows represents each regions and each column is a bin within the regions. This is very common method that you often see in the the literature.

Normally we use the ratio between the ChIP-seq data and the input, otherwise you may see the biases instead of the signal.

We use our own software for this but I will be happy to share it if you are interested.