Primera versió de la app. Mostra seqüència i alineament a una referència.

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",
+          fileInput("file1", "Sube fichero ab1", multiple = FALSE),
-                        "Number of bins:",
+          numericInput("thr", label = "Peak Threshold", value = 450),
-                        min = 1,
+          numericInput("dist1", label = "1st Distance Reduction", value = 4),
-                        max = 50,
+          numericInput("dist2", label = "2nd Distance Reduction", value = 7),
-                        value = 30)
+          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
   
-    output$distPlot <- renderPlot({
+  observe({
-        # generate bins based on input$bins from ui.R
+    if (!is.null(input$file1)){
-        x    <- faithful[, 2]
+      obj$aborig<-readsangerseq(input$file1$datapath)
-        bins <- seq(min(x), max(x), length.out = input$bins + 1)
+    }
+  })
+  
+  
+  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])
+      }
+
+      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