Threshold es defineix pel 50% dels pics. Es pot modificar a posteriori.

chromatoR/app.R

@@ -14,7 +14,8 @@ ui <- fluidPage(
         tabPanel("Detección de Bases",
           sidebarPanel(
             fileInput("file1", "Sube fichero ab1", multiple = FALSE),
-            numericInput("thr", label = "Peak Threshold", value = 450),
+            # numericInput("thr", label = "Peak Threshold", value = 450),
+            uiOutput('thr'),
             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),
@@ -60,6 +61,31 @@ server <- function(input, output) {
     }
   })
   
+  output$thr<-renderUI({
+    if (!is.null(obj$aborig)){
+      print(1)
+      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])
+      }
+      
+      Apeaks <- getpeaks(obj_ab@traceMatrix[,1])
+      Cpeaks <- getpeaks(obj_ab@traceMatrix[,2])
+      Gpeaks <- getpeaks(obj_ab@traceMatrix[,3])
+      Tpeaks <- getpeaks(obj_ab@traceMatrix[,4])
+      
+      peakCusMatrix<-rbind(Gpeaks, 
+                           Apeaks, 
+                           Tpeaks, 
+                           Cpeaks)
+      thr_calc<-quantile(peakCusMatrix[,2], 0.5)
+      numericInput("thr", label = "Peak Threshold", value = thr_calc)
+    }
+  })
   
   observeEvent(input$calab, {
     if (!is.null(obj$aborig)){
@@ -183,7 +209,12 @@ server <- function(input, output) {
       obj_ab@secondarySeqID <- "sangerseq package secondary basecalls"
       obj_ab@secondarySeq <- DNAString(paste(secondary, collapse=""))
       obj_ab<<-obj_ab
-      obj$seq<-T
+      if (is.null(obj$seq)){
+        obj$seq<-1
+      }else{
+        obj$seq<-obj$seq+1
+      }
+        
       print(as.character(obj_ab@primarySeq))
       progress$set(message = "Peaks Detected", value = 3)
       progress$close()