bitlab: chromeister/bin/compute

comparison chromeister/bin/compute_score.R @ 0:7fdf47a0bae8 draft

Uploaded

author	alvarofaure
date	Wed, 12 Dec 2018 07:18:40 -0500
parents
children

comparison

equal deleted inserted replaced

--1:000000000000
+:7fdf47a0bae8
+#!/usr/bin/env Rscript
+args = commandArgs(trailingOnly=TRUE)
+if(length(args) < 2){
+stop("USE: Rscript --vanilla plot.R <matrix> <matsize>")
+}
+growing_regions <- function(mat, reward = 6, penalty = 5, sidePenalty = 3, MAXHSPS = 500, TH = 10, WSIZE = 7){
+#write(mat, file ="data.txt", ncolumns = 200)
+len <- min(length(mat[1,]), length(mat[,1]))
+HSPS <- matrix(0, nrow=MAXHSPS, ncol=5)
+idx <- 1
+lH <- round(WSIZE/2) - 1
+rH <- round(WSIZE/2) + 1
+#print(paste("lH: ", lH, "rH: ", rH))
+if(WSIZE %% 2 == 0){
+print("WSIZE MUST BE ODD")
+stop()
+}
+i <- 1
+#readline(prompt="Press [enter] to continue")
+while(i < len){
+value <- max(mat[i,]) * reward
+if(value == 0) i <- i + 1
+pos <- which.max(mat[i,])
+# these two hold ending frag
+endfrag <- pos
+j <- i
+count_penalties <- 1
+#print("-----------------")
+while(value > 0 && j < len){
+#print(paste("took values ", j, endfrag, "which have score of ", mat[j, endfrag], "current value is: ", value))
+# Reset position used
+mat[max(1,j-1), max(1,endfrag-2)] <- 0
+mat[max(1,j-1), max(1,endfrag-1)] <- 0
+mat[max(1,j-1), endfrag] <- 0
+mat[max(1,j-1), min(len, endfrag+1)] <- 0
+mat[max(1,j-1), min(len, endfrag+2)] <- 0
+mat[j, max(1,endfrag-2)] <- 0
+mat[j, max(1,endfrag-1)] <- 0
+mat[j, endfrag] <- 0
+mat[j, min(len, endfrag+1)] <- 0
+mat[j, min(len, endfrag+2)] <- 0
+#print(paste("Erasing: (",max(1,j-1), max(1,endfrag-2),")(",max(1,j-1), max(1,endfrag-1),")(",max(1,j-1), endfrag,
+#            ")(",max(1,j-1), min(len, endfrag+1),")(",max(1,j-1), min(len, endfrag+2),")(",j, max(1,endfrag-2),
+#            ")(",j, max(1,endfrag-1),")(",j, endfrag,")(",j, min(len, endfrag+1),")(",j, min(len, endfrag+2),")"))
+# Go for next
+j <- j + 1
+# Check next, could be reverse or forward
+mleft <- max(1, endfrag-lH)
+mright <- min(len, endfrag+lH)
+window <- mat[j, mleft:mright]
+#print(paste("mleft: ", mleft, "mright", mright, "j is: ", j))
+#print(window)
+v <- max(window)
+selected <- which.max(window)
+# Make it rather go diagonally
+#print(paste("WIndow len: ", length(window)))
+chose_diagonal <- FALSE
+if(length(window) == WSIZE && v == window[lH]){
+selected <- lH
+chose_diagonal <- TRUE
+}
+if(length(window) == WSIZE && v == window[rH]){
+selected <- rH
+chose_diagonal <- TRUE
+}
+#print(paste("Selected value be like ", selected))
+if(v != 0){
+endfrag <- (mleft + selected ) # To make the indexing
+if(length(window) == WSIZE) endfrag <- endfrag - 1
+endfrag <- max(1, min(len, endfrag))
+#print(paste("\t", " endfragnew = ", endfrag, "max of window: ", max(window), "on position", which.max(window)))
+}
+#print(paste("Chose diagonal is ", chose_diagonal))
+# If no similarity is found
+if(v == 0){
+value <- value - count_penalties * penalty
+count_penalties <- count_penalties + 1
+#print("Got penalty #########")
+}else{
+# Similarity is found
+if(!chose_diagonal){
+value <- value + count_penalties * (-sidePenalty)
+count_penalties <- count_penalties + 1
+#print("Got SIDE penalty @@@ #########")
+}else{
+count_penalties <- 1
+value <- value + reward
+#print("Got reward thou #########")
+}
+}
+#readline(prompt="Press [enter] to continue")
+}
+# Check len of frag
+if(j-i > TH){
+# HSPS[idx, 1] <- i
+# HSPS[idx, 2] <- pos
+# HSPS[idx, 3] <- j
+# HSPS[idx, 4] <- endfrag
+# HSPS[idx, 5] <- abs(i-j)
+HSPS[idx, 1] <- pos
+HSPS[idx, 2] <- i
+HSPS[idx, 3] <- endfrag
+HSPS[idx, 4] <- j
+HSPS[idx, 5] <- abs(i-j)
+#print(paste("ACCEPT", i, pos, j, endfrag, "x-y", j-i, abs(endfrag-pos), sep = " "))
+idx <- idx + 1
+}else{
+#print("REJECT")
+}
+# This will prevent overlappign lines, I think
+#i <- j
+if(idx == MAXHSPS) break
+}
+return (HSPS)
+}
+detect_events <- function(HSPS, sampling){
+DIAG_SEPARATION <- 10
+# same as HSPS but adding the event
+output <- matrix(0, nrow=length(HSPS[,1]), ncol=1+length(HSPS[1,]))
+colnames(output) <- c("x1", "y1", "x2", "y2", "len", "event")
+j <- 0
+for(i in 1:(length(HSPS[,1]))){
+if(sum(HSPS[i,]) > 0){
+j <- j + 1
+is_inverted = FALSE
+is_diagonal = TRUE
+if(HSPS[i,1] > HSPS[i,3]) is_inverted = TRUE
+if(abs(HSPS[i,1] - HSPS[i,2]) > DIAG_SEPARATION && abs(HSPS[i,3] - HSPS[i,4]) > DIAG_SEPARATION) is_diagonal = FALSE
+output[i,1] <- HSPS[i,1] * sampling
+output[i,2] <- HSPS[i,2] * sampling
+output[i,3] <- HSPS[i,3] * sampling
+output[i,4] <- HSPS[i,4] * sampling
+output[i,5] <- HSPS[i,5] * sampling
+if(is_diagonal) output[i,6] <- "synteny block"
+if(is_diagonal && is_inverted) output[i,6] <- "inversion"
+if(!is_diagonal && !is_inverted) output[i,6] <- "transposition"
+if(!is_diagonal && is_inverted) output[i,6] <-"inverted transposition"
+}
+}
+return (output[1:j,])
+}
+paint_frags <- function(HSPS, l, sampling){
+plot(c(HSPS[1,1]*sampling, HSPS[1,3]*sampling), c(HSPS[1,2]*sampling, HSPS[1,4]*sampling), xlim = c(1,l*sampling), ylim = c(1,l*sampling), type="l", xlab="X-genome",ylab="Y-genome")
+for(i in 2:length(HSPS[,1])){
+if(sum(HSPS[i,]) > 0){
+lines(c(HSPS[i,1]*sampling, HSPS[i,3]*sampling), c(HSPS[i,2]*sampling, HSPS[i,4]*sampling))
+}
+}
+}
+supersample <- function(mat, upscale){
+l <- min(length(mat[1,]), length(mat[,1]))
+size <- round(l*upscale)
+m <- matrix(0, nrow=size, ncol=size)
+for(i in 1:size){
+for(j in 1:size){
+mleft <- max(1, i-1)
+mright <- min(l, i+1)
+mup <- max(1, j-1)
+mdown <- min(l, j+1)
+ri <- max(1, i/upscale)
+rj <- max(1, j/upscale)
+if(mat[ri, rj] > 0){
+m[i, j] <- 1
+}
+}
+}
+return (m)
+}
+downsample <- function(mat, downscale){
+l <- min(length(mat[1,]), length(mat[,1]))
+size <- round(l/downscale)
+m <- matrix(0, nrow=size, ncol=size)
+for(i in 1:l){
+for(j in 1:l){
+mup <- max(1, i-1)
+mdown <- min(l, i+1)
+mleft <- max(1, j-1)
+mright <- min(l, j+1)
+#print(paste("Matrix at ", i, j))
+#print(mat[mup:mdown, mleft:mright])
+if(sum(mat[mup:mdown, mleft:mright]) > 0){
+#print(paste("goes to ", round(i/downscale), round(j/downscale)))
+m[round(i/downscale), round(j/downscale)] <- 1
+}
+}
+}
+return (m)
+}
+path_mat = args[1]
+fancy_name <- strsplit(path_mat, "/")
+fancy_name <- (fancy_name[[1]])[length(fancy_name[[1]])]
+# Read sequence lenghts
+con <- file(path_mat,"r")
+seq_lengths <- readLines(con, n=2)
+seq_x_len <- as.numeric(seq_lengths[1])
+seq_y_len <- as.numeric(seq_lengths[2])
+close(con)
+data <- as.matrix(read.csv(path_mat, sep = " ", skip=2))
+# Get sequence names
+name_x <- strsplit(fancy_name, "-")[[1]][1]
+name_y <- strsplit(fancy_name, "-")[[1]][2]
+# Max of columns
+len_i <- as.numeric(args[2])
+len_j <- as.numeric(args[2])
+score_density <- data
+aux_density <- data
+pmax_pos <- which.max(aux_density[,1])
+for(i in 1:len_i){
+cmax_pos <- which.max(aux_density[i,]) # get max of row
+if((aux_density[i,cmax_pos]) > 0){ # if it has value
+#row_from_col <- which.max(aux_density[,cmax_pos]) # get max of the column pointed by maximum of row
+score_density[i,] <- 0 # put all others in row to 0 i.e. only use this max, EXCEPT for the maximum in the column
+#score_density[row_from_col, cmax_pos] <- 1
+score_density[i,cmax_pos] <- 1
+pmax_pos <- cmax_pos
+}
+}
+pmax_pos <- which.max(aux_density[1,])
+for(i in 1:len_i){
+cmax_pos <- which.max(aux_density[,i]) # get max of column
+if((aux_density[cmax_pos,i]) > 0){
+score_density[,i] <- 0
+score_density[cmax_pos,i] <- 1
+pmax_pos <- cmax_pos
+}
+}
+score_copy <- score_density
+diag_len <- 4
+for(i in 6:(len_i-5)){
+for(j in 6:(len_j-5)){
+value <- 0
+for(w in (-diag_len/2):(diag_len/2)){
+if(score_density[i+w,j+w] > 0){
+value <- value + 1
+}
+}
+if(value >= diag_len){
+for(k in 1:5){
+score_copy[i+k,j+k] <- 1
+}
+for(k in 1:5){
+score_copy[i-k,j-k] <- 1
+}
+}else if(score_copy[i,j]==0){
+score_copy[i,j] <- 0
+}
+}
+}
+for(i in 6:(len_i-5)){
+for(j in 6:(len_j-5)){
+value <- 0
+for(w in (-diag_len/2):(diag_len/2)){
+if(score_density[i-w,j+w] > 0){
+value <- value + 1
+}
+}
+if(value >= diag_len){
+for(k in 1:5){
+score_copy[i-k,j+k] <- 1
+}
+for(k in 1:5){
+score_copy[i+k,j-k] <- 1
+}
+}else if(score_copy[i,j]==0){
+score_copy[i,j] <- 0
+}
+}
+}
+# Kernel to remove single points
+for(i in 1:(length(score_copy[,1]))){
+for(j in 1:(length(score_copy[1,]))){
+value <- 0
+min_i <- max(1, i-1)
+max_i <- min(len_i, i+1)
+min_j <- max(1, j-1)
+max_j <- min(len_j, j+1)
+value <- sum(score_copy[min_i:max_i, min_j:max_j])
+if(value < 2) score_copy[i,j] <- 0
+}
+}
+# To compute the score
+score <- 0
+pmax_pos <- which.max(score_copy[,1])
+dist_th <- 1.5
+besti <- 1
+bestj <- pmax_pos
+dvec1 <- abs(which.max(score_copy[,2]) - which.max(score_copy[,1]))
+dvec2 <- abs(which.max(score_copy[,3]) - which.max(score_copy[,2]))
+dvec3 <- abs(which.max(score_copy[,4]) - which.max(score_copy[,3]))
+for(i in 5:len_i){
+#print(paste(paste(paste(dvec1, dvec2), dvec3), mean(c(dvec1, dvec2, dvec3))))
+distance <- mean(c(dvec1, dvec2, dvec3))
+dvec1 <- dvec2
+dvec2 <- dvec3
+dvec3 <- abs(which.max(score_copy[,i]) - which.max(score_copy[,i-1]))
+# If there is a 0 or we are too far away just add max distance!
+if(distance > dist_th || distance == 0){
+score <- score + len_i
+}
+}
+score <- (score/(len_i^2))
+print(score)
+sampling_value <- 5
+submat <- downsample(score_copy, sampling_value)
+m <- growing_regions((submat), WSIZE = 7, TH = 5, penalty = 15)
+events <- detect_events(m, sampling_value)
+events <- rbind(events, c(0,0,0,0,0,"Null event"))
+write(as.character(c(seq_x_len, seq_y_len)), file=paste(path_mat,".events.txt", sep=""), append = TRUE, sep =",", ncolumns=2)
+write(as.character(c("x1", "y1", "x2", "y2", "len", "event")), file=paste(path_mat,".events.txt", sep=""), append = TRUE, sep =",", ncolumns=6)
+for(i in 1:length(events[,1])){
+write(as.character(events[i,]), file=paste(path_mat,".events.txt", sep=""), append = TRUE, sep =",", ncolumns=6)
+}
+coords1 <- round(seq(from=0, to=1, by=0.2)*seq_x_len)
+coords2 <- round(seq(from=0, to=1, by=0.2)*seq_y_len)
+final_image <- apply((t(score_copy)), 2, rev)
+png(paste(path_mat, ".filt.png", sep=""), width = length(data[,1]), height = length(data[,1]))
+image(t(final_image), col = grey(seq(1, 0, length = 256)), xaxt='n', yaxt='n', main = paste(fancy_name, paste("filt. score=", score)), xlab = name_x, ylab = name_y, axes = FALSE)
+axis(1, tick = TRUE, labels = (coords1), at = c(0.0, 0.2, 0.4, 0.6, 0.8, 1))
+axis(2, tick = TRUE, labels = rev(coords2), at = c(0.0, 0.2, 0.4, 0.6, 0.8, 1))
+dev.off()
+# Output pixel coordinates of highly conserved signals
+# To clear it in case it existed
+write(c(paste("X", "Y")), file = paste("hits-XY-", paste(fancy_name, ".hits", sep=""), sep=""), append = FALSE, sep = "\n")
+for(i in 1:(length(score_copy[,1]))){
+for(j in 1:(length(score_copy[1,]))){
+if(score_copy[i,j] != 0){
+write(c(paste(i, j)), file = paste("hits-XY-", paste(fancy_name, ".hits", sep=""), sep=""), append = TRUE, sep = "\n")
+}
+}
+}

Mercurial > repos > alvarofaure > bitlab

comparison chromeister/bin/compute_score.R @ 0:7fdf47a0bae8 draft