# read information in the output file from hmmscan

data <- read.table("parse_hmmscan.out", sep = "\t", 
    header = TRUE, )

# return maximum sequence length

maxlen <- max(data$len)

# identify number of proteins

prots <- length(levels(data$protname))

plot(0, type = "n", xlim = c(-300, maxlen + 100), 
    ylim = c(0, prots + 2), xlab = "POS", ylab = "", yaxt = "n", 
    main = "Domain architecture of clotting proteins")

prevname <- ""
y <- 1
width <- 0.02
palette(rainbow(7))
lines <- length(data$protname)

for (i in (1:lines)) {
    
    test <- data$protname[i]
    if (test != prevname) {
        y <- y + 1
        # print protein names at y axis
        text(-10, y, data$protname[i], cex = 0.7, pos = 2)
        # draw lines corresponding to the protein lengths
        lines(c(1, data$len[i]), c(y, y), col = "grey", lw = 2)
    }
    prevname <- test
    
    
    # draw domain rectangles
    domlen <- length(levels(data$domname))
    # assign a colour to the domain
    for (k in (1:domlen)) {
        if (levels(data$domname)[k] == data$domname[i]) {
            color <- k
        }
    }
    ybottom <- y - width * color
    ytop <- y + width * color
    rect(data$begin[i], ybottom, data$end[i], ytop, col = color)
}

# Finally draw the domain information in a separate panel

x <- maxlen * 0.7
for (k in (1:domlen)) {
    pos <- 4 + k * 0.8
    rect(x, pos - 0.02 * k, x + 100, pos + 0.02 * k, col = k)
    text(x + 100, pos, levels(data$domname)[k], pos = 4, cex = 0.7)
    
}