comparison graph_coverage_multi_png.R @ 0:5558479694c7 draft

Uploaded script

0:5558479694c7

## draw coverage depth graphics for n references
## if n <= 8: return png
## if n >  8: return pdf

args = commandArgs(trailingOnly=TRUE)
library("reshape")
library("ggplot2")
library("gridExtra")  # Load some extensions to have multiplot pdf
library("stringr")    # To replace characters in strings

cov_file = args[1]
output = args[2]
coverage = read.table(cov_file,sep="\t", header=FALSE)
coverage = rename(coverage,c(V1="Chr", V2="locus",V3="depth"))

headers = unique(coverage$Chr)

##         round sup(length(each different ref)) / 8
## nb_pages = ceiling(length(headers) / 8)
nb_pages = floor(length(headers) / 8)

print("nb_images:")
print(nb_pages)

output_p = ''

## -------------------------------------------------------
## plot each page (i more than 1 page)
if(nb_pages > 0){
    for(n in 1:nb_pages-1){
        print(n) ## start with 0
        
        ## creates output file
        if(n == 0){
            output_p <- output
        }else{
            output_p <- str_replace(output, ".png", paste0("_", n ,".png") )
        }
        png(output_p)
        print(paste0("output file:", output_p))
        header_start_index = n*8 +1
        header_end_index = header_start_index + 7
        treated_headers = headers[header_start_index:header_end_index]
        coverage_eight_headers <- coverage[coverage$Chr %in% treated_headers,]
        print("header_start_index:")
        print(header_start_index)
        print("header_end_index:")        
        print(header_end_index)
        print("treated_headers:")
        print(treated_headers)
        maxlocus = max( coverage_eight_headers$locus )
                                        # # maxdepth = max( coverage_eight_headers$depth )
        print("maxlocus:")
        print(maxlocus)
        print(paste0(output_p," file will be created"))
        plot8cov <- ggplot(coverage_eight_headers, aes(x=locus, y=depth )) + geom_line(color="red", size=1) + facet_wrap(~coverage_eight_headers$Chr, scales="free", ncol=2, nrow=4) + scale_x_continuous(expand = c(0, 0), limits = c(0, NA))
        print(plot8cov)
        dev.off()
        print(paste0(output_p," file created"))
    }
}
## creates output file
if(nb_pages == 0){
    output_p <- output
}else{
    output_p <- str_replace(output, ".png", paste0("_", nb_pages ,".png") )
}
png(output_p)

nb_elem_last_page = length(headers) %% 8
if(nb_elem_last_page > 0){
    header_start_index = length(headers)-nb_elem_last_page+1
    header_end_index = length(headers)
    treated_headers = headers[header_start_index:header_end_index]
    print("treated_headers:")
    print(treated_headers)    
    coverage_last_headers= coverage[coverage$Chr %in% treated_headers,]
    print("header_start_index:")
    print(header_start_index)
    print("header_end_index:")        
    print(header_end_index)
    print("treated_headers:")
    print(treated_headers)
    maxlocus = max( coverage_last_headers$locus )
    ## maxdepth = pmax( coverage_last_headers$depth )
    print("maxlocus:")
    print(maxlocus)
    print(paste0(output_p," file will be created"))

    plot8cov <- ggplot(coverage_last_headers, aes(x=locus, y=depth )) + geom_line(color="red", size=1) + facet_wrap(~coverage_last_headers$Chr, scales="free", ncol=2, nrow=4) + scale_x_continuous(expand = c(0, 0), limits = c(0, NA))
    print(plot8cov)
    dev.off()
    print(paste0(output_p," file created"))
}
##-------------------------------------------------------