Mercurial > repos > kyost > atac_primer_tool

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

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
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])