comparison find_qPCR_regions.R @ 5:3cd53127a838 draft default tip

planemo upload commit 103cb51ec368438642504c3f98b76c363db478bb
author kyost
date Fri, 04 May 2018 02:24:47 -0400
parents fbfe7b785ea7
children
comparison
equal deleted inserted replaced
4:72571a30f17b 5:3cd53127a838
1 ## Command to run tool: 1 ## Command to run tool:
2 # Rscript --vanilla find_qPCR_regions.R o.coverage.bed f9.coverage.bed lib_sizes.txt cor_cutoff cov_cutoff output_file 2 # Rscript --vanilla find_qPCR_regions.R o.coverage.bed f9.coverage.bed lib_sizes.txt cor_cutoff cov_cutoff output_file
3 3
4 # Set up R error handling to go to stderr 4 options(warn=-1)
5 options(show.error.messages=F, error=function(){cat(geterrmessage(),file=stderr());q("no",1,F)})
6 5
7 # Avoid crashing Galaxy with an UTF8 error on German LC settings 6 # Avoid crashing Galaxy with an UTF8 error on German LC settings
8 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8") 7 loc <- Sys.setlocale("LC_MESSAGES", "en_US.UTF-8")
9 8
10 args <- commandArgs(TRUE) 9 args <- commandArgs(TRUE)
56 m <- ncol(asplit[[i]]) 55 m <- ncol(asplit[[i]])
57 56
58 a_split <- asplit[[i]] 57 a_split <- asplit[[i]]
59 b_split <- bsplit[[i]] 58 b_split <- bsplit[[i]]
60 59
61 if (sd(t(a_split[1,5:m])) == 0 || sd(t(b_split[-1,5:m])) == 0 ){
62 next
63 }
64 60
65 #calculate correlation between total peak reads and spanning reads in each window 61 #calculate correlation between total peak reads and spanning reads in each window
66 corb <- data.frame(cor(t(a_split[1,5:m]) 62 corb <- data.frame(cor(t(a_split[1,5:m])
67 , t(b_split[-1,5:m]) 63 , t(b_split[-1,5:m])
68 , method = "pearson")) 64 , method = "pearson"))
81 colnames(data2)<-c("position","averageReadDepth") 77 colnames(data2)<-c("position","averageReadDepth")
82 colnames(data3)<-c("position","averageReadDepth") 78 colnames(data3)<-c("position","averageReadDepth")
83 79
84 #keep windows for which correlation with total peak reads is above threshold and 80 #keep windows for which correlation with total peak reads is above threshold and
85 #number of spanning reads is above cutoff 81 #number of spanning reads is above cutoff
86 regions <- data[which(data$correlation > corr_cut & data3$averageReadDepth > cov_cut),] 82 regions <- data[which(data$correlation >= corr_cut & data3$averageReadDepth >= cov_cut),]
87 regions <- regions[,2:4] 83 regions <- regions[,2:4]
88 84
89 #collapse windows to non-overlappin regions 85 #collapse windows to non-overlappin regions
90 if (nrow(regions)>0){ 86 if (nrow(regions)>0){
91 regions$V4 <- paste(as.character(a_split[1,4]),"_region1",sep="") 87 regions$V4 <- paste(as.character(a_split[1,4]),"_region1",sep="")
92 collapsed_regions <- regions[1,] 88 collapsed_regions <- regions[1,]
93 nregion <- 1 89 nregion <- 1
94 if (nrow(regions)>=2){ 90 if (nrow(regions)>=2){
95 for (j in 2:nrow(regions)) { 91 for (j in 2:nrow(regions)) {
96 if (as.numeric(regions[j,2]) <= (as.numeric(regions[j-1,3]))) { 92 if (as.numeric(regions[j,2]) < (as.numeric(regions[j-1,3]))) {
97 collapsed_regions[nrow(collapsed_regions),3] <- regions[j,3] 93 collapsed_regions[nrow(collapsed_regions),3] <- regions[j,3]
98 }else{ 94 }else{
99 collapsed_regions <- rbind(collapsed_regions, regions[j,]) 95 collapsed_regions <- rbind(collapsed_regions, regions[j,])
100 nregion <- nregion + 1 96 nregion <- nregion + 1
101 collapsed_regions[nrow(collapsed_regions),4] <- paste(as.character(a_split[1,4]) 97 collapsed_regions[nrow(collapsed_regions),4] <- paste(as.character(a_split[1,4])