Canviar threshold total per filtre local i posar un màxim a la 2a agrupació de pics.

chromatoR/app.R

@@ -119,7 +119,7 @@ server <- function(input, output) {
         else return(c(max(region[,2], na.rm=TRUE), region[which.max(region[,2]),1]))
       }
       
-      redDiffs<-function(matrix, diff=5){
+      redDiffs<-function(matrix, diff=5, max=0){
         diffs<-diff(matrix[,1])
         pos<-matrix(nrow=0, ncol=2)
         i_temp<-matrix(nrow = 0, ncol = 2)
@@ -127,7 +127,7 @@ server <- function(input, output) {
         for(i in 1:nrow(matrix)){
           i_temp<-rbind(i_temp, matrix[i,])
           
-          if(diffs[i] >= diff | i == nrow(matrix)){
+          if(diffs[i] >= diff | i == nrow(matrix) | nrow(i_temp) == max){
             pos<-rbind(pos, i_temp[which(i_temp[,2] == max(i_temp[,2]))[1],])
             i_temp<-matrix(nrow = 0, ncol = 2)
           }
@@ -149,10 +149,29 @@ server <- function(input, output) {
       
       # Peaks joining and filtering by defined threshold
       peakthr<-input$thr
-      peakCusMatrix<-rbind(Gpeaks[Gpeaks[,2] > peakthr,], 
+      
-                           Apeaks[Apeaks[,2] > peakthr,], 
+      filt_peaks<-function(Mpeaks, prob=0.75){
-                           Tpeaks[Tpeaks[,2] > peakthr,], 
+        Mpeaks_filt<-matrix(nrow = 0, ncol=2)
-                           Cpeaks[Cpeaks[,2] > peakthr,])
+        means<-matrix(nrow=0, ncol=2)
+        for (i in 1:nrow(Mpeaks)){
+          pack<-Mpeaks[Mpeaks[,1] >= Mpeaks[i,1]-100 & Mpeaks[,1] <= Mpeaks[i,1]+100,]
+          if (Mpeaks[i,2] > quantile(pack[,2], na.rm=T, probs=prob)[[1]]){
+            Mpeaks_filt<-rbind(Mpeaks_filt, Mpeaks[i,])
+          }
+          means<-rbind(means, c(Mpeaks[i,1], quantile(pack[,2], na.rm=T, probs=prob)[[1]]))
+        }
+        return(list(Mpeaks_filt, means))
+      }
+
+      
+      # peakCusMatrix<-rbind(Gpeaks[Gpeaks[,2] > peakthr,],
+      #                      Apeaks[Apeaks[,2] > peakthr,],
+      #                      Tpeaks[Tpeaks[,2] > peakthr,],
+      #                      Cpeaks[Cpeaks[,2] > peakthr,])
+      
+      peakCusMatrix<-filt_peaks(rbind(Gpeaks,Apeaks,Tpeaks,Cpeaks),0.75)
+      peakmeans<<-peakCusMatrix[[2]]
+      peakCusMatrix<-peakCusMatrix[[1]]
 
       # In case there is more than one peak by position, we keep the higher
       peakCusMatrix<-peakCusMatrix %>% as.data.frame() %>% 
@@ -161,7 +180,7 @@ server <- function(input, output) {
 
       # Using the custom function, we aggregate peaks defining two distance thresholds
       pos1<-redDiffs(peakCusMatrix, input$dist1)
-      pos<-redDiffs(pos1, input$dist2)[,1]
+      pos<-redDiffs(pos1, input$dist2, max=2)[,1]
 
       progress$set(message = "Peaks Detected", value = 2)
       print("Peaks Detected")
@@ -341,6 +360,9 @@ server <- function(input, output) {
       
       ranPeaks<-range(picks$rows)
       
+      # pmeans<-peakmeans %>% as.data.frame() %>% rename(index=1, value=2)
+      # pmeans<-pmeans[]
+      
       obj$plot<-obj_ab@traceMatrix %>% as.data.frame() %>% 
         dplyr::rename(A=V1,C=V2,G=V3,T=V4) %>% add_column(rows=1:nrow(.)) %>% 
         filter(rows >= ranPeaks[1] & rows <= ranPeaks[2]) %>% gather(Base, Value, -rows) %>% 
@@ -351,7 +373,10 @@ server <- function(input, output) {
         geom_text(data=picks, aes(label=Base, y=3800, color=factor(Base,levels=c("G","A","T","C"))))+
         geom_text(data=picks, aes(label=num, y=4250), size=3, angle=90, hjust=0.5, vjust=0.5)+
         scale_color_manual(values=c("G"="black", "A"="#66CC00", "T"="red","C"="blue"))+
-        geom_hline(yintercept = input$thr)+
+        # geom_hline(yintercept = input$thr)+
+        geom_point(data=peakmeans %>% as.data.frame() %>% rename(index=1, value=2) %>%
+                     filter(index >= ranPeaks[1] & index <= ranPeaks[2]),
+                   aes(index, value))+
         theme_classic()
       print("Final Imagen")
     }