Load files
library(gtools)
kar <- read.table("~/Downloads/karyotype.human.hg19.txt")
head(kar, 2)
## V1 V2 V3 V4 V5 V6 V7
## 1 chr - hs1 hs1 0 249250621 chr1
## 2 chr - hs2 hs2 0 243199373 chr2
# this file was produce from a bam file using genomeCoverageBed (BedTools)
filebed <- c("~/Downloads/accepted_hits.uni.clus.prot.sat.bg")
bed <- read.table(filebed)
head(bed, 2)
## V1 V2 V3 V4 V5 V6
## 1 chr1 16218 16818 4 0 0
## 2 chr1 47080 47152 2 0 0
Prepare chromosomes
# get chromosomes
chr <- kar[kar$V1 == "chr", ]
# get bands
band <- kar[kar$V1 == "band", ]
# order chromosomes from chr1...chrY
band <- band[mixedorder(band$V2), ]
# index chromosome names to position in the graphic
bandchr <- c("hs1", "hs2", "hs3", "hs4", "hs5", "hs6", "hs7", "hs8",
"hs9", "hs10", "hs11", "hs12", "hs13", "hs14", "hs15", "hs16", "hs17", "hs18",
"hs19", "hs20", "hs21", "hs22", "hsx", "hsy")
bandpos <- seq(1, 96, 4)
names(bandpos) <- bandchr
# index bands to colors
bandtype <- as.character(unlist(unique(band$V7)))
bandcolor <- c("grey77", "grey75", "grey65", "grey55", "grey40",
"red3", "grey40", "grey55")
names(bandcolor) <- bandtype
# index bands to thickness (centromeres are thiner)
bandthin <- c(1, 1, 1, 1, 1, 0.1, 1, 1)
names(bandthin) <- bandtype
Create index between chromosome name and position on the y axis
bedchr <- c("chr1", "chr2", "chr3", "chr4", "chr5", "chr6", "chr7",
"chr8", "chr9", "chr10", "chr11", "chr12", "chr13", "chr14", "chr15", "chr16",
"chr17", "chr18", "chr19", "chr20", "chr21", "chr22", "chrX", "chrY")
bedpos <- bandpos
names(bedpos) <- bedchr
Some constants to plot the RNA-seq data
startpos <- 0.3 #plot the fragment over the chromosomes
vectors <- c(-1, 0.6, 1.5) #threshold to rank RNA-expression
blue <- "blue4" #fragment highly expressed will be blue
grey75 <- rgb(0.9, 0.9, 0.9, 0.5) #fragment low expressed will be transparent grey
Prepare the RNAseq coordinates
bed$V4 <- log10(bed$V4) #logarithm of the fragment coverage
maxe <- max(bed$V4) #get the maximum fragment coverage
ratio <- bed$V4/maxe #normalized by the maximum fragment coverage
ratio.high <- ratio[ratio >= 0.4] #remove low expressed genes
tmp <- bed[ratio >= 0.4, ] #remove low expressed genes
ratiocol <- cut(ratio.high, breaks = vectors, labels = c(grey75,
blue)) #assign colors
plot RNA-seq data and chromosomes
par(mar = c(1, 4, 1, 1))
plot(1, yaxt = "n", xaxt = "n", type = "n", xlim = c(0, max(band$V6)),
ylim = c(0, 94), xlab = "", ylab = "", lwd = 0, lty = 0, axes = FALSE) #create plot
axis(2, at = bandpos, labels = bedchr, las = 2, tick = 0, line = 0,
lty = 0, lwd = 0) #add axis
rect(band$V5, bandpos[as.character(band$V2)] - 0.2, band$V6, bandpos[as.character(band$V2)] +
0.2, col = bandcolor[as.character(band$V7)], lty = 0, lwd = bandthin[as.character(band$V7)]) #add chromosomes
rect(tmp$V2, bedpos[as.character(tmp$V1)], tmp$V3, bedpos[as.character(tmp$V1)] +
startpos + ratio.high, col = as.character(ratiocol), lty = 1, border = as.character(ratiocol),
lwd = 0.9) #add RNA-seq data