changeset 0:5558479694c7 draft

Uploaded script

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/graph_coverage_multi_png.R	Thu Jun 13 08:32:07 2024 +0000
@@ -0,0 +1,95 @@
+## 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"))
+}
+##-------------------------------------------------------