diff --git a/chromatoR/app.R b/chromatoR/app.R index 817b2b1..d9705fe 100644 --- a/chromatoR/app.R +++ b/chromatoR/app.R @@ -1,4 +1,7 @@ library(shiny) +library(tidyverse) +library(sangerseqR) +library(msa) # Define UI for application that draws a histogram ui <- fluidPage( @@ -9,31 +12,187 @@ ui <- fluidPage( # Sidebar with a slider input for number of bins sidebarLayout( sidebarPanel( - sliderInput("bins", - "Number of bins:", - min = 1, - max = 50, - value = 30) + fileInput("file1", "Sube fichero ab1", multiple = FALSE), + numericInput("thr", label = "Peak Threshold", value = 450), + numericInput("dist1", label = "1st Distance Reduction", value = 4), + numericInput("dist2", label = "2nd Distance Reduction", value = 7), + numericInput("ratio", label = "Ratio 2ary seq", value = 0.33), + textInput("old", label="Referencia"), + actionButton("calab", "Analizar") ), # Show a plot of the generated distribution mainPanel( - plotOutput("distPlot") + tags$head(tags$style(HTML("pre,.wrap { white-space: pre-wrap; word-break: break-all; }"))), + verbatimTextOutput("Primseq") %>% tagAppendAttributes(class="wrap"), + verbatimTextOutput("align") ) ) ) # Define server logic required to draw a histogram server <- function(input, output) { + # observe({ + # + # obj<<-readsangerseq(input$file1) + # } + # }) + obj<-reactiveValues() + obj$aborig<-NULL + obj_ab<-NULL + obj$seq<-NULL + + observe({ + if (!is.null(input$file1)){ + obj$aborig<-readsangerseq(input$file1$datapath) + } + }) + + + observeEvent(input$calab, { + if (!is.null(obj$aborig)){ + obj_ab<-obj$aborig + ## Functions + getpeaks <- function(trace) { + r <- rle(trace) + indexes <- which(rep(diff(sign(diff(c(-Inf, r$values, -Inf)))) == -2, + times = r$lengths)) + cbind(indexes, trace[indexes]) + } - output$distPlot <- renderPlot({ - # generate bins based on input$bins from ui.R - x <- faithful[, 2] - bins <- seq(min(x), max(x), length.out = input$bins + 1) + peakvalues <- function(x, pstart, pstop) { + region <- x[x[,1] > pstart & x[,1] < pstop, ,drop=FALSE] + if (length(region[,1]) == 0) return(c(0, NA)) + else return(c(max(region[,2], na.rm=TRUE), region[which.max(region[,2]),1])) + } + + redDiffs<-function(matrix, diff=5){ + diffs<-diff(matrix[,1]) + pos<-matrix(nrow=0, ncol=2) + i_temp<-matrix(nrow = 0, ncol = 2) + + for(i in 1:nrow(matrix)){ + i_temp<-rbind(i_temp, matrix[i,]) + + if(diffs[i] >= diff | i == nrow(matrix)){ + pos<-rbind(pos, i_temp[which(i_temp[,2] == max(i_temp[,2]))[1],]) + i_temp<-matrix(nrow = 0, ncol = 2) + } + } + return(pos) + } + print("Start") + + Apeaks <- getpeaks(obj_ab@traceMatrix[,1]) + Cpeaks <- getpeaks(obj_ab@traceMatrix[,2]) + Gpeaks <- getpeaks(obj_ab@traceMatrix[,3]) + Tpeaks <- getpeaks(obj_ab@traceMatrix[,4]) + + peakthr<-input$thr + peakCusMatrix<-rbind(Gpeaks[Gpeaks[,2] > peakthr,], + Apeaks[Apeaks[,2] > peakthr,], + Tpeaks[Tpeaks[,2] > peakthr,], + Cpeaks[Cpeaks[,2] > peakthr,]) + + peakCusMatrix<-peakCusMatrix %>% as.data.frame() %>% rename(Int=V2) %>% group_by(indexes) %>% summarise(Int=max(Int)[1]) + peakCusMatrix<-as.matrix(peakCusMatrix) + + pos1<-redDiffs(peakCusMatrix, input$dist1) + pos<-redDiffs(pos1, input$dist2) + print(nrow(pos)) + + pos1<-pos1[,1] + pos<-pos[,1] + + print("Peaks Detected") + + primarypeaks <- pos + diffs <- diff(c(0,primarypeaks)) + starts <- primarypeaks - 0.5*diffs + stops <- c(primarypeaks[1:(length(primarypeaks)-1)] + + 0.5*diffs[2:length(diffs)], + primarypeaks[length(diffs)] + 0.5*diffs[length(diffs)] + ) + + primary <- NULL + secondary <- NULL + tempPosMatrix <- matrix(nrow=length(starts), ncol=4) + tempAmpMatrix <- matrix(nrow=length(starts), ncol=4) + + ratio<-input$ratio + + length(starts) + + for(i in 1:length(starts)) { + Apeak <- peakvalues(Apeaks, starts[i], stops[i]) + Cpeak <- peakvalues(Cpeaks, starts[i], stops[i]) + Gpeak <- peakvalues(Gpeaks, starts[i], stops[i]) + Tpeak <- peakvalues(Tpeaks, starts[i], stops[i]) + if(is.na(Apeak[2]) & + is.na(Cpeak[2]) & + is.na(Gpeak[2]) & + is.na(Tpeak[2])) next #rare case where no peak found + signals <- c(Apeak[1], Cpeak[1], Gpeak[1], Tpeak[1]) + tempAmpMatrix[i,] <- signals + positions <- c(Apeak[2], Cpeak[2], Gpeak[2], Tpeak[2]) + tempPosMatrix[i,] <- positions + signalratios <- signals/max(signals, na.rm=TRUE) + Bases <- c("A", "C", "G", "T") + # Bases <- c("G", "A", "T", "G") + Bases[signalratios < ratio] <- NA + #sort by decreasing signal strength + Bases <- Bases[order(signals, decreasing=TRUE)] + positions <- positions[order(signals, decreasing=TRUE)] + if(length(Bases[!is.na(Bases)]) == 4 + | length(Bases[!is.na(Bases)]) == 0) { + primary <- c(primary, "N") + secondary <- c(secondary, "N") + } + else if(length(Bases[!is.na(Bases)]) > 1) { + primary <- c(primary, Bases[1]) + Bases2 <- Bases[2:4] + secondary <- c(secondary, + mergeIUPACLetters(paste(sort(Bases2[!is.na(Bases2)]), + collapse=""))) + } + else { + primary <- c(primary, Bases[1]) + secondary <- c(secondary, Bases[1]) + } + } + + obj_ab@peakPosMatrix <- tempPosMatrix[rowSums(!is.na(tempPosMatrix)) > 0,] + obj_ab@peakAmpMatrix <- tempAmpMatrix[rowSums(!is.na(tempPosMatrix)) > 0,] + obj_ab@primarySeqID <- "sangerseq package primary basecalls" + obj_ab@primarySeq <- DNAString(paste(primary, collapse="")) + obj_ab@secondarySeqID <- "sangerseq package secondary basecalls" + obj_ab@secondarySeq <- DNAString(paste(secondary, collapse="")) + obj_ab<<-obj_ab + obj$seq<-T + print(as.character(obj_ab@primarySeq)) + print("Finished") + } + }) + + output$Primseq <- renderText({ + observeEvent(obj$seq, {}) + + if (!is.null(obj$seq) & input$old != ""){ + print(1) + as.character(obj_ab@primarySeq) + } + }) + + output$align <- renderText({ + observeEvent(obj$seq, {}) + + if (!is.null(obj$seq) & input$old != ""){ + print(2) + alignament<-msaClustalW(c("Old"=toupper(input$old),"New"=obj_ab@primarySeq %>% as.character), type="dna") + paste(capture.output(print(alignament, show="complete")), collapse="\n") + } + }) - # draw the histogram with the specified number of bins - hist(x, breaks = bins, col = 'darkgray', border = 'white') - }) } # Run the application