---
title: "Accurate gene consensus at low nanopore coverage"
author:
- name: Rocio Espada
  affiliation: ESPCI Paris
package: single
output:
  BiocStyle::html_document
abstract: |
  This file contains the full code for replicating the analysis on our manuscript.
vignette: |
  %\VignetteIndexEntry{single_fullcode}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
df_print: kable
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, eval=FALSE)
knitr::opts_chunk$set(engine.opts = list(bash = "-l"))
# knitr::opts_knit$set(root.dir ="~/SINGLe/")

```

This code reproduces the analysis made on our manuscript "Accurate gene consensus at low nanopore coverage", from the raw fast5 files returned by the Nanopore Sequencer to final figures.

Here, we applied SINGLe to the gene of KlenTaq, a truncated variant of the well-known Taq polymerase, of approximately 1.7 kb in length. We first tested it on a small set of seven known mutants containing 2 to 9 point mutations (Small Set, or KT7), and later a larger library of approximately 1200 variants (KTLib).

Make sure that in the path were you are running this code there is:

* Folder 1_ReferenceFiles that has:
    - KT7_sequences_aligned.fasta  
    - KTrefamplicon.fasta
    - KTB_digested.fasta           
    - KTrefORF_1662.fasta
* Folder 2_NanoporeRawReads that has:
    - SmallSet_fast5 (a folder. If you have the tar, un tar it)
    - LargeSet_fast5 (a folder. If you have the tar, un tar it)
* Replace the path in the previous chunk by the correct root.dir


The chunks on the general inputs section load packages, functions and variables that are needed for different chunks below. They do not use a lot of memory, so I'd recommend running them all before running the particular section you are interested in. 

Also, sections are not independent between them. The clearest example is the single_train variable that it is obtained in the Single correction for the small library and then it is used for the correction of the large library long after. 

You'll need to have installed the following software:

* minimap2 (https://github.com/lh3/minimap2)
* ont_fast5_api (https://github.com/nanoporetech/ont_fast5_api)
* nanopolish (https://github.com/jts/nanopolish)
* racon (https://github.com/isovic/racon)
* medaka (https://github.com/nanoporetech/medaka)
* samtools (http://www.htslib.org/)

# General inputs

This includes loading needed packages, defining functions and loading data that will be used in some parts of the code. It's fast and light, so just run all chunks before going on.


Packages
```{r, eval=FALSE}
require(single)
suppressPackageStartupMessages(require(plyr))
suppressPackageStartupMessages(require(dplyr))
options(dplyr.summarise.inform=F)
suppressPackageStartupMessages(require(reshape2))
suppressPackageStartupMessages(require(foreach))
suppressPackageStartupMessages(require(ggplot2))
suppressPackageStartupMessages(require(ggbreak))
suppressPackageStartupMessages(require(grid))
suppressPackageStartupMessages(require(gridExtra))
suppressPackageStartupMessages(require(stringr))
suppressPackageStartupMessages(library(Biostrings))
suppressPackageStartupMessages(require(cowplot))
suppressPackageStartupMessages(library(e1071))
```

Functions
```{r, eval=FALSE}
detect_homopolymer_positions <- function(sequence){
  l <- length(sequence)
  homopolymer_positions <- list()
  aux_positions <- 1
  counter <- 0
  for(i in seq_along(sequence)[-l]){
    if (sequence[i+1]==sequence[i]){
      aux_positions <- c(aux_positions,i+1)
    }else{
      counter <- counter+1
      homopolymer_positions[[counter]] <- aux_positions
      aux_positions <- i+1
    }
  }
  if(sequence[i+1]==sequence[i]){
    counter <- counter+1
    homopolymer_positions[[counter]] <- aux_positions
  }
  return(homopolymer_positions)
}
biostr_to_char <- function(x){
    y <- strsplit(as.character(x), split="")[[1]]
    return(y)
}
detect_mutations <- function(sequence,reference){
    ind <- which(sequence != reference)
    muts <- paste0(reference[ind],ind,sequence[ind])
    return(muts)
}
replace_str_in_pos   <- function(string, pos, replacement){
  if(pos>nchar(string)){stop('Replacement intended outsife string')}
  string_vec   <- strsplit(string, split="")[[1]]
  string_vec[pos] <- replacement
  string_out <- paste(string_vec, collapse = "")
  return(string_out)
}
bc_one_mutation    <- function(x){
  bases <- c("A","C","G","T")
  y <- c()
  n <- 0
  for(i in 1:nchar(x)){
    n <- n+1
    non_wt <- setdiff(bases, substr(x, i,i))
    non_wt <- paste0(non_wt, collapse = "")
    non_wt <- paste0("[",non_wt, "]", collapse = "")
    y[n] <- replace_str_in_pos(x, i, non_wt)
  }
  for(i in 2:(nchar(x)-1)){
    n <- n+1
    y[n] <- paste0(substr(x,1,(i-1)),substr(x,(i+1), nchar(x)))
  }
  y <- unique(y)
  return(y)
}
dna_reverse_complement_vector <- function(x){
  y <- rev(dna_complement_vector(x))
  return(y)
}
dna_complement <- function(x){
  if(length(x)>1){stop("dna_complement: x has length > 1")}
  if(x=="A" | x=="a"){ y <- "T"}
  if(x=="C" | x=="c"){ y <- "G"}
  if(x=="G" | x=="g"){ y <- "C"}
  if(x=="T" | x=="t"){ y <- "A"}

  return(y)
}
dna_complement_vector <- function(x){
  y <- vapply(x, dna_complement, USE.NAMES = F)
  return(y)
}
dna_complement_string <- function(x){
  x <- strsplit(x, split="")[[1]]
  y <- dna_complement_vector(x)
  y <- paste(y, collapse = "")
  return(y)
}
dna_reverse_complement_string <- function(x){
  x <- strsplit(x, split="")[[1]]
  y <- dna_reverse_complement_vector(x)
  y <- paste(y, collapse = "")
  return(y)
}
```

Paths
```{r, eval=FALSE}
path_references  <- "1_ReferenceFiles/"
path_raw_data <- "2_NanoporeRawReads/"
path_save_figures <- "6_Figures/"
path_analysis_lib7 <- "3_Analysis/KT7_1700_2100/"
path_analysis_lib <- "3_Analysis/KTLib/"

if(!dir.exists(path_save_figures)){dir.create(path_save_figures)}
if(!dir.exists(path_analysis_lib7)){dir.create(path_analysis_lib7)}
if(!dir.exists(path_analysis_lib)){dir.create(path_analysis_lib)}
```

Colors
```{r, eval=FALSE}
red      <- rgb(218,0,0,maxColorValue = 255)
green    <- rgb(0,122,0,maxColorValue = 255)
red_tr   <- rgb(230,0,0,alpha = 200, maxColorValue = 255)
blue_tr  <- rgb(0,0,218,alpha = 250, maxColorValue = 255)
green_tr <- rgb(0,122,0,alpha = 200, maxColorValue = 255)
blue     <- rgb(0,0,250, maxColorValue = 255)

single_color     <- "#2727fa"  
medaka_color     <- "#26a026"  
guppy_color      <- "#fa2929"
nanopolish_color <- "#fa8c28"  
noweights_color  <- "#8c26fa"
racon_color      <- "#dcdc26"
nextpolish_color <- grey(.6)
colors_vector <- c(single_color,medaka_color,guppy_color,nanopolish_color,noweights_color, racon_color,nextpolish_color)
names(colors_vector) <- c("single","medaka","nanopore","nanopolish","no_weights", "racon","nextpolish")

colors_vector_capital <- colors_vector
names(colors_vector_capital) <- c("SINGLe","Medaka","Guppy","Nanopolish","No weights","Racon","NextPolish")

methods_capital <- setNames(c("Nanopolish","Medaka","Racon","SINGLe","Guppy","No weights","NextPolish"),
                                              c("nanopolish","medaka","racon","single","nanopore","no_weights","nextpolish"))
```

Plotting themes
```{r, eval=FALSE}
theme_plot <- theme_bw()+
    theme(line=element_line(colour = "black", 
                            size = 0.5, linetype = "dashed",lineend = "square"),
          rect= element_rect(fill = NULL, colour = NULL, 
                             size = 1,linetype = "solid",inherit.blank = FALSE),
          text=element_text(size = 10),title = element_text(size=10),
          axis.title = element_text(size=rel(1)),
          axis.text.x=element_text(angle=0, size=rel(.8)),
          axis.text.y=element_text(angle=0, size=rel(.8)),
          plot.tag = element_text(size=rel(1), face = "bold", vjust=-1),
          panel.grid = element_blank(), panel.border = element_rect(size=1),
          legend.background = element_blank(), 
          plot.margin=margin(t=0,r=3,l=0,b=0, unit="pt"),
          strip.background = element_rect(fill="white"),
          legend.text = element_text(size=rel(0.8))
    )

theme_set(theme_plot)

```

Reference sequence
```{r, eval=FALSE}
reference_file <-  "1_ReferenceFiles/KTrefORF_1662.fasta"
wildtypye_bstr <- readDNAStringSet(reference_file)
wildtype <- toupper(biostr_to_char(wildtypye_bstr))

wt_homopolymers <- detect_homopolymer_positions(wildtype)
wt_homopolymers <- wt_homopolymers[vapply(wt_homopolymers,length)>1]
pos_wt_homopolymers <- unlist(wt_homopolymers)


```

Experimental data for set of KlenTaq's 7 variants
```{r, eval=FALSE}
kt7_barcodes <- paste0("barcode",c("05","06","07","08","09","10","11"))
kt7_true_mutants <- data.frame(
                        Barcode = kt7_barcodes,
                        true_consensus = c("C867G G1120A C1214T G1316A",
                            "G23A C245T", 
                            "G94C A378G A905G G1023T A1381T",
                            "T425A G846A C851T G1403A", 
                            "G307A G490- T659A T675A G993A A1554G",
                            "A41T G331T G349A C772T C879T C1474A C1475T G1530-",
                            "G303A G376A C517G T857A T1056A G1550A"))

kt7_bc2mut <- c("#1 - 2sub","#2 - 4sub","#3 - 4sub","#4 - 5sub","#5 - 6sub","#6 - 7sub 1del","#7 - 5sub 1del")
names(kt7_bc2mut) <- paste0("barcode",c("06","05","08","07","11","10","09"))


mutant_sequences   <- readDNAStringSet("1_ReferenceFiles/KT7_sequences_aligned.fasta") #known sequences of mutants (different from reference_sequence)
kt_true_mutations <- list()
for(i in 1:length(mutant_sequences)){
    mutant_seq <- toupper(stringr::str_split(as.character(mutant_sequences[[i]]), pattern = "")[[1]])
    index <- which (mutant_seq!= wildtype)
    kt_true_mutations[[i]] <- data.frame(sequence_id=rep(names(mutant_sequences)[i], length(index)),
                                      position = index,  nucleotide =mutant_seq[index])
}
kt_true_mutations <- do.call(rbind,kt_true_mutations)

```


Others
```{r, eval=FALSE}
subsets       <- c(seq(3,9,by=1),seq(10,50,by=5))
n_repetitions <- 50
```


# Wild type and small set of 7 variants of KlenTaq

## Experiment

In a unique sequencing run we used a barcoding kit to measure several samples at the same time. The table below indicates the number of the barcode (that will be returned by guppy_barcoder) and it's content

| BarcodeID | VariantID | Mutations                                              |
|-----------|-----------|--------------------------------------------------------|
| BC01      |           | wildtype                                               |
| BC05      | Variant 2 | C867G G1120A C1214T C1316A                             |
| BC06      | Variant 1 | G23A C245T                                             |
| BC07      | Variant 4 | G94C A378G A905G G1023T A1381T                         |
| BC08      | Variant 3 | T425A G846A C851T G1403A                               |
| BC09      | Variant 7 | G307A T659A T675A G993A A1554G G490                    |
| BC10      | Variant 6 | A41T G331T G349A C772T A785G C879T C1474A C1475T G1530 |
| BC11      | Variant 5 | G303A G376A C517G T857A T1056A G1550A                  |



## Basecalling, demultiplexing and aligning to reference

Basecalling was done with the following command lines in bash (replace paths if needed)

```{bash, basecalling, eval= F}
## BASECALL USING GUPPY
guppy_basecaller -i 2_NanoporeRawReads/SmallSet_fast5/ -s 3_Analysis/SmallSet_SUP/ -c dna_r9.4.1_450bps_sup.cfg -x 'cuda:0' -r 

#Merge all files in a unique file
cat 3_Analysis/SmallSet_SUP/pass/*fastq > 3_Analysis/KT7.fastq
# for example: cat SmallSet_SUP/pass/*fastq > SmallSet.fastq

#Filter sequences by length (1700 to 2100 bp)
mkdir 3_Analysis/KT7_1700_2100
awk 'NR % 4 ==2 && length($0) > 1700 && length($0) < 2100 {print f  "\n" $0;getline;print;getline;print}  {f=$0}' 3_Analysis/KT7.fastq > 3_Analysis/KT7_1700_2100/KT7_1700_2100.fastq

#Demultiplexing by Guppy
guppy_barcoder -i 3_Analysis/KT7_1700_2100/ -s 3_Analysis/KT7_1700_2100

#Merge files to one file per barcode
cat 3_Analysis/KT7_1700_2100/barcode01/*fastq >> 3_Analysis/KT7_1700_2100/barcode01.fastq
cat 3_Analysis/KT7_1700_2100/barcode05/*fastq >> 3_Analysis/KT7_1700_2100/barcode05.fastq
cat 3_Analysis/KT7_1700_2100/barcode06/*fastq >> 3_Analysis/KT7_1700_2100/barcode06.fastq
cat 3_Analysis/KT7_1700_2100/barcode07/*fastq >> 3_Analysis/KT7_1700_2100/barcode07.fastq
cat 3_Analysis/KT7_1700_2100/barcode08/*fastq >> 3_Analysis/KT7_1700_2100/barcode08.fastq
cat 3_Analysis/KT7_1700_2100/barcode09/*fastq >> 3_Analysis/KT7_1700_2100/barcode09.fastq
cat 3_Analysis/KT7_1700_2100/barcode10/*fastq >> 3_Analysis/KT7_1700_2100/barcode10.fastq
cat 3_Analysis/KT7_1700_2100/barcode11/*fastq >> 3_Analysis/KT7_1700_2100/barcode11.fastq

```

Align using minimap2 and create count files.

```{bash, minimap2, eval= F}
#!/bin/bash
minimap2 -ax map-ont --sam-hit-only  1_ReferenceFiles/KTrefamplicon.fasta 3_Analysis/KT7_1700_2100/barcode01.fastq > 3_Analysis/KT7_1700_2100/barcode01.sam
samtools view -S -b 3_Analysis/KT7_1700_2100/barcode01.sam > 3_Analysis/KT7_1700_2100/barcode01.bam #TRANSFORM TO BAM
samtools sort 3_Analysis/KT7_1700_2100/barcode01.bam -o 3_Analysis/KT7_1700_2100/barcode01.sorted.bam #SORT BAM FILE

for (( counter=5; counter<12; counter++ ))
    do
if [ $counter -lt 10 ];
then
bc=0$counter
else
    bc=$counter
fi
name=barcode$bc
minimap2 -ax map-ont --sam-hit-only  1_ReferenceFiles/KTrefamplicon.fasta 3_Analysis/KT7_1700_2100/$name.fastq >3_Analysis/KT7_1700_2100/$name.sam
samtools view -S -b 3_Analysis/KT7_1700_2100/$name.sam > 3_Analysis/KT7_1700_2100/$name.bam # TRANSFORM TO BAM
samtools sort 3_Analysis/KT7_1700_2100/$name.bam -o 3_Analysis/KT7_1700_2100/$name.sorted.bam # SORT BAM FILE

echo -n "$bc "
done
printf "\n"
```

## SINGLe correction

```{r, eval=FALSE}
pos_start <- 60
pos_end <- 1721
```

```{r, warning=F, single, eval= F}
tic <- proc.time()

## TRAIN
#Forward
single_train <- single_train(bamfile= "3_Analysis/KT7_1700_2100/barcode01.sorted.bam",
               output = "3_Analysis/KT7_1700_2100/barcode01_5d4" ,
               mean.n.mutations = 5.4,verbose = T,
               rates.matrix = mutation_rate,
               refseq_fasta = reference_file,
               pos_start = pos_start, pos_end = pos_end, 
               save_partial=T, save_final=T)

## EVALUATE
cl <- parallel::makeForkCluster(7)
doParallel::registerDoParallel(cl)
foreach(n = 5:11)%dopar%{
    if(n<10){x=paste0("0",n)}else{x=n}
    single_evaluate(bamfile = paste0("3_Analysis/KT7_1700_2100/barcode",x,".sorted.bam"),
                     single_fits   = single_train,
                     output_file = paste0("3_Analysis/KT7_1700_2100/barcode",x,"_SINGLe.fastq"),
                     refseq_fasta= reference_file,
                     pos_start=pos_start,pos_end=pos_end,
                     verbose=T,gaps_weights="minimum",
                     save_original_scores = T, 
                     save=T)
}

toc <- proc.time()

cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)
```

------------------------------------------------------------------------

## Analysis of error rate
Compute error rate on wild type reads
```{r, eval= F}
tic <- proc.time()
reads_wt <- read.table(paste0(path_analysis_lib7,"barcode01_5d4_SINGLe_countsPNQ.txt"), header=T)
reads_wt <- reads_wt %>%
    mutate(isWT= nucleotide==wildtype[pos]) %>%
    dplyr::rename(Position=pos, Nucleotide=nucleotide, counts=count) 

total_error_rate <- reads_wt %>% group_by(isWT) %>% dplyr::summarise(counts=sum(counts))
total_error_perc <- total_error_rate[1,2]/sum(total_error_rate$counts)*100
cat("Error rate:", as.numeric(round(total_error_perc,2)),"%")

## Qscore distribution
reads_wt_qscore_per_position <- reads_wt %>%
    group_by(QUAL, isWT) %>%
    dplyr::summarise(Counts=sum(counts))

## Errors per position
reads_wt_errors_by_position <- reads_wt %>%
    group_by(Position, Nucleotide,isWT) %>%
    dplyr::summarise(Counts=sum(counts)) %>%
    mutate(n_plot =  floor(Position/300))
toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)
```

### Figure 1A
```{r, fig 1A, fig.width=3, fig.height=3, eval= F}
Figure_1A <- ggplot(reads_wt_qscore_per_position,aes(x=QUAL,y=Counts/(10^5),col=isWT))+
    geom_line(lwd=0.5)+
    xlab(expression(Q[score]))+ylab(bquote('Counts'~(x10^5)))+
    scale_color_manual(breaks=c(TRUE,FALSE), labels=c("Correct","Error"),values=c(green_tr,red_tr))+
    theme_plot+
    theme(legend.direction="vertical", 
          legend.position = c(.7,.8),
          legend.spacing.y = unit(0.05,"line"),
          legend.key.height = unit(0.5,"line"),
          legend.key.width = unit(0.6,"line"),
          legend.text=element_text(size=rel(.8)), 
          legend.title=element_blank(), 
          legend.background = element_blank(), 
          plot.margin = margin(t=0,r = 0,b = 0,l=0))+
    labs(tag = "A")

Figure_1A 
```

### Figure S1 (equivalent to Fig1A normalized)
```{r, fig S1, fig.width=3, fig.height=3 , eval= F}
Figure_1A_norm <- ggplot(reads_wt_qscore_per_position %>% dplyr::group_by(isWT) %>% dplyr::mutate(NormCounts=Counts/sum(Counts)),aes(x=QUAL,y=NormCounts,col=isWT))+
    geom_line(lwd=0.5)+
    xlab(expression(Q[score]))+ylab('Normalized counts')+
    scale_color_manual(breaks=c(TRUE,FALSE), labels=c("Correct","Error"),values=c(green_tr,red_tr))+
    theme_plot+
    theme(legend.direction="vertical", 
          legend.position = c(.7,.7),
          legend.spacing.y = unit(0.05,"line"),
          legend.key.height = unit(1,"line"),
          legend.text=element_text(size=rel(1)), 
          legend.title=element_blank(), 
          legend.background = element_blank())

Figure_1A_norm
ggsave(Figure_1A_norm, filename = paste0(path_save_figures,"Figure1Asuppl.png"),width =8 ,height = 8,dpi='print',units = "cm")
```

### Figure 1B
```{r, fig 1B, fig.width=3, fig.height=3, eval= F}
Figure_1B <- ggplot(reads_wt_errors_by_position %>% 
                        filter(!isWT & Nucleotide !="-" & Position >100 & Position < 150),
                    aes(x=Position,y=Counts,fill=Nucleotide))+
    guides(fill=guide_legend(ncol=2))+
    geom_col()+
    scale_x_continuous(breaks=seq(100,150,by=25))+
    theme_plot+
    theme(legend.direction="horizontal",
          legend.position = c(.3,.85),
          legend.key.size =  unit(.5,units = "line"),
          legend.box = "vertical" ,
          legend.background = element_blank(),
          legend.text=element_text(size=rel(.8)),
          legend.title=element_blank(),
          plot.margin = margin(t=0,r = 0,b = 0,l=0)
    )+
    labs(tag = "B")
Figure_1B
```

### Figure S2 (extended Fig 1B)
```{r, figS2, fig.width=7, fig.height=10, eval= F}
Figure_S2 <- ggplot(reads_wt_errors_by_position %>% filter(!isWT & Nucleotide !="-"),
                                     aes(x=Position, y=Counts,fill=Nucleotide))+
    geom_col()+facet_wrap(~n_plot, ncol=1, scales="free_x")+
    theme_plot+
    theme(strip.background = element_blank(), strip.text.x = element_blank(),
          legend.direction="horizontal", legend.position = "top",
          rect= element_rect(fill = NULL, colour = NULL, size = 1,linetype = "solid",inherit.blank = FALSE)
    )

Figure_S2
ggsave(Figure_S2, file=paste0(path_save_figures,"FigureS2.png"), dpi = 'print', width = 16, height=22,units="cm")
```

Choose example of fits  
```{r, fits examples, eval= F}
tic <- proc.time()
## Example of fit
data_fits <- read.table("3_Analysis/KT7_1700_2100/barcode01_5d4_SINGLe_fits.txt", header = T)  
data_mut <- read.table("3_Analysis/KT7_1700_2100/barcode01_5d4_SINGLe_data.txt", header = T)   

position <-  547             #   Choose position to plot, 547 in manuscript
mut.base <- "A"         #   Choose nucleotide to plot, A in manuscript
str <- "+"

#Filter all data by position and base
data_chosen_pos_base <- data_mut %>%
    filter(pos==position & nucleotide==mut.base & strand==str) %>%
    filter(count>0 | counts.wt>0) %>%
    mutate(ratio=counts.wt/(count+counts.wt),
           ratio_scaled =counts.wt.scaled/(counts.wt.scaled+counts.scaled))
wt  <- as.character(data_chosen_pos_base$wt.base[1])

#Filter all fits' data by position and base
data_fits_chosen_pos_base <- data_fits %>% 
    filter(pos==position & (nucleotide==mut.base | nucleotide==wt) & strand==str)
class_fit_prior <- vapply(data_fits_chosen_pos_base$priors, class)
if(any(class_fit_prior=="numeric")){
    ind_out <- which(class_fit_prior=="numeric");
    data_fits_chosen_pos_base <- data_fits_chosen_pos_base[-ind_out,]
}

#Construct fitting curves
quality_range <- 0:35   #* You'll draw the fitted curves for this Qscores (x-values)

curve_fitted_corrected <- data.frame(quality=quality_range,
                                     value=glm.predict.(quality_range, slope = data_fits_chosen_pos_base$prior_slope,
                                                        intercept =data_fits_chosen_pos_base$prior_intercept))

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

### Figure 1D
```{r,Fig1D, fig.width=3, fig.height=3, eval= F}
Figure_1D <- ggplot(data_chosen_pos_base)+
    geom_point(aes(x=QUAL, y=ratio_scaled), size=.5)+
    geom_line(data=curve_fitted_corrected, aes(x=quality, y=value), linetype="dashed", color=single_color)+
    scale_y_continuous(breaks=c(0,0.5,1), limits = c(0,1))+scale_x_continuous(limits=c(0,35))+
    xlab(expression(Q[score]))+ylab(expression(N[correct]))+ggtitle("SINGLe")+
    theme_plot+
    theme(plot.title = element_text(size=8))+
    labs(tag = "D")


Figure_1D
```

### Figure 1C 
```{r, Fig1C, fig.width=3, fig.height=3, eval= F}
fit_nopriors <- stats::glm(data_chosen_pos_base$ratio~ data_chosen_pos_base$QUAL,family = "quasibinomial")
fit_nopriors_coeff       <- stats::coefficients(fit_nopriors)
curve_fitted_corrected <- data.frame(quality=quality_range,
                                     value=glm.predict.(quality_range, slope = fit_nopriors_coeff[2],
                                                        intercept =fit_nopriors_coeff[1]))

Figure_1C <-  ggplot(data_chosen_pos_base, aes(x=QUAL, y=ratio))+geom_point()+
    geom_line(data=curve_fitted_corrected, aes(x=quality, y=value), linetype="dashed", color=guppy_color)+
    scale_y_continuous(breaks=c(0,0.5,1),limits = c(0,1))+scale_x_continuous(limits=c(0,35))+
    xlab(expression(Q[score]))+ylab(expression(N[correct]))+ggtitle("Guppy")+
    theme_plot+
    theme(plot.title = element_text(size=8), plot.margin = margin(t=0.5, l=0.5))+
    labs(tag = "C")
Figure_1C
```
    

### Figure 1 composition
```{r, Figure 1, eval= F}
Figure1 <- ( Figure_1A | Figure_1B ) / (Figure_1C | Figure_1D) 
Figure1
ggsave(Figure1, file=paste0(path_save_figures,"Figure1.png"),
       dpi = 'print', width = 8.5, height=8.5,units="cm")

```



## Strand bias analysis
```{r, strand bias, eval= F}
tic <- proc.time()
reads_wt <- read.table(paste0(path_analysis_lib7,"barcode01_5d4_SINGLe_countsPNQ.txt"), header=T) %>%
    dplyr::rename(Strand=strand) %>%
    mutate(isWT = nucleotide==wildtype[pos])

# mean Qscore is similar
tot_counts_for <- sum(reads_wt %>% filter(Strand=="+") %>% select(count))
mean_qscore_forward <- reads_wt %>% filter(Strand=="+") %>% summarise( sum(count*QUAL)/tot_counts_for)

tot_counts_rev <-sum(reads_wt %>% filter(Strand=="-") %>% select(count))
mean_qscore_reverse <- reads_wt %>% filter(Strand=="-") %>% summarise( sum(count*QUAL)/tot_counts_rev)

# Qscore distributions are similar:
counts_hist <- reads_wt %>%
    dplyr::group_by(QUAL,Strand) %>%
    dplyr::summarise(counts=sum(count)) %>%
    ungroup()%>%
    mutate(Proportion = counts/sum(counts)) %>%
    dplyr::rename(Qscore=QUAL) %>%
    mutate(Strand =ifelse(Strand=="+", "Forward","Reverse"))

# Percentage of errors are similar
perc_errors <- reads_wt %>%
    group_by(Strand, isWT) %>%
    summarise(counts = sum(count)) %>%
    mutate(isWT = c("yes","no")[2-as.numeric(isWT)])%>%
    reshape2::dcast(formula = Strand ~  isWT) %>%
    mutate(perc_error = no/(yes+no)*100)

# Percentages of errors per position are different
perc_errors_perpos <- reads_wt %>%
    group_by(Strand, isWT,pos) %>%
    summarise(counts = sum(count)) %>%
    mutate(isWT = c("yes","no")[2-as.numeric(isWT)])%>%
    reshape2::dcast(formula = Strand +pos ~  isWT) %>%
    mutate(perc_error = no/(yes+no)*100)

perc_errors_perpos_cast <- reshape2::dcast(perc_errors_perpos,formula=pos~Strand)%>%
    dplyr::rename(Forward="+", Reverse="-")

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)
```

### Figure S11
```{r, figS11, eval= F}
FigS11A <- ggplot(counts_hist, aes(x=Qscore, y=Proportion,linetype=Strand))+
    geom_line()+theme(legend.position = c(.7,.7))
FigS11B <- ggplot(perc_errors_perpos_cast, aes(x=Forward,y=Reverse))+
    geom_point(col=rgb(0,0,0,0.3))+
    scale_x_log10(limits=c(00.1,100))+
    scale_y_log10(limits=c(00.1,100))+
    xlab("Error (%) - forward strand")+
    ylab("Error (%) - reverse strand")


# Thus fits are different
all_fits <- read.table(paste0(path_analysis_lib7,"/barcode01_5d4_SINGLe_fits.txt"), header = T)
forward_fits <- all_fits %>% filter(strand=="+")
reverse_fits <- all_fits %>% filter(strand=="-")
qval <- seq(0,40,by=0.1)
par(mfrow=c(1,4), mar=c(4,4,1,1),mgp=c(2,.7,0), las=1)
first <- T
for(n in c(1278:1281)+4*7-1){
  y_for <- glm.predict.(x=qval, slope = forward_fits$prior_slope[n], intercept = forward_fits$prior_intercept[n])
  y_rev <- glm.predict.(x=qval, slope = reverse_fits$prior_slope[n], intercept = reverse_fits$prior_intercept[n])
  name <- paste(forward_fits$pos[n], forward_fits$nucleotide[n])
  if(first){
    df <- rbind(data.frame(Qscore=qval,Sense="Forward",Fit=y_for,n=name),
                data.frame(Qscore=qval,Sense="Reverse",Fit=y_rev,n=name))
    first <- F
  }else{
      df <- rbind(df,
                  rbind(data.frame(Qscore=qval,Sense="Forward",Fit=y_for,n=name),
                        data.frame(Qscore=qval,Sense="Reverse",Fit=y_rev,n=name))
                  )}
}
plot_fits_examples <- ggplot(df, aes(x=Qscore,y=Fit,lty=Sense))+geom_line()+facet_wrap(~n,nrow=1)

FigS11up <- plot_grid(FigS11A,FigS11B+theme(plot.margin=margin(l=10,r=5.5,t=7,b=5.5)),labels=c("A","B"),hjust = 0.1,vjust=1)
FigS11C <- plot_grid(plot_fits_examples,labels="C",hjust = 0.1)
FigS11 <-
    plot_grid(FigS11up, FigS11C, nrow=2,hjust = -1)
ggsave(paste0(path_save_figures,"FigS11_strandsense.png"),
       width = 7.62,height = 5.68,units = "in",dpi = 300)
FigS11
```

------------------------------------------------------------------------

## Compare signal to noise with p_right cut off

For Guppy and SINGLe, compare how does the signal and noise change if we accept as truth only the nucleotides with a p_right higher than a cut off.


Computation (takes \~ 1min). 
```{r inputs I, eval= F}
tic <- proc.time()
# Identify errors and correct reads (signal and noise)
n <- 0
for(bc_i in kt7_barcodes){
    message('Proccessing ', bc_i, '\n')
    n<- n+1
    #Barcode info
    actual_mutations <- kt_true_mutations%>%                   #mutations present in this barcode
        filter(sequence_id==bc_i)%>% select(nucleotide, position)
    actual_mutations <- actual_mutations[which(actual_mutations$nucleotide!="-"),]

    #Load corrected data
    seqs_single <- readQualityScaledDNAStringSet(paste0(path_analysis_lib7, bc_i,"_SINGLe.fastq"))
    aux_seqs <- vapply(as.character(seqs_single),strsplit, split="")
    aux_quals <- 1-as(quality(seqs_single),"NumericList")
    aux_pos <- lapply(aux_seqs, seq_along)
    aux_names <- c(vapply(names(seqs_single),rep, each=1662))
    data_single <- data.frame(nucleotide = unlist(aux_seqs),
                                p_SINGLe=unlist(aux_quals),
                                position=unlist(aux_pos) , 
                            seq_id=aux_names)
    rownames(data_single)<-NULL

    #Load original data
    bf <- Rsamtools::BamFile(paste0(path_analysis_lib7, bc_i,".bam"))
    reads <- Rsamtools::scanBam(bf)
    #keep sequences that start at least at pos_start
    index <- which(reads[[1]]$pos<=pos_start)
    reads[[1]] <- vapply(reads[[1]], function(x){x[index]})

    seqs_guppy <- GenomicAlignments::sequenceLayer(reads[[1]]$seq,
                                reads[[1]]$cigar,to = "reference")
    names(seqs_guppy) <- reads[[1]]$qname
    scores_guppy <- GenomicAlignments::sequenceLayer(reads[[1]]$qual,
                                reads[[1]]$cigar,to = "reference")
    names(scores_guppy) <- reads[[1]]$qname

    #Fill with gaps at the end of sequences shorter than pos_end
    index_short_sequences <- which(width(seqs_guppy)<pos_end)
    for(i in index_short_sequences){
        seqs_guppy[i] <- paste0(as.character(seqs_guppy[i]), 
                                   paste0(rep("-", pos_end-width(seqs_guppy[i])),
                                          collapse = ""),
                                   collapse = "")
        scores_guppy[i] <- paste0(as.character(scores_guppy[i]),
                                    paste0(rep("-", pos_end-width(scores_guppy[i])),
                                           collapse = ""),
                                    collapse = "")
    }
    seqs_guppy <- subseq(seqs_guppy, start=pos_start+1-reads[[1]]$pos,end=pos_end+1-reads[[1]]$pos)
    scores_guppy <- subseq(scores_guppy, start=pos_start+1-reads[[1]]$pos,end=pos_end+1-reads[[1]]$pos)
   
    aux_seqs <- vapply(as.character(seqs_guppy),strsplit, split="")
    aux_quals <- 1-as(scores_guppy,"NumericList")
    aux_pos <- lapply(aux_seqs, seq_along)
    aux_names <- c(vapply(names(seqs_guppy),rep, each=1662))
    data_guppy <- data.frame(seq_id=aux_names,
                            position=unlist(aux_pos), 
                            nucleotide = unlist(aux_seqs),
                            p_Guppy=unlist(aux_quals)
                            )
    rownames(data_guppy)<-NULL

    data <- full_join(data_single,data_guppy, by=c("seq_id","nucleotide","position")) %>%
        mutate(isWT = nucleotide==wildtype[position]) %>%
        filter(isWT==0 & nucleotide!="-") %>%
        group_by(position, nucleotide, p_Guppy,p_SINGLe) %>%
        dplyr::summarise(counts=n())
    
    signal_data_aux <- left_join(actual_mutations, data,by = c("nucleotide", "position"))
    noise_data_aux  <- anti_join(data, actual_mutations,by = c("nucleotide", "position"))

    if(nrow(signal_data_aux)+nrow(noise_data_aux) != nrow(data)){stop('Check joins of data')}

    #Store data
    if(n==1){
        signal_data <- signal_data_aux
        noise_data  <- noise_data_aux
    }else{
        signal_data <- rbind(signal_data,signal_data_aux)
        noise_data  <- rbind(noise_data, noise_data_aux)
    }
    rm(signal_data_aux, noise_data_aux, seqs_single, seqs_guppy, data, aux_seqs, aux_quals, aux_pos)
}


#Compute signal and noise for different cut-offs
p_cutoff_vec <- sort(c(0,10^seq(-10,-1, by=.5),seq(0.1,.99, by=0.01)))

signal_tp <- data.frame(matrix(NA, ncol=5, nrow=length(p_cutoff_vec)))
colnames(signal_tp) <- c("pcutoff", "Guppy", "SINGLe", "Guppy - weighted", "SINGLe - weighted")
signal_fn <- noise_tn <- noise_fp <- signal_tp

noise_data <- noise_data%>% ungroup()
signal_data <- signal_data %>% ungroup()
for(i in seq_along(p_cutoff_vec)){
    p_cutoff <- p_cutoff_vec[i]

    signal_tp_aux <-  signal_data %>%
        dplyr::summarise(counts_guppy = sum(counts[p_Guppy       >= p_cutoff]),
                         counts_single = sum(counts[p_SINGLe >= p_cutoff]),
                         counts_guppy_w = sum( (counts*p_Guppy)       [p_Guppy       >= p_cutoff]),
                         counts_single_w = sum( (counts*p_SINGLe) [p_SINGLe >= p_cutoff])
        )

    signal_fn_aux <-  signal_data %>%
        dplyr::summarise(counts_guppy = sum(counts[p_Guppy       < p_cutoff]),
                         counts_single = sum(counts[p_SINGLe < p_cutoff]),
                         counts_guppy_w = sum( (counts*p_Guppy)       [p_Guppy       < p_cutoff]),
                         counts_single_w = sum( (counts*p_SINGLe) [p_SINGLe < p_cutoff])
        )

    noise_fp_aux <-  noise_data %>%
        dplyr::summarise(counts_guppy = sum(counts[p_Guppy       >= p_cutoff]),
                         counts_single = sum(counts[p_SINGLe >= p_cutoff]),
                         counts_guppy_w = sum( (counts*p_Guppy)       [p_Guppy       >= p_cutoff]),
                         counts_single_w = sum( (counts*p_SINGLe) [p_SINGLe >= p_cutoff])
        )

    noise_tn_aux <-  noise_data %>%
        dplyr::summarise(counts_guppy = sum(counts[p_Guppy       < p_cutoff]),
                         counts_single = sum(counts[p_SINGLe < p_cutoff]),
                         counts_gupy_w = sum( (counts*p_Guppy)       [p_Guppy       < p_cutoff]),
                         counts_single_w = sum( (counts*p_SINGLe) [p_SINGLe < p_cutoff])
        )

    signal_tp[i,] <- c(p_cutoff, unlist(signal_tp_aux))
    noise_fp [i,] <- c(p_cutoff, unlist(noise_fp_aux))
    signal_fn[i,] <- c(p_cutoff, unlist(signal_fn_aux))
    noise_tn [i,] <- c(p_cutoff, unlist(noise_tn_aux))
}


aux_tp <- melt(signal_tp, id.vars = "pcutoff", value.name = "signal_tp")
aux_fp <- melt(noise_fp, id.vars = "pcutoff", value.name = "noise_fp")
aux_fn <- melt(signal_fn, id.vars = "pcutoff", value.name = "signal_fn")
aux_tn <- melt(noise_tn, id.vars = "pcutoff", value.name = "noise_tn")
rates <- aux_tp %>% full_join(aux_fp,by = c("pcutoff", "variable"))%>% full_join(aux_fn,by = c("pcutoff", "variable"))%>%full_join(aux_tn,by = c("pcutoff", "variable"))
rates <- rates %>%
    mutate(ratio = signal_tp / noise_fp)%>%
    mutate(tpr   = signal_tp / (signal_tp + signal_fn))%>%
    mutate(fpr   = noise_fp /  (noise_fp + noise_tn))

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```


### Figure 2A and 2B

```{r Fig2A 2B, fig.width=3, fig.height=3, eval= F}
Figure_2A <- ggplot(rates %>% filter(variable=="Guppy" | variable=="SINGLe"), aes(y=tpr,x=fpr, col=variable))+
    geom_line()+
    scale_color_manual(values=c(guppy_color,single_color), breaks = c("Guppy","SINGLe"))+
    scale_x_continuous(breaks=c(0,0.5,1), limits = c(0,1))+xlab("FPR")+
    scale_y_continuous(breaks=c(0,0.5,1), limits = c(0,1))+ylab("TPR")+
    theme_plot+
    theme(legend.title = element_blank(), legend.position = c(0.6,0.2),legend.spacing.y = unit(0,"cm"),
          legend.key.size = unit(0.5,"lines"), plot.margin = unit(c(0,5.5,5.5,5.5),"points"))+
    labs(tag = "A")

Figure_2B <-  ggplot(rates %>% filter(variable=="Guppy - weighted" | variable=="SINGLe - weighted") %>% mutate(variable = ifelse(variable == "Guppy - weighted", "Guppy", "SINGLe")), aes(pcutoff, y=ratio, col=variable))+
    geom_line()+
    scale_y_continuous(breaks=seq(1,4))+ylab("Signal / Noise")+
    scale_x_continuous(breaks=c(0,0.5,1), limits = c(0,1))+xlab(expression("Cut off on p"["right"]))+
    scale_color_manual(values=c(guppy_color,single_color), breaks=c("Guppy", "SINGLe"), name="")+
    theme_plot+theme(legend.position = c(0.3,0.85),legend.spacing.y = unit(0,"cm"),
                                legend.key.size = unit(0.5,"lines"), plot.margin = unit(c(0,5.5,5.5,5.5),"points"))+
    labs(tag = "B")

Figure_2A
Figure_2B
```

### Figure S5 (similar to Fig 2B)

```{r, FigS5, fig.width=3, fig.height=3, eval= F}
Figure_Sup3 <- ggplot(rates %>% filter(variable=="Guppy" | variable=="SINGLe"), aes(pcutoff, y=ratio, col=variable))+
    geom_line()+
    scale_y_continuous(breaks=seq(1,4))+ylab("Signal / Noise")+
    scale_x_continuous(breaks=c(0,0.5,1), limits = c(0,1))+xlab(expression("Cut off on p"["right"]))+
    scale_color_manual(values=c(guppy_color,single_color), breaks=c("Guppy", "SINGLe"), name="")+
    theme_plot+theme(legend.position = c(0.3,0.85),legend.spacing.y = unit(0,"cm"),
                                legend.key.size = unit(0.5,"lines"), plot.margin = unit(c(0,5.5,5.5,5.5),"points"))
Figure_Sup3
```

------------------------------------------------------------------------

## Convergence of consensus
Compute convergence of consensus. It takes long. 

General inputs
```{r, eval= F}
wildtype_homopolymers <-
    c(F, wildtype[2:length(wildtype)] == wildtype[1:(length(wildtype)-1)]) |
    c(wildtype[1:(length(wildtype)-1)] == wildtype[2:length(wildtype)] ,F)
pos_wt_homopolymers <- which(wildtype_homopolymers)
```



* Consensus by SINGLe, Guppy and No Weights, not modifying homopolymers
```{r, cons subsets single, eval= F}
tic <- proc.time()
for(i in 1:7){

   #Load corrected data
    seqs_single <- readQualityScaledDNAStringSet(paste0(path_analysis_lib7, kt7_barcodes[i],"_SINGLe.fastq"))
    aux_seqs <- vapply(as.character(seqs_single),strsplit, split="")
    aux_quals <- 1-as(quality(seqs_single),"NumericList")
    aux_pos <- lapply(aux_seqs, seq_along)
    aux_names <- c(vapply(names(seqs_single),rep, each=1662))
    data_single <- data.frame(nucleotide = unlist(aux_seqs),
                                p_SINGLe=unlist(aux_quals),
                                position=unlist(aux_pos) , 
                            SeqID=aux_names)
    rownames(data_single)<-NULL
    
    seqs_original <- readQualityScaledDNAStringSet(paste0(path_analysis_lib7, kt7_barcodes[i],"_SINGLe_original.fastq"))
    aux_seqs <- vapply(as.character(seqs_original),strsplit, split="")
    aux_quals <- 1-as(quality(seqs_original),"NumericList")
    aux_pos <- lapply(aux_seqs, seq_along)
    aux_names <- c(vapply(names(seqs_original),rep, each=1662))
    data_guppy <- data.frame(SeqID=aux_names,
                            position=unlist(aux_pos), 
                            nucleotide = unlist(aux_seqs),
                            p_Guppy=unlist(aux_quals)
                            )
    rownames(data_guppy) <- NULL

    data <- full_join(data_single,data_guppy, by=c("SeqID","nucleotide","position")) %>%
        mutate(isWT = nucleotide==wildtype[position]) %>%
        mutate(p_SINGLe=ifelse(is.na(p_SINGLe), 1, p_SINGLe))%>%
        mutate(p_noweights=1)

    #Compute consensus for each barcode using all sequences available and using the three available methods
    file_out_consensus<- paste0(path_analysis_lib7,kt7_barcodes[i],"_ConsensusBySINGLE.txt")
    if(file.exists(file_out_consensus)){file.remove(file_out_consensus)}
    cat("Barcode\tnseqs\trepetition\tmutation\n", file=file_out_consensus,append=F)
    pb <- txtProgressBar(0,length(subsets)*n_repetitions,style=3)
    counter <- 0
    bc_seqsID <- unique(data$SeqID)
    for( s in subsets){
        for (r in 1:n_repetitions) {
            subset_bc_reads <- data %>% 
                filter(SeqID %in% sample(bc_seqsID, s))

            cons_single <- weighted_consensus(df = subset_bc_reads %>%
                                                  select(nucleotide,p_SINGLe,position), 
                                           cutoff_prob = 0)
            cons_guppy  <- weighted_consensus(df = subset_bc_reads %>% 
                                                  select(nucleotide,p_Guppy,position),
                                           cutoff_prob = 0)
            cons_noweight <- weighted_consensus(df = subset_bc_reads %>%
                                                    select(nucleotide,p_noweights,position),
                                            cutoff_prob = 0)
            muts_single   <- detect_mutations(biostr_to_char(cons_single), wildtype)
            muts_guppy    <- detect_mutations(biostr_to_char(cons_guppy), wildtype)
            muts_noweight <- detect_mutations(biostr_to_char(cons_noweight), wildtype)

            n_single <- length(muts_single)
            n_guppy <- length(muts_guppy)
            n_noweight <- length(muts_noweight)
            n_tot <- n_single+n_guppy+n_noweight
            df <- data.frame(barcode = rep(kt7_barcodes[i],n_tot),
                            nseqs=s,repetition=r,
                            mutation=c(muts_single,muts_guppy,muts_noweight),
                            method = c(rep("SINGLe", n_single),
                                       rep("Guppy",n_guppy), 
                                       rep("No weights",n_noweight)))
            
            write.table(df, append = T, col.names = F,row.names = F, file = file_out_consensus)
            counter <- counter+1
            setTxtProgressBar(pb,counter)
        }
    }

}
rm(data)

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

* Consensus by SINGLe, Guppy and No Weights, sorting homopolymers
```{r,cons subsets single homopol, eval= F}
tic <- proc.time()
for(i in 1:7){

  #Load corrected data
    seqs_single <- readQualityScaledDNAStringSet(paste0(path_analysis_lib7, kt7_barcodes[i],"_SINGLe.fastq"))
    aux_seqs <- vapply(as.character(seqs_single),strsplit, split="")
    aux_quals <- 1-as(quality(seqs_single),"NumericList")
    aux_pos <- lapply(aux_seqs, seq_along)
    aux_names <- c(vapply(names(seqs_single),rep, each=1662))
    data_single <- data.frame(nucleotide = unlist(aux_seqs),
                                p_SINGLe=unlist(aux_quals),
                                position=unlist(aux_pos) , 
                            SeqID=aux_names)
    rownames(data_single)<-NULL
    
    seqs_original <- readQualityScaledDNAStringSet(paste0(path_analysis_lib7, kt7_barcodes[i],"_SINGLe_original.fastq"))
    aux_seqs <- vapply(as.character(seqs_original),strsplit, split="")
    aux_quals <- 1-as(quality(seqs_original),"NumericList")
    aux_pos <- lapply(aux_seqs, seq_along)
    aux_names <- c(vapply(names(seqs_original),rep, each=1662))
    data_guppy <- data.frame(SeqID=aux_names,
                            position=unlist(aux_pos), 
                            nucleotide = unlist(aux_seqs),
                            p_Guppy=unlist(aux_quals)
                            )
    rownames(data_guppy) <- NULL

    data <- full_join(data_single,data_guppy, by=c("SeqID","nucleotide","position")) %>%
        mutate(isWT = nucleotide==wildtype[position]) %>%
        mutate(p_SINGLe=ifelse(is.na(p_SINGLe), 1, p_SINGLe))%>%
        mutate(p_noweights=1)
    


    #Compute consensus for each barcode using all sequences available and using the three available methods
    file_out_consensus<- paste0(path_analysis_lib7,kt7_barcodes[i],"_ConsensusBySINGLE_rmHomopolforVCS.txt")
    if(file.exists(file_out_consensus)){file.remove(file_out_consensus)}
    cat("Barcode\tnseqs\trepetition\tmutation\n", file=file_out_consensus,append=F)
    pb <- txtProgressBar(0,length(subsets)*n_repetitions,style=3)

    counter <- 0
    bc_seqsID <- unique(data$SeqID)
    for( s in subsets){
        for (r in 1:n_repetitions) {
            subset_bc_reads <- data %>% 
                filter(SeqID %in% sample(bc_seqsID, s))

            index_homopolymers <- which(subset_bc_reads$position %in% pos_wt_homopolymers & subset_bc_reads$nucleotide=="-")
            n_memory <- c()
            for(n in index_homopolymers){
                if(n %in% n_memory){next()}
                hp_pos <- subset_bc_reads$position[n]
                hp_pos_ref <- wt_homopolymers[vapply(wt_homopolymers, function(x){hp_pos %in% x})][[1]]
                hp_vec_logical <- hp_pos_ref==hp_pos
                hp_vec_length <- length(hp_vec_logical)
                hp_vec <- rep(NA,hp_vec_length)
                ind_aux <- which(hp_vec_logical)
        
                hp_vec[1:ind_aux] <- rev(n+1-1:ind_aux)
                if(ind_aux<hp_vec_length){
                    hp_vec[(ind_aux+1):hp_vec_length] <- n+c(1:(hp_vec_length-ind_aux))
                }
                n_memory <- hp_vec
        
                data_aux <- subset_bc_reads[hp_vec,]
                if(all(data_aux$nucleotide=="-")){next()}
        
                pos <- data_aux$position
                nucl <- data_aux$nucleotide
                ind_order <- c(which(nucl!="-"),which(nucl=="-"))
        
                subset_bc_reads$position[hp_vec] <- pos[order(ind_order)]
        
                if(length( which(subset_bc_reads$position %in% pos_wt_homopolymers &
                                 subset_bc_reads$nucleotide=="-"))<length(index_homopolymers)){
                    cat("n=",n,"\n")
                    stop()
                }
            }
            subset_bc_reads <- subset_bc_reads %>% arrange(SeqID,position)
    
            cons_single <- weighted_consensus(df = subset_bc_reads %>%
                                                  select(nucleotide,p_SINGLe,position), 
                                           cutoff_prob = 0)
            cons_guppy  <- weighted_consensus(df = subset_bc_reads %>% 
                                                  select(nucleotide,p_Guppy,position),
                                           cutoff_prob = 0)
            cons_noweight <- weighted_consensus(df = subset_bc_reads %>%
                                                    select(nucleotide,p_noweights,position),
                                            cutoff_prob = 0)
            muts_single   <- detect_mutations(biostr_to_char(cons_single), wildtype)
            muts_guppy    <- detect_mutations(biostr_to_char(cons_guppy), wildtype)
            muts_noweight <- detect_mutations(biostr_to_char(cons_noweight), wildtype)

            n_single <- length(muts_single)
            n_guppy <- length(muts_guppy)
            n_noweight <- length(muts_noweight)
            n_tot <- n_single+n_guppy+n_noweight
          if(n_tot==0){
               df <- data.frame(barcode = rep(kt7_barcodes[i],3),
                            nseqs=s,
                            repetition=r,
                            mutation="wildtype",
                            method = c("SINGLe", "Guppy","No weights"))
       
          }else{
            df <- data.frame(barcode = rep(kt7_barcodes[i],n_tot),
                            nseqs=s,
                            repetition=r,
                            mutation=c(muts_single,muts_guppy,muts_noweight),
                            method = c(rep("SINGLe", n_single),rep("Guppy",n_guppy), rep("No weights",n_noweight)))
          }
            write.table(df, append = T, col.names = F,row.names = F, file = file_out_consensus)
            counter<- counter+1
            setTxtProgressBar(pb,counter)
        }
    }

}

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

* Consensus by Nanopolish and Racon+Medaka:

Prepare fast5 files for nanopolish. Please indicate proper paths indicated in capital letters
```{r, fast5 for nanopolish, eval= F}
tic <- proc.time()
PATH_TO_FAST5_FILES<- "2_NanoporeRawReads/SmallSet_fast5/" #PATH_TO_FAST5_FILES location of nanopore reads
dir.create(paste0(path_analysis_lib7,"fast5files/")) # individual fast5 files will be stored there temporally 

# Split to individual files
split_fast5 <- paste0("multi_to_single_fast5 -i ",PATH_TO_FAST5_FILES, " -s ",path_analysis_lib7, "fast5files/")  
system2(split_fast5)

input_table_barcodes <- read.table(paste0(path_analysis_lib7,"barcoding_summary.txt"),header=T)
input_table_barcodes <- input_table_barcodes %>%
    filter(barcode_arrangement %in% kt7_barcodes) %>%
    select(read_id,barcode_arrangement)


# Split individual files into folders according to barcode
individual_files <- dir(pattern="*.fast5",path=paste0(path_analysis_lib7,"fast5files/"),recursive=T)
individual_files <- data.frame(file =individual_files)
individual_files$read_id <- vapply(individual_files$file, function(x) { 
    sub(pattern = ".fast5",replacement = "",x=strsplit(x,split="/", fixed=T)[[1]][2])
    })
individual_files <- left_join(individual_files,input_table_barcodes)

for(bc in kt7_barcodes){
    dir.create(paste0(path_analysis_lib7,"fast5files/", bc))
    barcode_files <- individual_files %>% filter(barcode_arrangement==bc)
    
    commands <- paste0("mv ",path_analysis_lib7,"fast5files/", barcode_files$file," ", path_analysis_lib7,"fast5files/",bc,"/")
    vapply(commands,system2)
}

# Create one file per barcode
dir.create(paste0(path_analysis_lib7, "fast5files_grouped"))
for(bc in kt7_barcodes){
    dir.create(paste0(path_analysis_lib7, "fast5files_grouped/",bc))
    command <- paste0("single_to_multi_fast5 -i ",path_analysis_lib7,"fast5files/", bc, " -s ", path_analysis_lib7,"fast5files_grouped/",bc, " -f ", bc, " -n ", 1000)
    system2(command)
}

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

Create bash file to compute consensus by nanopolish and medaka
```{r, eval= F}
tic <- proc.time()
dir.create(paste0(path_analysis_lib7,"subsets/"))
wd <- getwd()
for(bc in kt7_barcodes){
    cat(bc,"\n")

    store_path <- paste0(path_analysis_lib7,"subsets/")
    fastq_file <- readLines(paste0(path_analysis_lib7, bc,".fastq"))
    name_lines <- grep(pattern="^@",x=fastq_file)
    bash_file <- paste0(path_analysis_lib7,bc, ".sh")
    cat("#!/bin/bash \n",file=bash_file, append=F)
    cat(paste0("cd ", wd,"\n"), file=bash_file,append=T)

    results_racon      <- paste0(path_analysis_lib7,bc,"_ConsensusbyRacon.txt")
    results_medaka     <- paste0(path_analysis_lib7,bc,"_ConsensusbyMedaka.txt")
    nanopolish_results <- paste0(path_analysis_lib7,bc,"_ConsensusbyNanopolish.txt")
    file.create(results_racon,overwrite=TRUE)
    file.create(results_medaka,overwrite=TRUE)
    file.create(nanopolish_results,overwrite=TRUE)
    
    for(ns in subsets){
        cat("\t subset size",ns,"\n")
        if(ns> length(name_lines)){next()}
        for(r in 1:n_repetitions){
            # cat("----------------",ns, r, "\n")
            basename <- paste0(bc,"_",ns,"_",r)
            filename <- paste0(store_path,basename)
            fastq <- paste0(filename,".fastq")

            ## Make file of subset of sequences
            choose_lines <- sample(name_lines, size= ns)
            choose_lines_all <- sort(c(choose_lines,choose_lines+1,choose_lines+2,choose_lines+3))
            writeLines(text = fastq_file[choose_lines_all], con = file(fastq))

            sam        <- paste0(filename,".sam")
            sorted     <- paste0(filename,".sorted.fasta")
            racon      <- paste0(filename,".racon.fasta")
            medaka     <- paste0(filename,"_medaka.fasta")
            nanopolish <- paste0(filename,"_nanopolish.txt")

            ## Run minimap2
            line_minimap <- paste0("minimap2 -ax map-ont ", reference_file, " ", fastq," > ", sam)

            ## Run racon
            line_racon <- paste0("racon ",fastq," ",sam, " ",reference_file, " >> ", racon)
            line_results_racon <- paste0("echo '>",basename, "' >> ",results_racon)
            line_results_racon2 <- paste0("tail -n 1 ",racon, " >> ",results_racon)
            line_remove_racon <- paste0("rm ", racon,".fai ",racon, ".mmi ",racon  )

            ## Run Medaka
            line_medaka <- paste0("medaka_consensus -i ",fastq," -d ",racon,"  -o ", medaka," -t 4 ")
            line_results_medaka <- paste0("echo '>",basename, "' >> ",results_medaka)
            line_results_medaka2 <- paste0("tail -n 1 ",filename,"_medaka.fasta/consensus.fasta", " >> ",results_medaka)
            line_remove_medaka <- paste0("rm -r ",filename,"_medaka.fasta/")

            ## Run nanopolish
            line_nanopolish_index  <- paste0("nanopolish index -d ", path_analysis_lib7, "fast5files_grouped/",bc," ",fastq)
            line_nanpolish_sort    <- paste0(" samtools sort ",sam," -o ",filename,".sorted.fastq -T temporal.tmp ")
            line_nanopolish_index2 <- paste0("samtools index ", filename,".sorted.fastq")
            line_nanopolish_variants <- paste0("nanopolish variants --consensus -o ", nanopolish, " -r ",fastq," -b ",filename,".sorted.fastq -g ", reference_file)
            line_results_nanopolish <- paste0("awk '!/^#/ {print \"",filename,"\\t\" $0}' ",nanopolish, " >> ",nanopolish_results)
            line_remove_nanopolish <- paste0("rm ", filename,".fastq.index* ", filename,".sorted* ", nanopolish, " ", filename,".sam ")

            cat(line_minimap,
                line_racon,
                line_results_racon,
                line_results_racon2,
                line_medaka,
                line_results_medaka,
                line_results_medaka2,
                line_nanopolish_index,
                line_nanpolish_sort,
                line_nanopolish_index2,
                line_nanopolish_variants,
                line_results_nanopolish,
                line_remove_racon,
                line_remove_medaka,
                line_remove_nanopolish,
                file=bash_file, sep="\n", append=T)
        }
    }
}

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

Run barcodeXX.sh bash files in a terminal with conda activate medaka:

```{bash, eval= F}
cd 3_Analysis/KT7_1700_2100
chmod +x *sh

conda activate medaka

nohup ./barcode05.sh > bc05_out.log 2>bc05_err.log &
nohup ./barcode06.sh > bc06_out.log 2>bc06_err.log &
nohup ./barcode07.sh > bc07_out.log 2>bc07_err.log &
nohup ./barcode08.sh > bc08_out.log 2>bc08_err.log &
nohup ./barcode09.sh > bc09_out.log 2>bc09_err.log &
nohup ./barcode10.sh > bc10_out.log 2>bc10_err.log &
nohup ./barcode11.sh > bc11_out.log 2>bc11_err.log &

```

Run NextPolish

```{r, nextpolish, eval= F}
tic <- proc.time()
create_run_cfg <- function(run_cfg_file, ref_fasta, workdir){
    cat("[General]\n", file=run_cfg_file, append=FALSE)
    cat('job_type = local\n', file=run_cfg_file, append=TRUE)
    cat('job_prefix = nextPolish\n', file=run_cfg_file, append=TRUE)
    cat('task = best\n', file=run_cfg_file, append=TRUE)
    cat('rewrite = yes\n', file=run_cfg_file, append=TRUE)
    cat('rerun = 3\n', file=run_cfg_file, append=TRUE)
    cat('parallel_jobs = 2\n', file=run_cfg_file, append=TRUE)
    cat('multithread_jobs = 4\n', file=run_cfg_file, append=TRUE)
    cat('genome =', ref_fasta,' #genome file\n', file=run_cfg_file, append=TRUE)
    cat('genome_size = auto\n', file=run_cfg_file, append=TRUE)
    cat('workdir =',workdir, "\n",file=run_cfg_file, append=TRUE, sep="")
    cat('   polish_options = -p {multithread_jobs}\n', file=run_cfg_file, append=TRUE)

    cat('[lgs_option]\n', file=run_cfg_file, append=TRUE)
    cat('lgs_fofn =lgs.fofn\n', file=run_cfg_file, append=TRUE)
    cat('lgs_options = -min_read_len 1k -max_depth 100\n', file=run_cfg_file, append=TRUE)
    cat('lgs_minimap2_options = -x map-ont\n', file=run_cfg_file, append=TRUE)
    return(0)
}

dir.create(paste0(path_analysis_lib7,"nextpolish/"))
system2(paste0("cp ", reference_file, " ",path_analysis_lib7,"nextpolish/"))
for(bc in kt7_barcodes){
    save_file <- paste0(path_analysis_lib7,bc, "_ConsensusbyNextPolish.fasta")

    for(ns in subsets){
        cat("\t",ns,"\n")
        if(ns> length(name_lines)){next()}
        for(r in 1:n_repetitions){
            cat("----------------",ns, r, "\n")

            #Create file with subset of sequences
            subset_name <- paste0(bc,"_",ns, "_",r)
            filename_subset <- paste0(run_folder,subset_name,".fastq")
           
            ## Run nextpolish
            system2(paste0("mv ", path_analysis_lib7,"subsets/", subset_name, ".fastq ", path_analysis_lib7,"nextpolish/"))
            system2(paste0("echo '", subset_name,".fastq' >", path_analysis_lib7,"nextpolish/lgs.fofn"))
            create_run_cfg(run_cfg_file=paste0(path_analysis_lib7,"nextpolish/run.cfg"),
                           ref_fasta="KTrefORF_1662.fasta",
                           workdir=paste0(subset_name,"/"))
            system2(paste0("cd ",path_analysis_lib7,"nextpolish/ ; nextPolish run.cfg"))
            system2(paste0("sed -i '1s/.*/>",subset_name,"'/ ",path_analysis_lib7,"nextpolish/",subset_name, "/genome.nextpolish.fasta "))
            system2(paste0("cat ",path_analysis_lib7,"nextpolish/", subset_name, "/genome.nextpolish.fasta >>", save_file))
            system2(paste0("rm -r ", path_analysis_lib7,"nextpolish/", subset_name))
            system2(paste0("rm -r ", path_analysis_lib7,"nextpolish/", subset_name,".fastq"))
            system2(paste0("rm -r ", path_analysis_lib7,"nextpolish/lgs.fofn"))
            system2(paste0("rm -r ", path_analysis_lib7,"nextpolish/run.cfg"))
        }
    }
}
toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

Identify mutations
```{r,identify mutations, eval= F}
tic <- proc.time()
df_nextpolish <- data.frame(
    Barcode=NA,
    Nreads=NA,
    repetition=NA,
    mutation=NA
)
for(bc in kt7_barcodes){
    np<- readDNAStringSet(paste0(path_analysis_lib7,bc, "_ConsensusbyNextPolish.fasta"))
    
    ali <- pairwiseAlignment(pattern = np, subject = wildtypye_bstr)
    for(i in seq_along(ali)){
        ref <- biostr_to_char(subject(ali[i]))
        seq <- biostr_to_char(pattern(ali[i]))
        ind_noI <- which(ref!="-")
        ref <- ref[ind_noI]
        seq <- seq[ind_noI]
        mutations <- detect_mutations(seq,ref)
        ids <- str_split(names(np)[i], "_")[[1]]
        if(length(mutations)==0){next()}
        df_aux <- data.frame(Barcode=ids[1], 
                             Nreads=as.numeric(str_remove(ids[2],"s")), 
                             repetition=as.numeric(str_remove(ids[3],"r")), 
                             mutation = mutations)
        df_nextpolish <- rbind(df_nextpolish,df_aux)
    }
}
df_nextpolish <- df_nextpolish %>% filter(!is.na(Barcode))
toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```


Re-shape results to a unique matrix
```{r, eval= F}
tic <- proc.time()
# Racon
racon_mismatches <- data.frame(Barcode=NA,Nreads=NA,repetition=NA,Position=NA,base.original=NA,base.mutated=NA,mutation=NA)
for(bc in kt7_barcodes){
    racon_consensus <- readDNAStringSet(paste0(path_analysis_lib7,bc,"_ConsensusbyRacon.txt"))
    for(i in seq_along(racon_consensus)){
        name <- names(racon_consensus)[i]
        info <- strsplit(split="_",x=name)[[1]]

        racon <- racon_consensus[[i]]
        ali <- pairwiseAlignment(wildtype_bstr,racon)
        ref <- strsplit(as.character(pattern(ali)), split="")[[1]]
        read <- strsplit(as.character(subject(ali)),split="")[[1]]
        index <- which(ref!=read)
        mismatches <- paste0(ref[index],index,read[index])
        if(length(mismatches)>0){
            df_aux <- data.frame(Barcode=bc,Nreads=info[2],repetition=info[3],Position=index,
                                 base.original=ref[index],base.mutated=read[index],mutation=mismatches)
            racon_mismatches <- racon_mismatches %>% rbind(df_aux)
        }
    }
}
racon_mismatches <- racon_mismatches %>% filter(!is.na(Barcode)) %>% mutate(method="Racon",homopolymers=NA)

# Medaka
medaka_mismatches <- data.frame(Barcode=NA,Nreads=NA,repetition=NA,Position=NA,base.original=NA,base.mutated=NA,mutation=NA)
for(bc in kt7_barcodes){
    medaka_consensus <- readDNAStringSet(paste0(path_analysis_lib7,bc,"_ConsensusbyMedaka.txt"))
    for(i in seq_along(medaka_consensus)){
        name <- names(medaka_consensus)[i]
        info <- strsplit(split="_",x=name)[[1]]

        medaka <- medaka_consensus[[i]]
        ali <- pairwiseAlignment(wildtype_bstr,medaka)
        ref <- strsplit(as.character(pattern(ali)), split="")[[1]]
        read <- strsplit(as.character(subject(ali)),split="")[[1]]
        index <- which(ref!=read)
        mismatches <- paste0(ref[index],index,read[index])
        if(length(mismatches)>0){
            df_aux <- data.frame(Barcode=bc,Nreads=info[2],repetition=info[3],Position=index,
                                 base.original=ref[index],base.mutated=read[index],mutation=mismatches)
            medaka_mismatches <- medaka_mismatches %>% rbind(df_aux)
        }
    }
}
medaka_mismatches <- medaka_mismatches %>% filter(!is.na(Barcode))  %>% mutate(method="Medaka",homopolymers=NA)

# Nanopolish
nanopolish_mismatches <- data.frame(Barcode=NA,Nreads=NA,repetition=NA,Position=NA,base.original=NA,base.mutated=NA)
for(bc in kt7_barcodes){

    results_nanopolish <- system2(paste0("awk \'/^[^#]/ \' ",bc,"_ConsensusbyNanopolish.txt"), intern = T)
    if(length(results_nanopolish)==0){next()}
    aux_nanopolish <- do.call(rbind,vapply(results_nanopolish, strsplit, split="\t"))
    colnames(aux_nanopolish) <- rownames(aux_nanopolish) <- NULL
    aux_nanopolish <- aux_nanopolish[,c(1,3,5,6)]
    info <- vapply(aux_nanopolish[,1],strsplit,split="_")

    if(nrow(aux_nanopolish)>0){
        aux_nanopolish <- data.frame(Barcode=bc,
                                     Nreads=vapply(info, function(x){x[2]}),
                                     repetition=vapply(info, function(x){x[3]}),
                                     Position=aux_nanopolish[,2],
                                     base.original=aux_nanopolish[,3],base.mutated=aux_nanopolish[,4])
        nanopolish_mismatches <- nanopolish_mismatches %>% rbind(aux_nanopolish)
    }
}

nanopolish_mismatches <- nanopolish_mismatches %>% filter(!is.na(Barcode)) %>%
    mutate(Nreads=as.numeric(Nreads), repetition=as.numeric(repetition), Position=as.numeric(Position))

#remove insertions
ind_insertions <- nchar(nanopolish_mismatches$base.mutated)>1
nanopolish_mismatches$base.mutated[ind_insertions] <- vapply(nanopolish_mismatches$base.mutated[ind_insertions] , substr,1,1)
nanopolish_mismatches <- nanopolish_mismatches %>% filter(base.mutated!=base.original)

#split deletions
index_aux <- which(nchar(nanopolish_mismatches$base.original)==2)
nanopolish_mismatches$base.original[index_aux] <- substr(nanopolish_mismatches$base.original[index_aux],2,2)
nanopolish_mismatches$base.mutated[index_aux] <- "-"
nanopolish_mismatches$Position[index_aux] <- as.numeric(nanopolish_mismatches$Position[index_aux])+1

index_aux_3 <- which(nchar(nanopolish_mismatches$base.original)==3)
aux_df_1 <- nanopolish_mismatches[index_aux_3,]
aux_df_1$Position <- as.numeric(aux_df_1$Position)+1
aux_df_1$base.original <- substr(aux_df_1$base.original, 2,2)
aux_df_1$base.mutated <- "-"

aux_df_2 <- nanopolish_mismatches[index_aux_3,]
aux_df_2$Position <- as.numeric(aux_df_2$Position)+2
aux_df_2$base.original <- substr(aux_df_2$base.original, 3,3)
aux_df_2$base.mutated <- "-"

nanopolish_mismatches <- nanopolish_mismatches[-index_aux_3,]
nanopolish_mismatches <- nanopolish_mismatches %>% rbind(aux_df_1) %>% rbind(aux_df_2) %>% arrange(Barcode,Nreads,repetition,Position) %>% mutate(method="Nanopolish",homopolymers=NA, mutation=paste0(base.original,Position,base.mutated))

#SINGLe

single_mismatches <- lapply(kt7_barcodes,
                              function(x){
                                  name=paste0(x,"_ConsensusBySINGLE.txt")
                                  y = read.table(name, header=T,row.names=NULL)
                                  return(y)})
colnames(single_mismatches) <- c("Barcode","Nreads","repetition","mutation","method")
single_mismatches$homopolymers <- "original"


single_mismatches_sortHP <- lapply(kt7_barcodes,
                                   function(x){
                                       name=paste0(x,"_ConsensusBySINGLE_50reps_rmHomopolforVCS.txt")
                                       y = read.table(name, header=T,row.names=NULL)
                                       return(y)})
colnames(single_mismatches_sortHP) <- c("Barcode","Nreads","repetition","mutation","method")
single_mismatches_sortHP$homopolymers = "sorted_at_VCS"

single_mismatches <- rbind(single_mismatches_sortHP,single_mismatches) %>%
    mutate(base.original = str_sub(mutation,1,1),
           base.mutated=str_sub(mutation,-1,-1),
           position=str_sub(mutation,start=2,end=-2)) %>%
    mutate(position = as.numeric(position)) %>% dplyr::rename(Position=position)  %>%
    mutate(method = recode(method,!!!setNames(c("SINGLe","Guppy","No weights"),c("single","nanopore","no_weights"))))


df_nextpolish  <- df_nextpolish %>%
    mutate(method="NextPolish") %>%
    mutate(homopolymers=NA) %>%
    mutate(base.original = str_sub(mutation, 1,1)) %>%
    mutate(base.mutated = str_sub(mutation, -1,-1)) %>%
    mutate(Position = str_sub(mutation,2,-2)) 

#Save
kt7_mismatches_subsets_allMethods <- rbind(medaka_mismatches,nanopolish_mismatches,single_mismatches, df_nextpolish)
write.table(kt7_mismatches_subsets_allMethods, file=paste0(path_analysis_lib7, "KT7_ConsensusConvergence.txt"))

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

* Consensus by NextPolish

Load results of previous chunks
```{r, eval= F}
kt7_mismatches_subsets_allMethods <- read.table(paste0(path_analysis_lib7, "KT7_ConsensusConvergence.txt"), header=T)
```

Compute averages
```{r, eval= F}
tic <- proc.time()
consensus_by_method <- kt7_mismatches_subsets_allMethods %>%
    mutate(variant = recode(Barcode, !!!kt7_bc2mut))

#correct equivalent deletion
ind_aux <- which(consensus_by_method$mutation=="G489-")
consensus_by_method$mutation[ind_aux] <- "G490-"
consensus_by_method$Position[ind_aux] <- 490

consensusStrand_by_method <- consensus_by_method  %>%
    group_by(Barcode,variant,Nreads,repetition,method,homopolymers) %>%
    summarise(consensus = paste(mutation, collapse=" ")) %>%
    left_join(kt7_true_mutants, by="Barcode") %>%
    mutate(correct = consensus==true_consensus)

n_correct_consensus_average <- consensusStrand_by_method %>%
    group_by(Barcode,variant,Nreads,method,homopolymers) %>%
    summarise(total_correct = sum(correct)) %>%
    mutate(success_rate = total_correct/50)


nmismatches_by_method <- consensus_by_method %>%
    group_by(Barcode,variant,Nreads,repetition, method,homopolymers) %>%
    summarise(n_mismatches=n())

average_nmismatches_by_method <- nmismatches_by_method %>%
    group_by(Barcode,variant,Nreads,method,homopolymers) %>%
    summarise(mean=sum(n_mismatches)/50,median=median(n_mismatches))

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)
```

### Figure 2C

```{r, fig.width=3, fig.height=3, eval= F}

linetype_vector <- c("solid","solid","dashed","dashed","solid","solid")
names(linetype_vector) <- c("Medaka","Nanopolish","Guppy","No weights","SINGLe","NextPolish")

average_mismatch_to_plot <- average_nmismatches_by_method %>% 
    filter(method !="Racon") %>%
    filter((method=="SINGLe" & homopolymers=="sorted_at_VCS") |
               is.na(homopolymers) | 
               (method!="SINGLe" & homopolymers=="original") )

n_correct_to_plot <- n_correct_consensus_average %>% filter(method !="Racon") %>%
    filter((method=="SINGLe" & homopolymers=="sorted_at_VCS") |is.na(homopolymers) | (method!="SINGLe" & homopolymers=="original") )

Figure_2C <- ggplot(average_mismatch_to_plot %>%filter(Barcode=="barcode09"),
                           aes(x=Nreads, y=mean, col=method,lty=method))+
    geom_point(size=0.5)+geom_line()+
    ylab("Mismatches")+
    xlab("Reads")+xlim(c(0,20))+
    scale_color_manual(values=colors_vector_capital)+
    scale_linetype_manual(values=linetype_vector, guide="none")+
    theme_plot+
    theme(legend.position = "none",plot.margin = unit(c(0,5.5,5.5,5.5),"points"))+
    labs(tag = "C")

Figure_2C
```

### Figure 2D
```{r, eval= F}
Figure_2D <- ggplot(n_correct_to_plot ,
                     aes(x=Nreads,y=success_rate,col=method, lty=method))+
    geom_point(size=1)+
    geom_line()+
    guides(col=guide_legend(ncol=1))+
    geom_hline(aes(yintercept=0.9),col=grey(.5),linetype="dashed")+
    scale_color_manual(values=colors_vector_capital[-6])+
    scale_linetype_manual(values=linetype_vector, guide="none")+
    theme_plot+
    scale_y_continuous(breaks=c(0,.5,1), limits = c(0,1))+
    scale_x_continuous(breaks=c(10,30,50))+
    theme(
          legend.position = c(0.8,0.5),
          legend.box.margin = margin(0,0,0,0),
          legend.box.spacing=unit(1,"pt"))+
    facet_wrap(~variant, ncol=3, dir = "v")+
    labs(tag = "D",x="Reads",y="Success rate",col="Method")

Figure_2D
```

### Figure 2 composition

```{r, fig.width=7, fig.height=10, eval= F}
Figure2 <- grid.arrange(Figure_2A , Figure_2B ,Figure_2C, Figure_2D, 
             layout_matrix=matrix(c(1,2,3,4,4,4), byrow = T,nrow=2), 
             heights=c(1,3))
Figure2
ggsave(Figure2, file=paste0(path_save_figures,"Figure2.png"), dpi = 'print', width = 16, height=20,units="cm")

```

## Sensibility and sensitivity of consensus by method

Compute TP/FP/TN/FN by method
```{r, eval= F}
tic <- proc.time()
kt7_mismatches_subsets_allMethods <- read.table(paste0(path_analysis_lib7, "KT7_ConsensusConvergence.txt"), header=T)

consensus_by_method <- kt7_mismatches_subsets_allMethods %>%
    mutate(variant = recode(Barcode, !!!kt7_bc2mut))

data_mutations_local <- kt_true_mutations %>%
    mutate(wt = wildtype[position]) %>%
    mutate(mutation =paste0(wt,position,nucleotide)) %>%
    dplyr::rename(Barcode=sequence_id)%>%
    select(Barcode, mutation) %>%
    mutate(true_mut =1) %>%
    group_by(Barcode) %>%
    mutate(n_muts = n())

a <- consensus_by_method %>%
    filter(method!="Racon") %>%
    filter((method=="SINGLe" & homopolymers=="sorted_at_VCS") |is.na(homopolymers) | (method!="SINGLe" & homopolymers=="original") )%>%
    select(Barcode, Nreads, repetition, mutation, method) %>%
    left_join(data_mutations_local %>% select(-n_muts), by=c("Barcode", "mutation")) %>%
    mutate(true_mut = ifelse(is.na(true_mut), 0,true_mut)) %>%
    left_join(data_mutations_local %>% select(-true_mut, - mutation) %>% unique(), by=c("Barcode"))%>%
    mutate(variant = recode(Barcode, !!!kt7_bc2mut)) %>%
    select(-Barcode) %>%
    dplyr::rename(Method=method)

class_vec <- c("True mutations", "False mutations", "False wild type", "True wild type", "Detected mutations")
names(class_vec) <- c("true_pos","false_pos","false_neg","true_neg", "total_mut")

sum_a_melt_av <- a %>%
    group_by(variant,Nreads, repetition,Method) %>%
    summarise(total_mut = n(),
              true_pos = sum(true_mut),
              false_pos =sum(1-true_mut),
              n_muts=mean(n_muts)) %>%
    mutate(false_neg = n_muts-true_pos, true_neg=1662-total_mut-false_neg) %>%
    reshape2::melt(c("variant","Nreads","repetition","Method")) %>%
    group_by(variant,Nreads, Method, variable) %>%
    dplyr::summarise(mean_value = mean(value))%>%
    mutate(variable = recode(variable, !!!class_vec))

sum_sens <- sum_a_melt_av %>%
    group_by(Nreads, variant, Method) %>%
    reshape2::dcast(Nreads+ variant+Method~variable) %>%
    mutate(sensitivity = `True mutations`/(`True mutations` + `False wild type`) ) %>%
    mutate(specificity = `True wild type`/(`True wild type` + `False mutations`) )

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

### Figure S6
```{r, fig.width=7, fig.height=10, eval= F}
Figure_S6A <- ggplot(sum_a_melt_av %>% filter(Nreads <10 & variable != "n_muts" ),
       aes(x=Nreads, y=mean_value, col=Method))+
    geom_line()+
    facet_grid(variable~variant, scales = "free_y", labeller = label_wrap_gen(width=10))+
    xlab("Reads")+ylab("Mean of counts")+
    theme_plot+
    scale_color_manual(values=colors_vector_capital[-6])+
    theme(legend.position = "none")+
    ggtitle(label="A")

Figure_S6B <-ggplot(sum_sens %>% filter(Nreads <10  ),
       aes(x=specificity, y=sensitivity, col=Method))+
    facet_grid(~variant,  labeller = label_wrap_gen(width=10))+
    geom_point()+geom_line()+
    theme_plot+
    scale_color_manual(values=colors_vector_capital[-6])+
    theme(legend.position = "bottom", strip.background = element_blank(),
          strip.text.x = element_blank(),plot.margin = margin(t = 0,r=35,b = 0,l = 10))+
    ggtitle(label="B")+
    scale_x_continuous(breaks=c(0.990,1))+
    scale_y_continuous(breaks=c(0.2,0.6,1.0))


Figure_S6 <- grid.arrange(Figure_S6A , Figure_S6B,
                        layout_matrix=matrix(c(1,2), byrow = T,nrow=2),
                        heights=c(2.5,1))

ggsave(Figure_S6,filename=paste0(path_save_figures,"FigureS6.png"), width=16, height = 18, units = "cm", dpi='print')


```


------------------------------------------------------------------------
### Figure S4 (homopolymers)

Homopolymers analysis
```{r, eval= F}
Figure_S4 <- ggplot(n_correct_consensus_average %>% filter(method %in% c("SINGLe","No weights","Guppy")),
                                    aes(x=Nreads,y=success_rate,col=method,lty=homopolymers))+
    geom_line()+facet_grid(method~variant)+
    ylab("Success rate")+geom_hline(aes(yintercept=0.9),col=grey(.5), lty="dashed")+
    xlab("Reads")+
    scale_color_manual(values=colors_vector_capital[c(1,3,5)],name="Method")+
    scale_linetype_manual(values=c("solid","dashed"),name="Homopolymers",labels=c("Original","Sorted"))+
    scale_y_continuous(breaks=c(0,.5,1), limits = c(0,1))+scale_x_continuous(breaks=c(10,20,30,40,50))+
    theme_plot+
    theme(legend.position = "bottom",legend.spacing.y = unit(0,"line"),
          legend.key.size = unit(1.5,"line"),
          plot.margin = unit(c(0,5.5,5.5,5.5),"points"),
          strip.background = element_rect(fill= "white"))
Figure_S4

ggsave(Figure_S4,filename=paste0(path_save_figures,"FigureS4.png"), width=17, height = 10, units = "cm", dpi='print')

```




------------------------------------------------------------------------
# KlenTaq large library 


### Filenames
```{r}
file_out_consensus<- paste0(path_analysis_lib,"/KTLib_ConsensusSINGLe.txt")
file_out_consensus_topten  <- paste0(path_analysis_lib,"KTLib_ConsensusbySINGLe_toptenBC.txt")
file_consensus_KTLib <- paste0(path_analysis_lib, "KTLib_consensus_table.txt")
```


## Basecall, align, detect barcodes, SINGLe
Basecalling
```{bash LIBbasecall, eval=F}
guppy_basecaller -i 2_NanoporeRawReads/LargeSet_fast5 -s 3_Analysis/LargeSet_SUP -c dna_r9.4.1_450bps_sup.cfg -x 'cuda:0' -r
#Filter sequences by length (1700 to 2100 bp)
mkdir 3_Analysis/KTLib/
awk 'NR % 4 ==2 && length($0) > 1700 && length($0) <2100 {print f  "\n" $0;getline;print;getline;print}  {f=$0}' 3_Analysis/LargeSet_SUP/pass/*fastq > 3_Analysis/KTLib/KTLib_1700_2100.fastq
```

Keep reads with mean Qscore\>15
```{r LIBfilter Q, eval=F}
seqsum_ktlib <- read.table("3_Analysis/LargeSet_SUP/sequencing_summary.txt", header=T) %>%
    filter(passes_filtering & mean_qscore_template>15) %>% 
    select(read_id)

file_ktlib_fastq <- "3_Analysis/KTLib/KTLib_1700_2100.fastq"#LargeSet_1700_2100.fastq"
file_ktlib_q15_fastq  <- "3_Analysis/KTLib/KTLib_1700_2100_Q15.fastq"
nlines <- as.numeric(str_split(system2(paste0("wc -l ",file_ktlib_fastq), intern=T), pattern=" ")[[1]][1])
connectionFile <- file(file_ktlib_fastq,"r")
n<-0
pb <- txtProgressBar(0,nlines,style = 3)
while(n < nlines){
  line <- readLines(connectionFile,n=1)
  id <- substr(line,2,37)
  if(id %in% seqsum_ktlib$read_id){
    cat(line,"\n",file = file_ktlib_q15_fastq, append=T, sep="")
    line <- readLines(connectionFile,n=1)
    cat(line,"\n",file = file_ktlib_q15_fastq, append=T, sep="")
    line <- readLines(connectionFile,n=1)
    cat(line,"\n",file = file_ktlib_q15_fastq, append=T, sep="")
    line <- readLines(connectionFile,n=1)
    cat(line,"\n",file = file_ktlib_q15_fastq, append=T, sep="")
  }else{
    line <- readLines(connectionFile,n=3)
  }
    n <- n+4
  setTxtProgressBar(pb,n)
}
 close(connectionFile)
```

Alignment using minimap2
```{bash LIBminimap2, eval=F}
minimap2 -ax map-ont 1_ReferenceFiles/KTB_digested.fasta 3_Analysis/KTLib/KTLib_1700_2100_Q15.fastq > 3_Analysis/KTLib/KTLib_1700_2100_Q15.sam

samtools view -S -b 3_Analysis/KTLib/KTLib_1700_2100_Q15.sam > 3_Analysis/KTLib/KTLib_1700_2100_Q15.bam #TRANSFORM TO BAM
samtools sort 3_Analysis/KTLib/KTLib_1700_2100_Q15.bam -o 3_Analysis/KTLib/KTLib_1700_2100_Q15.sorted.bam #SORT BAM FILE

```

SINGLe correction
```{r LIBSINGLe, eval=F}
tic <- proc.time()
pos_start <- 103
pos_end <- 1764

single_evaluate(bamfile    = "3_Analysis/KTLib/KTLib_1700_2100_Q15.sorted.bam",
                single_fits   = single_train,
                output_file = "3_Analysis/KTLib/KTLib_1700_2100_Q15_SINGLe.fastq",
                refseq_fasta = reference_file,
                pos_start=pos_start,pos_end=pos_end,
                verbose=T,gaps_weights="minimum",save=TRUE, 
                save_original_scores=TRUE)

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```


Barcodes detection. First create bash file.
```{r LIB detect barcodes, eval=F}
tic <-  proc.time()
## Inputs
bc_before   <- "actggc"
bc_after    <- "aagcc"
bc_before   <- toupper(bc_before)
bc_after    <- toupper(bc_after)
l           <- nchar(bc_before)
length_bc   <- 36
length_deltabc <- 2

inFile    <- "KTLib_1700_2100_Q15.fastq"
out.file  <- "KTLib_1700_2100_barcodes.txt"
code.file <- "3_Analysis/KTLib/FindBarcodes.sh"

bc_min <- length_bc-length_deltabc
bc_max <- length_bc+length_deltabc

## All possible variations of barcodes
bc_before_compl <- dna_reverse_complement_string(bc_before)
bc_after_compl <- dna_reverse_complement_string(bc_after)

bc_before_all       <- bc_one_mutation(bc_before)
bc_before_compl_all <- bc_one_mutation(bc_before_compl)
bc_after_all        <- bc_one_mutation(bc_after)
bc_after_compl_all  <- bc_one_mutation(bc_after_compl)

bc_combinations <- rbind(data.frame(before=bc_before,after= bc_after, forw=TRUE),
                         data.frame(before=bc_before_all, after=bc_after, forw=TRUE),
                         data.frame(before=bc_before, after=bc_after_all, forw=TRUE), 
                         data.frame(before=bc_after_compl, after=bc_before_compl, forw=FALSE),
                         data.frame(before=bc_after_compl, after=bc_before_compl_all, forw=FALSE),
                         data.frame(before=bc_after_compl_all, after=bc_before_compl, forw=FALSE))


## Main
if(file.exists(code.file)){file.remove(code.file)} ; file.create(code.file)
cat("#!/bin/bash\n\n",file=code.file)
cat("echo \"Name\tLineInFile\tBCsequence\tBCposition\tDirectionRead\" > ",out.file," \n", file=code.file, append=FALSE)

#2. Sequences with barcodes 
for(i in 1:nrow(bc_combinations)){
  bc_before <- as.character(bc_combinations$before[i])
  bc_after <- as.character(bc_combinations$after[i])
  forward  <- bc_combinations$forw[i]
  
  awk_in  <- paste0("awk '")
  awk_match <- paste0("NR%4==2 && match($0,/", bc_before, "([ACGT]{", bc_min,",",bc_max,"})", bc_after,"/)")
  if(forward){
    awk_print <- paste0(" {print a \"\t\" NR \"\t\" substr($0,RSTART+",l,",RLENGTH-",2*l,") \"\t\" RSTART+",l," \"\t+\" } {a=$0}' ",inFile, ">> ", out.file)
  }else{
    awk_print <- paste0(" {print a \"\t\" NR \"\t\" substr($0,RSTART+",l,",RLENGTH-",2*l,") \"\t\" RSTART+",l," \"\t-\" } {a=$0}' ",inFile, ">> ", out.file)
  }
  y_awk   <- paste0(awk_in, awk_match, awk_print)
  cat(y_awk,"\n", file=code.file, append=TRUE)
}

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)
```

Then run bash file.
```{bash, eval=F}
cd 3_Analysis/KTLib/
chmod +x FindBarcodes.sh
./FindBarcodes.sh
```

Filter barcodes that appear more than 3 times
```{r LIBfilter BC>3, eval=F}
tic <- proc.time()
barcodes_table  <- read.table(paste0(path_analysis_lib,"KTLib_1700_2100_barcodes.txt"),header=T, sep="\t")
barcodes_counts <- barcodes_table %>% group_by(BCsequence) %>% summarise(counts=n()) %>% filter(counts>3)
barcodesID      <- unique(barcodes_counts$BCsequence); n_bc <- length(barcodesID)
barcodes_table_fourcounts <- barcodes_table %>% 
    filter(BCsequence %in% barcodesID)%>%
    mutate(readID=str_sub(Name, 2, 37)) ;
nrow_bc <- nrow(barcodes_table_fourcounts)
write.table(barcodes_table_fourcounts,file=paste0(path_analysis_lib,"KTLib_1700_2100_barcodes_4reads.txt"))
toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

## Consensus for all barcode by all methods

Load data 
```{r, eval=FALSE}
barcodes_4reads <- read.table(file=paste0(path_analysis_lib,"KTLib_1700_2100_barcodes_4reads.txt")) %>%
    select(readID, BCsequence) %>% arrange(BCsequence)

barcodes_unique <- unique(barcodes_4reads$BCsequence)

```



Generate fasta and fastq files that we need
```{r create fast files, eval=FALSE}
tic <- proc.time()
dir.create( paste0(path_analysis_lib, "fastq_per_barcode"))
dir.create( paste0(path_analysis_lib, "fasta_per_barcode"))

KTLib_Q15_seqs <- readQualityScaledDNAStringSet(paste0(path_analysis_lib,"KTLib_1700_2100_Q15.fastq"))
names(KTLib_Q15_seqs) <- str_sub(names(KTLib_Q15_seqs) , 1, 36)
for(s in 1:nrow(barcodes_4reads)){
    seq_name <- barcodes_4reads$readID[s]
    bc_name <- barcodes_4reads$BCsequence[s]
    file_name_q <- paste0(path_analysis_lib, "fastq_per_barcode/", bc_name,".fastq")
    file_name_a <- paste0(path_analysis_lib, "fasta_per_barcode/", bc_name,".fasta")
    # which(names(KTLib_Q15_seqs)==seq_name)
    writeQualityScaledXStringSet(KTLib_Q15_seqs[seq_name], filepath = file_name_q, append = T)
    writeXStringSet(KTLib_Q15_seqs[seq_name], filepath = file_name_a, append = T)
}
toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```


Load corrected sequences
```{r LIBconsensus load seqs, eval=F}

tic <- proc.time()
seqs_single_ktlib <- readQualityScaledDNAStringSet(paste0(path_analysis_lib,"/KTLib_1700_2100_Q15_SINGLe.fastq"))
aux_seqs <- vapply(as.character(seqs_single_ktlib),strsplit, split="")
aux_quals <- 1-as(quality(seqs_single_ktlib),"NumericList")
aux_pos <- lapply(aux_seqs, seq)
aux_names <- c(vapply(names(seqs_single_ktlib),rep, each=1662))
data_single_ktlib <- data.frame(nucleotide = unlist(aux_seqs),
                        p_SINGLe=unlist(aux_quals),
                        position=unlist(aux_pos),
                        readID=aux_names) 
rownames(data_single_ktlib)<-NULL

seqs_guppy_ktlib <- readQualityScaledDNAStringSet(paste0(path_analysis_lib,"/KTLib_1700_2100_Q15_SINGLe_original.fastq"))
aux_seqs <- vapply(as.character(seqs_single_ktlib),strsplit, split="")
aux_quals <- 1-as(quality(seqs_single_ktlib),"NumericList")
aux_pos <- lapply(aux_seqs, seq)
aux_names <- c(vapply(names(seqs_single_ktlib),rep, each=1662))
data_guppy_ktlib <- data.frame(nucleotide = unlist(aux_seqs),
                        p_Guppy=unlist(aux_quals),
                        position=unlist(aux_pos),
                        readID=aux_names)
rownames(data_guppy_ktlib)<-NULL


data_single_ktlib <- data_single_ktlib %>% 
    left_join(data_guppy_ktlib, by = c("nucleotide", "position", "readID"))%>%
    left_join(barcodes_4reads, by = "readID")%>%
    filter(!is.na(BCsequence))  %>%
    mutate(isWT = nucleotide==wildtype[position]) %>%
    mutate(p_SINGLe=ifelse(is.na(p_SINGLe), 1, p_SINGLe))%>%
    mutate(p_noweights=1)


toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

barcodes_unique <- unique(data_single_ktlib$BCsequence)
n_bc <- length(barcodes_unique)

# save.image(file='myEnvironment.RData')
```

```{r, eval=F}
# load('myEnvironment.RData')
```

By SINGLe (sorting homopolymers)
```{r LIBconsensus by single, eval=F}
tic <- proc.time()
#Compute consensus for each barcode using all sequences available and using the three available methods

cat("Barcode;Nreads;mutations_by_single;mutations_by_nanopore;mutations_no_weights\n", file=file_out_consensus, append=F)

pb <- txtProgressBar(0,n_bc,style=3)
counter<-0
wt_homopolymers <- detect_homopolymer_positions(wildtype)
wt_homopolymers <- wt_homopolymers[vapply(wt_homopolymers,length)>1]
pos_wt_homopolymers <- unlist(wt_homopolymers)

for( bc in barcodes_unique){
  bc_reads <- data_single_ktlib %>%
    filter(BCsequence==bc) 
  
  index_homopolymers <- which(bc_reads$position %in% pos_wt_homopolymers & bc_reads$nucleotide=="-")
  n_memory <- c()
  for(n in index_homopolymers){
    if(n %in% n_memory){next()}
    hp_pos <- bc_reads$position[n]
    hp_pos_ref <- wt_homopolymers[vapply(wt_homopolymers, function(x){hp_pos %in% x})][[1]]
    hp_vec_logical <- hp_pos_ref==hp_pos
    hp_vec_length <- length(hp_vec_logical)
    hp_vec <- rep(NA,hp_vec_length)
    ind_aux <- which(hp_vec_logical)

    hp_vec[1:ind_aux] <- rev(n+1-1:ind_aux)
    if(ind_aux<hp_vec_length){
      hp_vec[(ind_aux+1):hp_vec_length] <- n+c(1:(hp_vec_length-ind_aux))
    }
    n_memory <- hp_vec
    
    data_aux <- bc_reads[hp_vec,]
    if(all(data_aux$nucleotide=="-")){next()}
    
    pos <- data_aux$position
    nucl <- data_aux$nucleotide
    ind_order <- c(which(nucl!="-"),which(nucl=="-"))
    
    bc_reads$position[hp_vec] <- pos[order(ind_order)]
    
    if(length( which(bc_reads$position %in% pos_wt_homopolymers & bc_reads$nucleotide=="-"))<length(index_homopolymers)){
      cat("n=",n,"\n")
      stop()
    }
  }
  bc_reads <- bc_reads %>% arrange(readID,position)

  cons_single <- weighted_consensus(df = bc_reads %>%
                          select(nucleotide,p_SINGLe,position), 
                           cutoff_prob = 0)
  cons_nanopore <- weighted_consensus(df = bc_reads %>%
                          select(nucleotide,p_Guppy,position), 
                           cutoff_prob = 0)
  cons_noweight <- weighted_consensus(df = bc_reads %>%
                          select(nucleotide,p_noweights,position), 
                           cutoff_prob = 0)

  muts_single <- c(detect_mutations(biostr_to_char(cons_single),  wildtype))
  muts_nanopore <- c(detect_mutations( biostr_to_char(cons_nanopore), wildtype))
  muts_noweight <- c(detect_mutations( biostr_to_char(cons_noweight),  wildtype))

  n_single <- length(muts_single)
  n_nanopore <- length(muts_nanopore)
  n_noweight <- length(muts_noweight)
  
  if(n_single==0){n_single=1;muts_single="wildtype"}
  if(n_nanopore==0){n_nanopore=1;muts_nanopore="wildtype"}
  if(n_noweight==0){n_noweight=1;muts_noweight="wildtype"}
  n_tot = n_single+n_nanopore+n_noweight
      
  df <- data.frame(barcode = rep(bc,n_tot),nseqs=length(unique(bc_reads$readID)),
                  mutation=c(muts_single,muts_nanopore,muts_noweight),
                  method = c(rep("single", n_single),
                             rep("nanopore",n_nanopore),
                             rep("no_weights",n_noweight)))
  
  write.table(df, append = T, col.names = F,row.names = F, file = file_out_consensus)
  counter <- counter+1
  setTxtProgressBar(pb,counter)
}
toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

# save.image(file='myEnvironment.RData')

```


Prepare fast5 files for Nanopolish. Please indicate proper paths indicated in capital letters
```{r LIBprepare fast5, eval=F}
tic <- proc.time()
dir.create(paste0(path_analysis_lib,"/fast5files/") )# individual fast5 files will be stored there temporally 

split_fast5 <- paste0("multi_to_single_fast5 -i 2_NanoporeRawReads/LargeSet_fast5 -s 3_Analysis/KTLib/fast5files/")  #PATH_TO_FAST5_FILES location of nanopore reads
system2(split_fast5)

# Place them in folders according to barcode
single_files       <- dir(paste0(path_analysis_lib,"/fast5files/"), full.names = T, recursive = T)
single_files_name  <- dir(paste0(path_analysis_lib,"/fast5files/"), full.names = F, recursive = T)

single_files <- data.frame(file=single_files, readID=single_files_name) %>%
  mutate(readID = sub(readID, pattern=".fast5",replacement="")) %>%
  mutate(readID = sub(readID, pattern="[[:digit:]]+/", replacement=""))

barcodes_table_fourcounts <- barcodes_table %>% left_join(single_files, by="readID")
nrow_bc <- nrow(barcodes_table_fourcounts)
pb <- txtProgressBar(min = 0, max = nrow_bc, style = 3)
for(i in 1:nrow_bc){
  setTxtProgressBar(pb,i)
  seqid <- barcodes_table_fourcounts$readID[i]
  barcode <- barcodes_table_fourcounts$BCsequence[i]
  barcode_dir <- paste0(path_analysis_lib,"fast5files/",barcode)
  file <- barcodes_table_fourcounts$file[i]
  if(!dir.exists(barcode_dir)){ dir.create(barcode_dir)}
  command <- paste0("cp ", file, " ",barcode_dir)
  system2(command)
}

# Create unique fast5 file per barcode
pb <- txtProgressBar(min = 0, max = nrow_bc, style = 3)
for(bc in barcodesID){
  setTxtProgressBar(pb,bc)
  command <- paste0("single_to_multi_fast5 -i 3_Analysis/KTLib/fast5files/",bc," -s 3_Analysis/KTLib/fast5files_byBC/ -f ",bc)
  system2(command)
}

barcodes_table_fourcounts <- barcodes_table_fourcounts %>% select(-file)
rm(single_files,path_fast5,path_fast5_individual,path_fast5_individual_out)

### Create fastq file for interesting barcodes
dir.create(paste0(path_analysis_lib,"/fastq/"))
fastq_file <- readLines(paste0(path_analysis_lib,"/KTLib_1700_2100_Q15.fastq"))
n_fastq <- length(fastq_file)
for(i in seq(1,n_fastq, by=4)){
  ind <- grep(pattern = substr(fastq_file[i],2,37),x=barcodes_table_fourcounts$SeqID)
  if(length(ind)==0){next(i,"not found \n")}
  write(fastq_file[c(i,i+1,i+2,i+3)],
        file=paste0(path_analysis_lib,"/fastq/",barcodes_table_fourcounts$BCsequence[ind],".fastq"),
        append=T)
}


### Create fasta from fastq
fastq_files <- dir(paste0(path_analysis_lib,"/fastq/"), full.names = T, pattern="*.fastq")
fasta_files <- str_replace_all(fastq_files, "fastq", "fasta")
dir.create(paste0(path_analysis_lib,"/fasta/"))
commands <- paste0("seqtk seq -a ", fastq_files, " > ",fasta_files)
vapply(commands, system2)

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

Compute using nanopolish and medaka. first generate bash file
```{r LIBconsensus others, eval=FALSE}

# Load data
tic <- proc.time()
#outputs
bash_file <- "ConsensusByMethods.sh"
results_medaka     <- paste0(path_analysis_lib,"ConsensusbyMedaka.txt")
nanopolish_results <- paste0(path_analysis_lib,"Consensus_by_Nanopolish.txt")

for(bc in barcodes_unique){
  #all filenames
  filename <- bc
  fasta      <- paste0(path_analysis_lib,"fasta/",bc,".fasta")
  fastq      <- paste0(path_analysis_lib,"fastq/",bc,".fastq")
  sam        <- paste0(path_analysis_lib,bc,".sam")
  sorted     <- paste0(path_analysis_lib,bc,".sorted.fasta")
  racon      <- paste0(path_analysis_lib,bc,".racon.fasta")
  medaka     <- paste0(path_analysis_lib,bc,"_medaka.fasta")
  nanopolish <- paste0(path_analysis_lib,bc,"_nanopolish.txt")
      
  ## Run minimap2
  line_minimap <- paste0("minimap2 -ax map-ont ", reference_file, " ", fastq," > ", sam)
  
  ## Run racon
  line_racon <- paste0("racon ",fastq," ",sam, " ",reference_file, " >> ", racon)
  line_remove_racon <- paste0("rm ", racon,".fai ",racon, ".mmi ",racon  )
  
  ## Run Medaka
  line_medaka <- paste0("medaka_consensus -i ",fastq," -d ",racon,"  -o ", medaka," -t 4 ")
  line_results_medaka <- paste0("echo '>",filename, "' >> ",results_medaka)
  line_results_medaka2 <- paste0("tail -n 1 ",filename,"_medaka.fasta/consensus.fasta", " >> ",results_medaka)
  line_remove_medaka <- paste0("rm -r ",filename,"_medaka.fasta/")
  
  ## Run nanopolish
  line_nanopolish_index  <- paste0("nanopolish index -d ",bc,"/ ",fastq)
  line_nanopolish_sort    <- paste0("samtools sort ",sam," -o ",sorted," -T temporal.tmp ")
  line_nanopolish_index2 <- paste0("samtools index ", sorted)
  line_nanopolish_variants <- paste0("nanopolish variants --consensus -o ", nanopolish, " -r ",fastq," -b ",sorted," -g ", reference_file)
  line_results_nanopolish <- paste0("awk '!/^#/ {print \"",filename,"\\t\" $0}' ",nanopolish, " >> ",nanopolish_results)
  line_remove_nanopolish <- paste0("rm ", fastq,".index* ", filename,".sorted* ", nanopolish, " ", sam)

  cat(line_minimap,
      line_racon,
      line_medaka,
      line_results_medaka,
      line_results_medaka2,
      line_nanopolish_index,
      line_nanopolish_sort,
      line_nanopolish_index2,
      line_nanopolish_variants,
      line_results_nanopolish,
      line_remove_racon,
      line_remove_medaka,
      line_remove_nanopolish,
      file=bash_file, sep="\n", append=T)
}
toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

Then run it
```{bash LIBconsensus others bash, eval=F}
chmod +x ConsensusByMethods.sh
./ConsensusByMethods.sh
```


Align consensus sequences to reference
```{bash, eval=F}
minimap2 -ava-ont 1_ReferenceFiles/KTrefORF_1662.fasta ConsensusbyMedaka.fasta > ConsensusbyMedaka_aligned.sam
minimap2 -ava-ont 1_ReferenceFiles/KTrefORF_1662.fasta ConsensusbyRacon.fasta > ConsensusbyRacon_aligned.sam
```


Parse results into a table consensus_mismatches

```{r LIBconsensus parse, eval=F}
tic <- proc.time()
pos_initial <- 1
pos_final <- length(wildtype)

cigar_reference_counts   <- c("M","D","N","=","X")
cigar_query_counts <- c("M","I","S","=","X")

sam_data_medaka <- read.table(paste0(path_analysis_lib,"ConsensusbyMedaka_aligned.sam"), skip=2)
medaka_mismatches     <- data.frame(Barcode=NA,Position=NA,base.original=NA,base.mutated=NA)
for(i in 1:nrow(sam_data_medaka)){

  line_data       <- unlist(sam_data_medaka[i,])
  pos_ini_aliref  <-  as.numeric(line_data[4])
  cigar           <-  line_data[6]
  sequence        <-  strsplit(line_data[10],split="")[[1]]

  #Obtain CIGAR vector
  cigar.table      <- data.frame(n=as.numeric(strsplit(cigar, split="[MIDNSHP=X]")[[1]]),Op= strsplit(cigar, split="[0-9]+")[[1]][-1])
  cigar_vec        <- unlist(apply(cigar.table,1, function(x){rep(x[2],x[1])}))
  names(cigar_vec) <- NULL
  if(pos_ini_aliref>1){cigar_vec <- c(rep("D",pos_ini_aliref-1),cigar_vec)}
  
  #Data frame for the original sequence, and reference for the model (full)
  mismatches            <- data.frame(Barcode =rep(line_data[1], length(cigar_vec)),
                                      cigar=cigar_vec, 
                                      base.mutated=NA,pos_in_samalig=NA,Position=NA,base.original=NA)
  
  rows_with_query     <- mismatches$cigar%in%cigar_query_counts
  rows_with_reference <- mismatches$cigar%in%cigar_reference_counts
  
  mismatches$base.mutated[rows_with_query]  <- sequence     #Nucleotide
  mismatches$base.mutated[is.na(mismatches$base.mutated)] <- "-"

  mismatches$pos_in_samalig[rows_with_reference] <- pos_ini_aliref + c(0:(sum(rows_with_reference)-1))  #position - according to aligning reference
  
  mismatches$Position[rows_with_reference] <- 1:length(wildtype) #position - according to single reference
  mismatches$base.original[rows_with_reference]   <- wildtype

  index <- which(mismatches$base.mutated!=mismatches$base.original)
  
  mismatches <- mismatches[index,c("Barcode", "Position","base.original","base.mutated")]
  medaka_mismatches <- medaka_mismatches %>% rbind(mismatches) 

}
medaka_mismatches <- medaka_mismatches %>% filter(!is.na(Barcode))  %>% mutate(method="Medaka",mutation=paste0(base.original,Position,base.mutated))

sam_data_racon <- read.table("ConsensusbyRacon_aligned.sam", skip=2)
racon_mismatches     <- data.frame(Barcode=NA,Position=NA,base.original=NA,base.mutated=NA)
for(i in 1:nrow(sam_data_racon)){

  line_data       <-  unlist(sam_data_racon[i,])
  pos_ini_aliref  <-  as.numeric(line_data[4])
  cigar           <-  line_data[6]
  sequence        <-  strsplit(line_data[10],split="")[[1]]

  #Obtain CIGAR vector
  cigar.table      <- data.frame(n=as.numeric(strsplit(cigar, split="[MIDNSHP=X]")[[1]]),Op= strsplit(cigar, split="[0-9]+")[[1]][-1])
  cigar_vec        <- unlist(apply(cigar.table,1, function(x){rep(x[2],x[1])}))
  names(cigar_vec) <- NULL
  if(pos_ini_aliref>1){cigar_vec <- c(rep("D", pos_ini_aliref-1), cigar_vec)}
  
  #Data frame for the original sequence, and reference for the model (full)
  mismatches            <- data.frame(Barcode =rep(line_data[1], length(cigar_vec)),
                                      cigar=cigar_vec, 
                                      base.mutated=NA,pos_in_samalig=NA,Position=NA,base.original=NA)
  
  rows_with_query     <- mismatches$cigar%in%cigar_query_counts
  rows_with_reference <- mismatches$cigar%in%cigar_reference_counts
  
  mismatches$base.mutated[rows_with_query]  <- sequence     #Nucleotide
  mismatches$base.mutated[is.na(mismatches$base.mutated)] <- "-"
  
  mismatches$pos_in_samalig[rows_with_reference] <- pos_ini_aliref + c(0:(sum(rows_with_reference)-1))  #position - according to aligning reference
  
  mismatches$Position[rows_with_reference] <- 1:length(wildtype) #position - according to single reference
  mismatches$base.original[rows_with_reference]   <- wildtype

  index <- which(mismatches$base.mutated!=mismatches$base.original)
  
  mismatches <- mismatches[index,c("Barcode", "Position","base.original","base.mutated")]
  racon_mismatches <- racon_mismatches %>% rbind(mismatches) 

}
racon_mismatches <- racon_mismatches %>% filter(!is.na(Barcode)) %>% mutate(method="Racon", mutation=paste0(base.original,Position,base.mutated))

results_nanopolish <- readLines("ConsensusbyNanopolish.txt")
results_nanopolish <- do.call(rbind,vapply(results_nanopolish, strsplit, split="\t"))
colnames(results_nanopolish) <- rownames(results_nanopolish) <- NULL
results_nanopolish <- results_nanopolish[,c(1,3,5,6)]
colnames(results_nanopolish) <- c("Barcode","Position","base.original","base.mutated")
results_nanopolish <- as.data.frame(results_nanopolish)
results_nanopolish$Position <- as.numeric(results_nanopolish$Position)

results_nanopolish <- results_nanopolish %>% filter(!is.na(Barcode)) 

#remove insertions
ind_insertions <- nchar(results_nanopolish$base.mutated)>1
results_nanopolish$base.mutated[ind_insertions] <- vapply(results_nanopolish$base.mutated[ind_insertions] , substr,1,1)
results_nanopolish <- results_nanopolish %>% filter(base.mutated!=base.original)

#Split deletions
index_aux <- which(nchar(results_nanopolish$base.original)==2)
results_nanopolish$base.original[index_aux] <- substr(results_nanopolish$base.original[index_aux],2,2)
results_nanopolish$base.mutated[index_aux] <- "-"
results_nanopolish$Position[index_aux] <- as.numeric(results_nanopolish$Position[index_aux])+1

index_aux_3 <- which(nchar(results_nanopolish$base.original)==3)
aux_df_1 <- results_nanopolish[index_aux_3,]
aux_df_1$Position <- as.numeric(aux_df_1$Position)+1
aux_df_1$base.original <- substr(aux_df_1$base.original, 2,2)
aux_df_1$base.mutated <- "-"

aux_df_2 <- results_nanopolish[index_aux_3,]
aux_df_2$Position <- as.numeric(aux_df_2$Position)+2
aux_df_2$base.original <- substr(aux_df_2$base.original, 3,3)
aux_df_2$base.mutated <- "-"

results_nanopolish <- results_nanopolish[-index_aux_3,]
results_nanopolish <- results_nanopolish %>% rbind(aux_df_1) %>% rbind(aux_df_2)
rm(aux_df_1,aux_df_2)

index_aux_4 <- which(nchar(results_nanopolish$base.original)==4)
aux_df_1 <- results_nanopolish[index_aux_4,]
aux_df_1$Position <- as.numeric(aux_df_1$Position)+1
aux_df_1$base.original <- substr(aux_df_1$base.original, 2,2)
aux_df_1$base.mutated <- "-"

aux_df_2 <- results_nanopolish[index_aux_4,]
aux_df_2$Position <- as.numeric(aux_df_2$Position)+2
aux_df_2$base.original <- substr(aux_df_2$base.original, 3,3)
aux_df_2$base.mutated <- "-"

aux_df_3 <- results_nanopolish[index_aux_4,]
aux_df_3$Position <- as.numeric(aux_df_3$Position)+3
aux_df_3$base.original <- substr(aux_df_3$base.original, 4,4)
aux_df_3$base.mutated <- "-"

results_nanopolish <- results_nanopolish[-index_aux_3,]
results_nanopolish <- results_nanopolish %>% rbind(aux_df_1) %>% rbind(aux_df_2)%>% rbind(aux_df_3) %>% arrange(Barcode,Position) %>% mutate(method="Nanopolish", mutation=paste0(base.original,Position,base.mutated))


save(results_nanopolish, medaka_mismatches, racon_mismatches, file=paste0(path_consensus,"Mismatches_by_Method.Rdata"))


single_mismatches             <- read.table("ConsensusBySINGLE.txt", header=F,skip=1)
colnames(single_mismatches) <- colnames(single_mismatches_sortHPinVCS) <- c("Barcode","Nreads","mutation","method")

single_mismatches <- single_mismatches %>% 
  mutate(base.original = stringr::str_sub(mutation,1,1))     %>%
  mutate(base.mutated = stringr::str_sub(mutation,-1,-1))    %>%
  mutate(Position = stringr::str_sub(mutation,2,-2))         %>%
  select(Barcode,base.original,base.mutated,Position,method,mutation) 


KTLib_mutations <- rbind(single_mismatches,single_mismatches_sortHPinVCS, results_nanopolish, medaka_mismatches,racon_mismatches)
write.table(KTLib_mutations, file=file_consensus_KTLib)

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```


## Analyse convergence of 10 most frequent barcodes

Identify ten most frequent barcodes
```{r LIB10BCs, eval=F}
tic <- proc.time()
KTLib_topBC <- barcodes_4reads %>% 
    group_by(BCsequence) %>% 
    summarise(n=n()) %>%
    arrange(-n)%>%
    head(n=10)
KTLib_topBC_reads<- barcodes_4reads %>% 
    filter(BCsequence %in% KTLib_topBC$BCsequence) 
toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```


Create a fastq file for each of them
```{r LIB10BCs create fastq, eval=F}
tic <- proc.time()
for(i in unique(KTLib_topBC$BCsequence)){
    reads_bc_i <- KTLib_topBC_reads %>% filter(BCsequence==i)
    sequences_bc_i <- seqs_single_ktlib [names(seqs_single_ktlib) %in% reads_bc_i$readID]
    writeQualityScaledXStringSet(sequences_bc_i, file=paste0(path_analysis_lib, i, ".fastq"))  
}
toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

Compute consensus on subsets of reads by SINGLE
```{r LIB10BCs consensus single, eval=F}
tic <- proc.time()
#Load corrected data
data_single_ktlib_topten <- data_single_ktlib %>%
    filter(Name %in% KTLib_topBC_reads$Name ) 

# # Load data
# reads_corrected <- read.table("3_PreprocessedData/SINGLed_data/LargeSet_F_corrected.txt", header=T)
# reads_corrected <- reads_corrected %>% filter(SeqID %in% KTLib_BC_reads$SeqID )
# reads_corrected_rev <- read.table("3_PreprocessedData/SINGLed_data/LargeSet_R_corrected.txt", header=T)
# reads_corrected_rev <- reads_corrected_rev %>% filter(SeqID %in% KTLib_BC_reads$SeqID )
# reads_corrected_all     <- rbind(reads_corrected,reads_corrected_rev)
#                          
# reads_corrected <- left_join(reads_corrected_all, KTLib_topBC_reads %>% select(SeqID,BCsequence) , by="SeqID")
# reads_corrected <- reads_corrected %>% filter(!is.na(reads_corrected$BCsequence))
# reads_corrected$p_right_priors_model[is.na(reads_corrected$p_right_priors_model)]<-1


#Compute consensus for each barcode using all sequences available and using the three available methods
cat("Barcode\tNreads\trepetition\tmutation\tmethod\n", file=file_out_consensus_topten, append=F)

barcodes_unique <- unique(data_single_ktlib_topten$BCsequence)
n_bc <- length(barcodes_unique)

for(bc in barcodes_unique){
    cat(bc, "\n")
  bc_reads <- data_single_ktlib_topten %>%
    filter(BCsequence==bc)
  
  #compute consensus on subsets
  bc_seqsID <- unique(bc_reads$readID)
  for( s in subsets){
    cat(s, "\t")
    counter<- 0
    pb <- txtProgressBar(0,n_repetitions,style=3)
    for (r in 1:n_repetitions) {
      counter <- counter+1
      setTxtProgressBar(pb,counter)
  
      subset_bc_reads <- bc_reads %>% 
          filter(readID %in% sample(bc_seqsID, s))
      
    #sort homopolymers
    index_homopolymers <- which(subset_bc_reads$position %in% pos_wt_homopolymers
                                & subset_bc_reads$nucleotide=="-")
    n_memory <- c()
    for(n in index_homopolymers){
        if(n %in% n_memory){next()}
        hp_pos <- subset_bc_reads$position[n]
        hp_pos_ref <- wt_homopolymers[vapply(wt_homopolymers, function(x){hp_pos %in% x})][[1]]
        hp_vec_logical <- hp_pos_ref==hp_pos
        hp_vec_length <- length(hp_vec_logical)
        hp_vec <- rep(NA,hp_vec_length)
        ind_aux <- which(hp_vec_logical)
        
        hp_vec[1:ind_aux] <- rev(n+1-1:ind_aux)
        if(ind_aux<hp_vec_length){
          hp_vec[(ind_aux+1):hp_vec_length] <- n+c(1:(hp_vec_length-ind_aux))
        }
        n_memory <- hp_vec
        
        data_aux <- subset_bc_reads[hp_vec,]
        if(all(data_aux$nucleotide=="-")){next()}
        
        pos <- data_aux$position
        nucl <- data_aux$nucleotide
        ind_order <- c(which(nucl!="-"),which(nucl=="-"))
        
        subset_bc_reads$position[hp_vec] <- pos[order(ind_order)]
        
        if(length( which(subset_bc_reads$position %in% pos_wt_homopolymers & subset_bc_reads$nucleotide=="-"))<length(index_homopolymers)){
          cat("n=",n,"\n")
          stop()
        }
    }
    subset_bc_reads <- subset_bc_reads %>% arrange(readID,position)
    
    cons_single <- weighted_consensus(df = subset_bc_reads %>%
                                          select(nucleotide,p_SINGLe,position), 
                                   cutoff_prob = 0)
    cons_guppy  <- weighted_consensus(df = subset_bc_reads %>% 
                                          select(nucleotide,p_Guppy,position),
                                   cutoff_prob = 0)
    cons_noweight <- weighted_consensus(df = subset_bc_reads %>%
                                            select(nucleotide,p_noweights,position),
                                    cutoff_prob = 0)
    muts_single   <- detect_mutations(biostr_to_char(cons_single), wildtype)
    muts_guppy    <- detect_mutations(biostr_to_char(cons_guppy), wildtype)
    muts_noweight <- detect_mutations(biostr_to_char(cons_noweight), wildtype)

      n_single <- length(muts_single)
      n_guppy <- length(muts_guppy)
      n_noweight <- length(muts_noweight)
      n_tot <- n_single+n_guppy+n_noweight
      if(n_tot==0){next()}
      df <- data.frame(barcode = rep(bc,n_tot),nseqs=s,repetition=r,
                      mutation=c(muts_single,muts_guppy,muts_noweight),
                      method = c(rep("SINGLe", n_single),
                                 rep("Guppy",n_guppy), 
                                 rep("No weights",n_noweight)))
      write.table(df, append = T, col.names = F,row.names = F, file = file_out_consensus_topten)
    }
  }

}

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)
```

Compute consensus on subsets of reads by Nanopolish and Medaka
```{r LIB10BCs consensus others, eval=F}
tic <- proc.time()
dir.create(paste0(path_analysis_lib,"/subsets/"))
file.copy(reference_file, paste0(path_analysis_lib,"/subsets/"))
## Make fastq files with subsets
for(bc in unique(KTLib_topBC_reads$BCsequence)){
  
  fastq_file <- readLines(paste0(path_analysis_lib, bc,".fastq"))
  name_lines <- grep(pattern="^@",x=fastq_file)

  for(ns in subsets){
    if(ns> length(name_lines)){next()}
    for(r in 1:n_repetitions){
      fastq_subset_out    <- paste0(path_analysis_lib,"subsets/",bc,"_",ns,"_",r,".fastq")

      ## Make file of subset of sequences
      choose_lines <- sample(name_lines, size= ns)
      choose_lines_all <- sort(c(choose_lines,choose_lines+1,choose_lines+2,choose_lines+3))
      writeLines(text = fastq_file[choose_lines_all], con = file(fastq_subset_out))
    }
  }
}

in_files <- 
    dir(pattern=".fastq", path=paste0(path_analysis_lib,"/subsets/"))

## .sh for minimap2
for(i in seq_along(in_files)){
  fastq_file <- in_files[i]
  sam_file   <- sub(pattern=".fastq", replacement = ".sam", x=fastq_file)
  line_minimap <- paste0("minimap2 -ax map-ont KTrefORF_1662.fasta ",fastq_file," > ", sam_file)
  cat(line_minimap,file=paste0(path_analysis_lib,"/subsets/Minimap_allBC.sh"), sep="\n", append=i!=1)
}

# .sh for racon
results_racon      <- "../ConsensusbyRacon_subsets.txt"
for(i in seq_along(in_files)){
  fastq_file <- in_files[i]
  sam_file   <- sub(pattern=".fastq", replacement = ".sam", x=fastq_file)
  racon      <- sub(pattern=".fastq", replacement = ".racon.fasta", x=fastq_file)
  filename     <- sub(pattern=".fastq", replacement = "", x=fastq_file)
  
  ## Run racon
  line_racon <- paste0("racon ",fastq_file," ",sam_file, " KTrefORF_1662.fasta >> ", racon)
  line_results_racon <- paste0("echo '>",filename, "' >> ",results_racon)
  line_results_racon2 <- paste0("tail -n 1 ",racon, " >> ",results_racon)

  cat(line_racon,line_results_racon,line_results_racon2,file=paste0(path_analysis_lib,"/subsets/Racon_allBC.sh"), sep="\n", append= i!=1) 
}

# .sh for Medaka
for(i in seq_along(in_files)){
  fastq_file <- in_files[i]
  racon      <- sub(pattern=".fastq", replacement = ".racon.fasta", x=fastq_file)
  filename   <- sub(pattern=".fastq", replacement = "", x=fastq_file)
  medaka     <- sub(pattern=".fastq", replacement = "_medaka.fasta", x=fastq_file)
  results_medaka <- paste0("Consensus_",filename,".fasta")

  ## Run Medaka
  line_medaka <- paste0("medaka_consensus -i ",fastq_file," -d ",racon,"  -o ", medaka," -t 4 ")
  line_results_medaka <- paste0("echo '>",filename, "' >> ",results_medaka)
  line_results_medaka2 <- paste0("tail -n 1 ",filename,"_medaka.fasta/consensus.fasta", " >> ",results_medaka)
  line_remove_medaka <- paste0("rm -r ",medaka)
  
  cat(line_medaka,line_results_medaka,line_results_medaka2,line_remove_medaka,file=paste0(path_analysis_lib,"/subsets/Medaka_",filename,".sh"), sep="\n", append=F) 
}

# .sh for nanopolish
for(i in seq_along(in_files)){
  fastq_file         <- in_files[i]
  sam_file           <- sub(pattern=".fastq", replacement = ".sam", x=fastq_file)
  nanopolish_file    <- sub(pattern=".fastq", replacement = "_nanopolish.txt", x=fastq_file)
  filename           <- sub(pattern=".fastq", replacement = "", x=fastq_file)
  nanopolish_results <- paste0(filename,"_ConsensusbyNanopolish_subsets.txt")
  bc <- strsplit(filename, split="_")[[1]][1]
  
  fast5_dir <- paste0("../fast5/individualfiles/",bc)

  ## lines nanopolish
  line_nanopolish_index  <- paste0("nanopolish index -d  ../fast5/ subsets/",fastq_file)
  line_nanpolish_sort    <- paste0("samtools sort subsets/",sam_file," -o ",filename,".sorted.fasta -T temporal.tmp ")
  line_nanopolish_index2 <- paste0("samtools index ", filename,".sorted.fasta")
  line_nanopolish_variants <- paste0("nanopolish variants --consensus -o ", nanopolish_file, " -r BCtopten/",fastq_file," -b ",filename,".sorted.fasta -g ", reference_file)
  line_results_nanopolish <- paste0("awk '!/^#/ {print \"",filename,"\\t\" $0}' ",nanopolish_file, " >> ",nanopolish_results)
  line_remove_nanopolish <- paste0("rm ", filename,"_nanopolish.txt ", filename,".sorted* ")

  cat("cd ../ \n", 
      line_nanopolish_index,
      line_nanpolish_sort,
      line_nanopolish_index2,
      line_nanopolish_variants,
      line_results_nanopolish,
      line_remove_nanopolish,
      file=paste0("subsets/Nanopolish_",filename,".sh"), sep="\n", append=F)
}
toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

```{bash, eval=F}
cd 3_Analysis/KTLib/subsets
chmod +x Minimap_allBC.sh
./Minimap_allBC.sh

chmod +x Racon_allBC.sh
./Racon_allBC.sh

chmod +x Medaka*
.Medaka*  

chmod +x Nanopolish*
./Nanopolish*  

cat *_ConsensusbyRacon.txt > ConsensusbyRacon_subsets.txt
cat *_ConsensusbyMedaka.txt > ConsensusbyMedaka_subsets.txt
cat *_ConsensusbyNanopolish.txt > ConsensusbyNanopolish_subsets.txt

minimap2 -ava-ont ../../1_ReferenceFiles/KTrefORF_1662.fasta ConsensusbyRacon_subsets.fasta > ConsensusbyRacon_subsets_aligned.sam
minimap2 -ava-ont ../../1_ReferenceFiles/KTrefORF_1662.fasta ConsensusbyMedaka_subsets.fasta > ConsensusbyMedaka_subsets_aligned.sam
```

Identify mutations and save table
```{r}
file_consensus_top10_inLIB <- paste0(path_analysis_lib,"Convergence_10top_mismatches.txt" )

```

```{r LIB10BCs consensus to muts, eval=FALSE}
tic <- proc.time()
pos_initial <- 1
pos_final <- length(wildtype)

cigar_reference_counts   <- c("M","D","N","=","X")
cigar_query_counts <- c("M","I","S","=","X")

sam_data_medaka <- read.table("subsets/ConsensusbyMedaka_subsets_aligned.sam", skip=2)
medaka_mismatches_subsets     <- data.frame(Barcode=NA,Position=NA,base.original=NA,base.mutated=NA)
for(i in 1:nrow(sam_data_medaka)){
  line_data       <- unlist(sam_data_medaka[i,])
  pos_ini_aliref  <-  as.numeric(line_data[4])
  cigar           <-  line_data[6]
  sequence        <-  strsplit(line_data[10],split="")[[1]]

  #Obtain CIGAR vector
  cigar.table      <- data.frame(n=as.numeric(strsplit(cigar, split="[MIDNSHP=X]")[[1]]),Op= strsplit(cigar, split="[0-9]+")[[1]][-1])
  cigar_vec        <- unlist(apply(cigar.table,1, function(x){rep(x[2],x[1])}))
  names(cigar_vec) <- NULL
  if(pos_ini_aliref>1){cigar_vec <- c(rep("D",pos_ini_aliref-1),cigar_vec)}
  
  #Data frame for the original sequence, and reference for the model (full)
  mismatches            <- data.frame(Barcode =rep(line_data[1], length(cigar_vec)),
                                      cigar=cigar_vec, 
                                      base.mutated=NA,pos_in_samalig=NA,Position=NA,base.original=NA)
  
  rows_with_query     <- mismatches$cigar%in%cigar_query_counts
  rows_with_reference <- mismatches$cigar%in%cigar_reference_counts
  
  mismatches$base.mutated[rows_with_query]  <- sequence     #Nucleotide
  mismatches$base.mutated[is.na(mismatches$base.mutated)] <- "-"

  mismatches$pos_in_samalig[rows_with_reference] <- pos_ini_aliref + c(0:(sum(rows_with_reference)-1))  #position - according to aligning reference
  
  mismatches$Position[rows_with_reference] <- 1:length(wildtype) #position - according to single reference
  mismatches$base.original[rows_with_reference]   <- wildtype

  index <- which(mismatches$base.mutated!=mismatches$base.original)
  
  mismatches <- mismatches[index,c("Barcode", "Position","base.original","base.mutated")]
  medaka_mismatches_subsets <- medaka_mismatches_subsets %>% rbind(mismatches) 
}

medaka_mismatches_subsets <- medaka_mismatches_subsets %>% filter(!is.na(Barcode))
medaka_aux <- do.call(rbind,vapply(medaka_mismatches_subsets[,1],strsplit,split="_"))
medaka_mismatches_subsets <- medaka_mismatches_subsets %>% mutate(Barcode = medaka_aux[,1], Nreads=as.numeric(medaka_aux[,2]), repetition=as.numeric(medaka_aux[,3]))
medaka_mismatches_subsets <- medaka_mismatches_subsets  %>% mutate(Method="Medaka", homopolymers=NA, mutation=paste0(base.original,Position,base.mutated)) %>% dplyr::rename(repetition=repetitions) 

nanopolish_mismatches_subsets = data.frame(barcode=NA,nseqs=NA,repetition=NA,position=NA,base.original=NA,base.mutated=NA)
results_nanopolish <- system2(paste0("awk \'/^[^#]/ \' BCtopten/nanopolish_consensus/*"), intern = T)
results_nanopolish <- do.call(rbind,vapply(results_nanopolish, strsplit, split="\t"))
results_nanopolish <- results_nanopolish[,c(1,3,5,6)]
nanopolish_aux <- do.call(rbind,vapply(results_nanopolish[,1],strsplit,split="_"))
results_nanopolish <- cbind(nanopolish_aux,results_nanopolish[,-1])
colnames(results_nanopolish) <- c("Barcode","Nreads","repetition", "Position","base.original","base.mutated" )
rownames(results_nanopolish) <- NULL
nanopolish_mismatches_subsets <- data.frame(Barcode = results_nanopolish[,1], 
                                 Nreads=as.numeric(results_nanopolish[,2]),
                                 repetition = as.numeric(results_nanopolish[,3]),
                                 Position = as.numeric(results_nanopolish[,4]), 
                                 base.original=results_nanopolish[,5],
                                 base.mutated = results_nanopolish[,6])

#remove insertions
ind_insertions <- nchar(nanopolish_mismatches_subsets$base.mutated)>1
nanopolish_mismatches_subsets$base.mutated[ind_insertions] <- vapply(nanopolish_mismatches_subsets$base.mutated[ind_insertions] , substr,1,1)
nanopolish_mismatches_subsets <- nanopolish_mismatches_subsets %>% filter(base.mutated!=base.original)

#Split deletions
index_aux <- which(nchar(nanopolish_mismatches_subsets$base.original)==2)
nanopolish_mismatches_subsets$base.original[index_aux] <- substr(nanopolish_mismatches_subsets$base.original[index_aux],2,2)
nanopolish_mismatches_subsets$base.mutated[index_aux] <- "-"
nanopolish_mismatches_subsets$Position[index_aux] <- as.numeric(nanopolish_mismatches_subsets$Position[index_aux])+1

index_aux_3 <- which(nchar(nanopolish_mismatches_subsets$base.original)==3)
aux_df_1 <- nanopolish_mismatches_subsets[index_aux_3,]
aux_df_1$Position <- as.numeric(aux_df_1$Position)+1
aux_df_1$base.original <- substr(aux_df_1$base.original, 2,2)
aux_df_1$base.mutated <- "-"

aux_df_2 <- nanopolish_mismatches_subsets[index_aux_3,]
aux_df_2$Position <- as.numeric(aux_df_2$Position)+2
aux_df_2$base.original <- substr(aux_df_2$base.original, 3,3)
aux_df_2$base.mutated <- "-"

nanopolish_mismatches_subsets <- nanopolish_mismatches_subsets[-index_aux_3,]
nanopolish_mismatches_subsets <- nanopolish_mismatches_subsets %>% rbind(aux_df_1) %>% rbind(aux_df_2)
rm(aux_df_1,aux_df_2)

index_aux_4 <- which(nchar(nanopolish_mismatches_subsets$base.original)==4)
aux_df_1 <- nanopolish_mismatches_subsets[index_aux_4,]
aux_df_1$Position <- as.numeric(aux_df_1$Position)+1
aux_df_1$base.original <- substr(aux_df_1$base.original, 2,2)
aux_df_1$base.mutated <- "-"

aux_df_2 <- nanopolish_mismatches_subsets[index_aux_4,]
aux_df_2$Position <- as.numeric(aux_df_2$Position)+2
aux_df_2$base.original <- substr(aux_df_2$base.original, 3,3)
aux_df_2$base.mutated <- "-"

aux_df_3 <- nanopolish_mismatches_subsets[index_aux_4,]
aux_df_3$Position <- as.numeric(aux_df_3$Position)+3
aux_df_3$base.original <- substr(aux_df_3$base.original, 4,4)
aux_df_3$base.mutated <- "-"

nanopolish_mismatches_subsets <- nanopolish_mismatches_subsets[-index_aux_3,]
nanopolish_mismatches_subsets <- nanopolish_mismatches_subsets %>% rbind(aux_df_1) %>% rbind(aux_df_2)%>% rbind(aux_df_3) %>% arrange(Barcode,Nreads,repetition,Position) %>% mutate(Method="Nanopolish",homopolymers=NA, mutation=paste0(base.original,Position,base.mutated))


single_mismatches <- read.table(paste0(path_analysis_lib,"ConsensusSINGLE_toptenBC.txt"), header=F,row.names=NULL, skip=1)
colnames(single_mismatches)<-c("Barcode","Nreads","repetition","mutation","Method")
single_mismatches <- single_mismatches %>%
  mutate(base.original = str_sub(mutation,1,1), 
         base.mutated=str_sub(mutation,-1,-1),
         Position=str_sub(mutation,start=2,end=-2)) %>%
  mutate(Position = as.numeric(Position), homopolymers="original") #%>% select(-mutation)


write.table(rbind(nanopolish_mismatches_subsets %>% mutate(Method="Nanopolish"),
                  medaka_mismatches_subsets %>% mutate(Method="Medaka"),
                  single_mismatches, file=file_consensus_top10_inLIB))
            
toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

Compute true mutants from all reads according to each method
```{r, eval=FALSE}
tic <- proc.time()

KTLib_topBC_true_mutations <- read.table("3_PreprocessedData/Fig3_LargeSet_Consensus.txt", header=T) %>%
    dplyr::rename(Method=method) %>%
    mutate(Method = recode(Method,!!!methods_capital))%>% 
    filter(Barcode %in% KTLib_topBC$BCsequence) %>%
    filter(Method!="Racon") %>%
    filter((Method=="SINGLe" & homopolymers=="sorted_in_VCS") |
               is.na(homopolymers) | 
               (Method!="SINGLe" & homopolymers=="original") )%>%
    mutate(mutation=paste0(base.original,Position,base.mutated)) %>%
    select(Barcode,Method,mutation)
rownames(KTLib_topBC_true_mutations) <- NULL

#manually replace equivalent mutations
KTLib_topBC_true_mutations <- KTLib_topBC_true_mutations %>%
    mutate(mutation=ifelse(mutation=="G677-","G676-", mutation)) %>%
    mutate(mutation=ifelse(mutation=="G419-","G416-", mutation)) %>%
    mutate(mutation=ifelse(mutation=="G63-","G62-", mutation)) %>%
    mutate(mutation=ifelse(mutation=="G429-","G428-", mutation)) %>%
    mutate(mutation=ifelse(mutation=="A1545-","A1543-", mutation)) %>%
    mutate(mutation=ifelse(mutation=="G244-","G243-", mutation)) %>%
    mutate(mutation=ifelse(mutation=="G263-","G261-", mutation)) %>%
    mutate(mutation=ifelse(mutation=="G677-","G676-", mutation))

aux_nanopolish <- data.frame(Barcode= KTLib_topBC_true_mutations$Barcode[KTLib_topBC_true_mutations$mutation=="GGAC1306G"], 
                             Method="Nanopolish", 
                             mutation=c("G1307-","A1308-", "C1309-"))
    
KTLib_topBC_true_mutations <- KTLib_topBC_true_mutations %>% 
    filter(mutation!="GGAC1306G") %>%
    rbind(aux_nanopolish)

compare_mutations <- reshape2::dcast(KTLib_topBC_true_mutations %>%
                                         select(Barcode, Method, mutation),
                                     Barcode+mutation~Method)

KTLib_topBC_true_mutants <- KTLib_topBC_true_mutations %>%
    group_by(Barcode,Method) %>% 
    summarise(n_mutations=n(),mutant = paste0(mutation, collapse = " ")) %>%
    ungroup()

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

Compute success rate for consensus calling by method, as a function of reads used
```{r LIB10BCs success rate, eval=FALSE}
tic <- proc.time()
# consensus_by_method
KTLib_topBC_subsets_mutations <- read.table("3_PreprocessedData/Fig3_ConvergenceConsensus_top10BC.txt", header=T) %>%
    filter(Method!="racon") %>%
    filter((Method=="single" & homopolymers=="sorted_for_VCS") |
               is.na(homopolymers) | 
               (Method!="single" & homopolymers=="original") )  %>%
    select(-homopolymers,-variant) %>%
    select(-base.original,-Position,-base.mutated)%>%
    mutate(Method = recode(Method,!!!methods_capital))

KTLib_topBC_subsets_mutations <- KTLib_topBC_subsets_mutations %>%
    mutate(mutation=ifelse(mutation=="G677-","G676-", mutation)) %>%
    mutate(mutation=ifelse(mutation=="G419-","G416-", mutation)) %>%
    mutate(mutation=ifelse(mutation=="G63-","G62-", mutation)) %>%
    mutate(mutation=ifelse(mutation=="G429-","G428-", mutation)) %>%
    mutate(mutation=ifelse(mutation=="A1545-","A1543-", mutation)) %>%
    mutate(mutation=ifelse(mutation=="G244-","G243-", mutation)) %>%
    mutate(mutation=ifelse(mutation=="G263-","G261-", mutation)) %>%
    mutate(mutation=ifelse(mutation=="G677-","G676-", mutation))

indx <- which(KTLib_topBC_subsets_mutations$mutation=="GGAC1306G")
aux_nanopolish <- data.frame(Barcode= rep(KTLib_topBC_subsets_mutations$Barcode[indx],3),
                            Nreads = rep(KTLib_topBC_subsets_mutations$Nreads[indx], each=3),
                            repetition = rep(KTLib_topBC_subsets_mutations$repetition[indx], each=3),
                             Method="Nanopolish",
                             mutation=rep(c("G1307-","A1308-", "C1309-"), lenght.out=length(indx)))
    
KTLib_topBC_subsets_mutations <- KTLib_topBC_subsets_mutations %>% 
    filter(mutation!="GGAC1306G") %>%
    rbind(aux_nanopolish)

KTLib_topBC_subsets_mutant <- KTLib_topBC_subsets_mutations %>%
    select(Barcode,Nreads,repetition,Method,mutation) %>%
    group_by(Barcode,Nreads,repetition,Method) %>%
    summarise(n_mutations=n(),mutant=paste(mutation, collapse = " "))

KTLib_topBC_subsets_success_rate <- left_join(KTLib_topBC_subsets_mutant, 
                                         KTLib_topBC_true_mutants,
                                         by=c("Barcode","Method"),
                                         suffix=c(".subset",".true")) %>%
    mutate(equal = mutant.subset==mutant.true ) %>% 
    filter(!is.na(n_mutations.true)) %>%
    group_by(Barcode,Nreads,Method) %>% 
    dplyr::summarise(success_rate= sum(equal)/n()) 

ggplot(KTLib_topBC_subsets_success_rate,aes(x=Nreads,y=success_rate,col=Method) )+geom_point()+geom_line()+
  facet_wrap(~Barcode,ncol=2)+scale_color_manual(values=colors_vector_capital[-6])+
  ylab("Success rate")  

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```



### Figure 3A
```{r Fig3A, fig.width=3, fig.height=3, eval=FALSE}
fig_3a <- ggplot(KTLib_topBC_subsets_success_rate %>% 
                    filter((Method =="SINGLe"|Method=="Medaka") & Barcode=="TCAGTTATTACGTCCTGCGGACAGTAAGGTCCGAG"),
                aes(x=Nreads,y=success_rate,col=Method) )+
    geom_point()+
    geom_line()+
    scale_color_manual(values=colors_vector_capital[1:2])+
    scale_y_continuous(limits=c(0,1),breaks = c(0,.25,.5,.75,1))+
    geom_hline(aes(yintercept=0.9), col=grey(0.5),lty="dashed")+
    theme_plot+
  theme(legend.position = c(.9,0.4),
        legend.justification = c("right","top"),
        legend.title = element_blank(),
        legend.spacing.y = unit(0,"line"))+
    labs(tag = "A",x="Reads", y="Success rate")

fig_3a
```

### Figure S7: top 10 bc convergence

```{r FigS7, fig.width=7, fig.height=10, eval=F}
Fig_Sup7<- ggplot(KTLib_topBC_subsets_success_rate %>% filter(Method!= "Racon"),
                  aes(x=Nreads,y=success_rate,col=Method) )+
    geom_point()+geom_line()+
  facet_wrap(~Barcode,ncol=2)+scale_color_manual(values=colors_vector_capital[-6])+
  scale_y_continuous(breaks=c(0,.5,1))+
  geom_hline(aes(yintercept=0.9), col=grey(.5),linetype="dashed")+
  ylab("Success rate")  +theme_plot+
  theme(plot.margin = unit(c(0,5.5,5.5,5.5),"points"),
        strip.background = element_rect(fill= "white"),
        strip.text = element_text(size=8), 
        legend.position = "bottom")

Fig_Sup7
ggsave(Fig_Sup7, file=paste0(path_save_figures,"FigureS7.png"), 
       dpi = 'print', width = 17, height=18,units="cm")

```


### Load consensus and count reads per barcode
```{r, eval=F}
tic <- proc.time()
KTLib_mutations <- read.table("3_PreprocessedData/Fig3_LargeSet_Consensus.txt") %>% 
    filter(method!="Racon") %>% 
    filter((method=="single" & homopolymers=="sorted_in_VCS") |
               is.na(homopolymers) | 
               (method!="single" & homopolymers=="original") )

aux_nanopolish <- KTLib_mutations %>%
        filter(nchar(base.original)>1) %>%
    mutate(base.mutated="-")
aux_nanopolish <- (aux_nanopolish %>% mutate(base.original=str_sub(base.original, 2,2)) ) %>%
    rbind(aux_nanopolish%>% mutate(base.original=str_sub(base.original, 3,3))) %>% 
    rbind(aux_nanopolish%>% mutate(base.original=str_sub(base.original, 4,4)))

KTLib_mutations <- KTLib_mutations %>%
    filter(!nchar(base.original)>1) %>%
    rbind(aux_nanopolish)



KTLib_nreads_per_barcode <- barcodes_4reads %>% 
    group_by(BCsequence) %>% dplyr::summarise(Nreads=n()) %>% 
    mutate(Barcode = BCsequence) %>% select(-BCsequence) %>% select(Barcode, Nreads) %>% 
    arrange(Barcode)

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

### Compute the actual number of mutations in the library by SINGLe

```{r, fig.width=3, fig.height=3, eval=F}
tic <- proc.time()
n_mutations <- KTLib_mutations %>%
    filter(method=="single") %>%
    group_by(Barcode) %>%
    dplyr::summarise(n_mutations = n()) %>% 
    left_join(KTLib_nreads_per_barcode) %>% filter(Nreads>5)

plot_nmutations <- ggplot(n_mutations, aes(x=n_mutations))+
    geom_histogram()+
    theme_plot+
    labs(x="Number of mutations", y="Counts")
mean_n_mutations <- mean(n_mutations$n_mutations)
median_n_mutations <- median(n_mutations$n_mutations)
sd_n_mutations <- sd(n_mutations$n_mutations)

mutation_rate_measured <- mean_n_mutations/1662*1000
 sd_n_mutations/1662*1000
 
 
plot_nmutations
cat("Mean N mutations: ", round(mean_n_mutations,2), 
 "\nMedian N mutations: ", round(median_n_mutations,2),
 "\nSD N mutations: ", round(sd_n_mutations,2), 
 "\nMutation rate measured: ", round(mutation_rate_measured,2), "muts/kb")

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

## Compare number of mutations by SINGLe and Medaka
```{r, eval=F}
tic <- proc.time()
consensus_mismatches_ms <- KTLib_mutations %>% 
    filter((method=="Medaka" )| (method=="single")) 

nmismatches_ms <- consensus_mismatches_ms %>%
  group_by(Barcode, method) %>%
  summarise(n_mismatches=n())



nmismatches_ms_difference <- left_join(nmismatches_ms %>% filter(method=="Medaka" ), 
                                       nmismatches_ms %>% filter(method=="single" ), 
                                       by="Barcode", suffix=c("",".single")) %>%
  mutate(difference = n_mismatches.single-n_mismatches) %>%
  left_join(KTLib_nreads_per_barcode, by = "Barcode") %>% 
  ungroup() %>% 
  mutate(Reads =  ifelse(Nreads<=5,"5 or less reads", ifelse(Nreads<=10,"6 to 10 reads","10 or more reads"))) %>%
  mutate(Reads=factor(Reads, levels=c("5 or less reads","6 to 10 reads","10 or more reads")))%>%
  select(-Nreads) %>%  ungroup()  %>% filter(!is.na(n_mismatches.single))


nmismatches_ms_difference_hist <- nmismatches_ms_difference %>% 
    group_by(method,difference,Reads) %>% 
    summarise(counts=n()) 

nmismatches_ms_difference_hist <- nmismatches_ms_difference %>% 
    ungroup()%>%
    mutate(difference = ifelse(difference< -6, -6, difference))%>%
    group_by(method,difference,Reads) %>% 
    summarise(counts=n()) %>%
     group_by(Reads) %>% mutate(perc = counts/sum(counts)*100) 

table_mutations_compare <- nmismatches_ms_difference_hist %>% 
    mutate(diff = ifelse(difference<0,-1,ifelse(difference==0,0,1)))%>% 
    group_by(Reads, diff) %>% summarise(n=sum(counts)) %>%
    group_by(Reads) %>% 
    mutate(perc=n/sum(n)*100)
toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

### Figure 3B
```{r Fig 3B, fig.width=5, fig.height=5, eval=F}
fig_3b <-  ggplot(table_mutations_compare , aes(x=diff,y=perc))+
    geom_bar(position = position_dodge2(preserve="single"),stat="identity")+
    xlab("Method that reports more mutations")+ylab("Percentage")+
    facet_wrap(~Reads, ncol=1, strip.position = "right", labeller = label_wrap_gen(width=12))+
    scale_y_continuous(limits = c(0,110), breaks=c(50,100))+
    scale_x_continuous(breaks=seq(-1,1,by=1), labels=c("Medaka", "Both match", "SINGLe"))+
    theme(legend.position = c(0,1),
            legend.justification = c("left","top"),
            legend.title = element_text(size = rel(1)),
            legend.text = element_text(size=rel(1)),
            legend.spacing.y = unit(0.5,"line"),
            legend.key.size = unit(2,"pt"),
            plot.margin = unit(c(0,5.5,5.5,5.5),"points"), 
          strip.background = element_rect(fill=grey(.95)))+
    labs(tag = "B")+
    geom_text(aes(y=perc,label=paste0(round(perc),"%"),vjust=-.2))

fig_3b
```


## Distribution of mutations and skewness 
```{r skewness, eval=F}
tic <- proc.time()
results_for_entropy<- KTLib_mutations %>%  
  filter(base.original!="wildtype" ) %>%
  filter(Position >0 & Position <=1662) %>%
  left_join(KTLib_nreads_per_barcode, by="Barcode") %>%
  mutate(group_Nreads_barcode = ifelse(Nreads>10, "> 10 reads", 
                                      ifelse(Nreads>5,"6 to 10 reads","5 or less reads"))) %>%
  mutate(group_Nreads_barcode=factor(group_Nreads_barcode, levels=c("5 or less reads","6 to 10 reads", "> 10 reads")))%>%
  group_by(method,Position,group_Nreads_barcode, homopolymers) %>%
  dplyr::summarise(n=n())%>% ungroup() %>%
    mutate(Method = recode(method,!!!setNames(c("Nanopolish","Medaka","SINGLe","Guppy","No weights"),c("nanopolish","medaka","single","nanopore","no_weights")))) %>%
    select(-method)


skewness <- results_for_entropy%>%
  group_by(Method, group_Nreads_barcode) %>%
  dplyr::summarise(skewness=skewness(n))

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

### Figure 3C
```{r Fig3C, fig.width=5, fig.height=5, eval=F}
fig_3c <- ggplot(results_for_entropy %>% filter(Method=="Medaka" | Method == "SINGLe"),
                 aes(x=Position, y=n, col=Method))+
    geom_point(size=0.5)+
    scale_color_manual(values=colors_vector_capital)+
    ylab("Counts")+
    facet_grid(Method~group_Nreads_barcode)+
    geom_text(
        data    = skewness  %>% filter(Method=="Medaka" | Method == "SINGLe"),
        mapping = aes(x = -Inf, y = -Inf, label = paste0("sk=", round(skewness,2))),
        hjust   = 0,
        vjust   = 1, 
        y=rep(140,6), x=rep(80,6), 
        col="black"
    )+
    theme_plot+
    theme(legend.position = "none",plot.margin = unit(c(1,5.5,5.5,5.5),"points"),)+
    labs(tag = "C")

fig_3c
```

### Figure S8

```{r FigS8, eval=F}
Figure_S8 <- ggplot(results_for_entropy ,
                 aes(x=Position, y=n, col=Method))+
  geom_point(size=0.5)+scale_color_manual(values=colors_vector_capital[-6])+geom_line()+
  ylab("Counts")+
  facet_grid(Method~group_Nreads_barcode)+
   geom_text(
  data    = skewness,
  mapping = aes(x = -Inf, y = -Inf, label = paste0("sk=", round(skewness,2))),
  hjust   = 0,
  vjust   = 1, 
  y=rep(200,nrow(skewness)), x=rep(80,nrow(skewness)), 
  col="black"
)+theme_plot+  theme(legend.position = "none",plot.margin = unit(c(0,5.5,5.5,5.5),"points"),strip.background = element_rect(fill="white"),strip.text = element_text(size=8))

Figure_S8
ggsave(Figure_S8, filename="6_Figures/FigureS8.png",dpi = 'print', width = 17, units = "cm", height=20)

```

## Figure 3

```{r Fig3, eval=F}
Figure3 <- ( fig_3a | fig_3b) / fig_3c 
Figure3
ggsave(Figure3, file=paste0(path_save_figures,"Figure3.png"),
       dpi = 'print', width = 16, height=16,units="cm")

```






## How is the mutation rate predicted compared to the one reported by Agilent?
### Figure S9
```{r FigS9, eval=F}
kit_rates <- data.frame(
    Mutation=c("Deletion", "Ts AG TC" ,"Ts GA CT"  ,"Tv AT TA","Tv AC TG" ,"Tv GC CG", "Tv GT CA"),
    true_perc=c(4.8,17.5,25.5,28.5,4.7,4.1,14.1))
aux <- KTLib_mutations %>% 
    left_join(KTLib_nreads_per_barcode,by="Barcode") %>%
    mutate(reads_number = case_when(
        Nreads<6 ~ "4 or 5 reads", 
        Nreads <11 ~ "6 to 10 reads", 
        Nreads>10 ~ "more than 10 reads") )%>%
    mutate(reads_number= factor(reads_number, levels=c("4 or 5 reads", "6 to 10 reads", "more than 10 reads")))%>%
    mutate(Mutation=case_when(
        base.mutated=="-"~ "Deletion",
        base.original=="A" & base.mutated=="G" ~ "Ts AG TC",
        base.original=="T" & base.mutated=="C" ~ "Ts AG TC",
        base.original=="G" & base.mutated=="A" ~ "Ts GA CT",
        base.original=="C" & base.mutated=="T" ~ "Ts GA CT",
        base.original=="A" & base.mutated=="T" ~ "Tv AT TA",
        base.original=="T" & base.mutated=="A" ~ "Tv AT TA",
        base.original=="A" & base.mutated=="C" ~ "Tv AC TG",
        base.original=="T" & base.mutated=="G" ~ "Tv AC TG",
        base.original=="G" & base.mutated=="C" ~ "Tv GC CG",
        base.original=="C" & base.mutated=="G" ~ "Tv GC CG",
        base.original=="G" & base.mutated=="T" ~ "Tv GT CA",
        base.original=="C" & base.mutated=="A" ~ "Tv GT CA"
    )) %>%
    group_by(method,reads_number,Mutation)%>%
    summarise(counts = n()) %>%
    group_by(method,reads_number)  %>%
    mutate(counts=counts/sum(counts)*100) %>%
    left_join(kit_rates, by="Mutation")%>%
    ungroup()%>%
    mutate(Method = recode(method,!!!setNames(c("Nanopolish","Medaka","SINGLe","Guppy","No weights"),c("nanopolish","medaka","single","nanopore","no_weights")))) %>%
    select(-method)

plot_mr_xy <- ggplot(aux, aes(col=Method, x=true_perc, y=counts, shape=Mutation,label=Mutation))+
    geom_point()+
    geom_abline(aes(slope=1,intercept=0), lty="dashed", col=grey(.5))+
    facet_wrap(~reads_number,ncol=1)+xlim(c(0,30))+scale_shape_manual(values=1:7)+theme_plot+
    xlab("By manufacturer (%)")+
    ylab("Observed (%)")+theme_plot+ggtitle("Mutation rate")


plot_mr_cor <- ggplot(aux %>% group_by(Method, reads_number) %>% 
           summarise(cor=cor(counts, true_perc)) %>% 
           ungroup(), aes(x=reads_number, y=cor, col=Method))+
    geom_point()+
    geom_line(aes(x=as.numeric(reads_number)))+
    ylab("Correlation")+xlab("")+
    # ggtitle("Compare observed mutation rate vs. manufacturer's report")+
    scale_y_continuous(breaks=c(0.8,0.9,1), limits = c(0.8,1))+theme_plot+
    theme(axis.text.x = element_text(angle=90,hjust = 1,vjust=0.5))


FigureS9 <- grid.arrange(plot_mr_xy , plot_mr_cor+theme(legend.position = "none"), ncol=2) 
             
ggsave(FigureS9, file=paste0(path_save_figures,"FigureSMutRate.png"), 
       dpi = 'print', width = 16, height=12,units="cm")

```



# How many reads needed to fit?

```{r , eval=F}
tic <- proc.time()
a <- load("/media/rocio/Data10TB/Work/Bibliography/RondelezLab/2022_SINGLe/GigaScience_submission/Espada_et_al/5_Revision/01_Subsets/01_KlenTaq/barcode01_allseqs_countspnq.Rdata")
full_fits <- read.table("/media/rocio/Data10TB/Work/Bibliography/RondelezLab/2022_SINGLe/GigaScience_submission/Espada_et_al/5_Revision/01_Subsets/01_KlenTaq/barcode01_allreads_fits.txt", header=T)

full_fits <- full_fits %>%
    mutate(rqs = 1-deviance/null_deviance)

error_rate <- data_mut %>%
    group_by(pos,strand,nucleotide) %>%
    dplyr::summarise(error_rate = sum(count)/(sum(count)+sum(counts.wt)),
                     counts_mut = sum(count))

full_fits<- full_fits %>%
    left_join(error_rate, by=c("strand","pos","nucleotide"))

## Full fits mutated positions
mutationstested <- kt_true_mutations %>% 
    select(position, nucleotide) %>%
    dplyr::rename(pos=position)
full_fits_mut <- left_join(mutationstested, full_fits)
hex <- scales::hue_pal()(4)
names(hex)<- c(303,331,879,1316)

full_fits_mut_examples <- full_fits_mut %>%
    filter( (pos==1316 & nucleotide=="A" & strand=="+")|
                (pos==879  & nucleotide=="T" & strand=="-")|
                (pos==303  & nucleotide=="A" & strand=="-")|
                (pos==331  & nucleotide=="T" & strand=="+"))

plot_counts_RSQ <- ggplot(full_fits,aes(x=rqs,y=counts_mut))+
    geom_point(col=rgb(0,0,0,0.1), stroke=0, size=2)+
    scale_y_log10()+
    theme_plot+
    geom_point(data=full_fits_mut_examples,
               aes(x=rqs,y=counts_mut, col=as.factor(pos)),
               size=2,shape="square")+
    scale_color_manual(values=hex)+
    labs(x="Quality (Q)", y="Counts(C)")+theme(legend.position = "none")

xlim <- 1:50
plots <- list()
for(i in 1:nrow(full_fits_mut_examples)){
    nuc <- full_fits_mut_examples$nuc[i]
    pot <- full_fits_mut_examples$pos[i]
    str  <- full_fits_mut_examples$strand[i]
    aux_df <- data_mut %>%
        dplyr::filter(pos==pot & nucleotide ==nuc & strand==str)%>%
        dplyr::select(strand,QUAL,count, counts.wt,counts.scaled,counts.wt.scaled) %>%
        dplyr::mutate(tot.scaled=counts.scaled+counts.wt.scaled) %>%
        dplyr::mutate(proportion.wt.scaled=counts.wt.scaled/tot.scaled) %>%
        arrange(QUAL)
    qval  <- aux_df$QUAL
    yvals <- aux_df$proportion.wt.scaled
    no <- !is.na(yvals)
    yvals <-yvals[no]
    qval <- qval[no]

    fitt  <- stats::glm(yvals~ qval,family = "binomial")

    prediction <- data.frame(QUAL=xlim, fitted = predict(fitt,list(qval=xlim), type = "response"))
    # lines(qval,predict(fitt,list(qval=qval),type="response"), col=2, lwd=2)

    qq<- format((full_fits_mut_examples %>%
                     filter(pos==pot & nucleotide ==nuc & strand==str))$rqs,scientific=F, digits=2)
    cc<- (full_fits_mut_examples %>%
              filter(pos==pot & nucleotide ==nuc & strand==str))$counts

    plots[[i]] <- ggplot(aux_df, aes(x=QUAL,y=proportion.wt.scaled))+
        geom_point()+
   
        geom_line(data=prediction, aes(x=QUAL, y=fitted),
                  col=hex[as.character(pot)])+
        ggtitle(paste(pot,nuc,str),
                subtitle=paste("C:", cc, " QF:",qq))+
        labs(x="Qscore", y="ncorrect")+
        xlim(range(xlim))+
        scale_y_continuous(breaks=c(0,1), limits = c(0,1))+
        theme(plot.title = element_text(margin = margin(0,0,0,0)),
              plot.subtitle =  element_text(margin = margin(0,0,0,0)))

}


toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

### Figure S10
```{r FigS10, eval=F}
FigureS10 <- grid.arrange(plot_counts_RSQ , plots[[1]], plots[[4]], plots[[3]], plots[[2]] ,
                       layout_matrix=matrix(c(1,1,2,3,4,5), byrow = F,nrow=2),
                       widths=c(2.5,1,1))
ggsave(FigureS10, file=paste0(path_save_figures,"FigureExtra.png"), dpi = 'print',
       height = 8.7, width=20,units="cm")

FigureS10
```

# Correlation neighbouring Qscore
```{r Qscore correlation, eval=F}
tic <- proc.time()
path_11 <- "/media/rocio/Data10TB/minIONreads/KlenTaq/2019_11_KT/KT7_fullpipeline/KT7/BasecallSuperior/pass/out/Demultiplex/minimap2/SINGLE/"

data_11 <- read.table(paste0(path_11, "barcode05.for_corrected.txt"), header=T)
data_11 <- data_11
error_lines <- which(!data_11$isWT&
                         !(data_11$position==867 & data_11$nucleotide=="G")&
                         !(data_11$position==1120 & data_11$nucleotide=="A") &
                         !(data_11$position==1214 & data_11$nucleotide=="T") &
                         !(data_11$position==1316 & data_11$nucleotide=="A") &
                         !(data_11$nucleotide=="-"))
correct_lines <- which(data_11$isWT &
                       !(data_11$position==867 & data_11$nucleotide=="G")&
                       !(data_11$position==1120 & data_11$nucleotide=="A") &
                       !(data_11$position==1214 & data_11$nucleotide=="T") &
                       !(data_11$position==1316 & data_11$nucleotide=="A") &
                       !(data_11$nucleotide=="-"))


rand_pos <- sample(1:nrow(data_11), 200)
df11 <- rbind(
    data.frame(type="Errors",
               difference = abs(data_11$original_quality[error_lines]-data_11$original_quality[error_lines+1]), wt = data_11$isWT[error_lines+1]),
    data.frame(type="Correct",
               difference = abs(data_11$original_quality[correct_lines]-data_11$original_quality[correct_lines+1]), wt = data_11$isWT[correct_lines+1])) %>%
    rbind( data.frame(type="Random",
               difference = abs(data_11$original_quality[rand_pos]-data_11$original_quality[sample(1:nrow(data_11), 200)]), wt = data_11$isWT[rand_pos])
    ) %>%
    mutate(difference = ifelse(difference>20, 20, difference))

toc <- proc.time()
cat("it took" , (toc-tic)[3], "s to run this chunk\n"); rm(tic, toc)

```

### Figure S3
```{r FigS3, fig.width=7, fig.height=3, eval=F}
FigureS3 <- ggplot(df11, aes(x=difference))+
    geom_histogram(binwidth = 1,aes(y=..density.., text=..density..))+
    xlab("Qscore difference with next nucleotide")+
    ylab("Counts")+
    labs(tag="C")+
    facet_wrap(~type, scales = "free")

ggsave(FigureS3, file=paste0(path_save_figures,"FigureS3.png"), dpi = 'print',
       height = 5, width=15,units="cm")

FigureS3
```



Depricated nextpolish:
```{r, eval=F}
# create_run_cfg <- function(run_cfg_file, ref_fasta, workdir){
#     cat("[General]\n", file=run_cfg_file, append=FALSE)
#     cat('job_type = local\n', file=run_cfg_file, append=TRUE)
#     cat('job_prefix = nextPolish\n', file=run_cfg_file, append=TRUE)
#     cat('task = best\n', file=run_cfg_file, append=TRUE)
#     cat('rewrite = yes\n', file=run_cfg_file, append=TRUE)
#     cat('rerun = 3\n', file=run_cfg_file, append=TRUE)
#     cat('parallel_jobs = 2\n', file=run_cfg_file, append=TRUE)
#     cat('multithread_jobs = 4\n', file=run_cfg_file, append=TRUE)
#     cat('genome =', ref_fasta,' #genome file\n', file=run_cfg_file, append=TRUE)
#     cat('genome_size = auto\n', file=run_cfg_file, append=TRUE)
#     cat('workdir =',workdir, "\n",file=run_cfg_file, append=TRUE, sep="")
#     cat('   polish_options = -p {multithread_jobs}\n', file=run_cfg_file, append=TRUE)
# 
#     cat('[lgs_option]\n', file=run_cfg_file, append=TRUE)
#     cat('lgs_fofn =lgs.fofn\n', file=run_cfg_file, append=TRUE)
#     cat('lgs_options = -min_read_len 1k -max_depth 100\n', file=run_cfg_file, append=TRUE)
#     cat('lgs_minimap2_options = -x map-ont\n', file=run_cfg_file, append=TRUE)
#     return(0)
# }
# 
# 
# path0 <- "/home/respada/22nextpolish/bc06/"
# setwd(path0)
# run_folder  <- "RunFolder/"
# save_folder <- ""
# data_folder <- ""
# 
# for(bc in kt7_barcodes){
# 
#     #Load fastq file with all reads
#     all_sequences <- readLines(paste0(path_analysis_lib7, bc,".fastq"))
#     name_lines <- grep(pattern="^@",x=all_sequences)
#     save_file <- paste0(save_folder,bc, "_ConsensusByNextPolish.fasta")
# 
#     for(ns in subsets){
#         cat("\t",ns,"\n")
#         if(ns> length(name_lines)){next()}
#         for(r in 1:n_repetitions){
#             cat("----------------",ns, r, "\n")
# 
#             #Create file with subset of sequences
#             subset_name = paste0(bc,"_s",ns, "_r",r)
#             filename_subset <- paste0(run_folder,subset_name,".fastq")
#             choose_lines <- sample(name_lines, size= ns)
#             choose_lines_all <- sort(c(choose_lines,choose_lines+1,choose_lines+2,choose_lines+3))
#             writeLines(text = all_sequences[choose_lines_all], con = file(filename_subset))
# 
#             ## Run nextpolish
#             system2(paste0("echo '", subset_name,".fastq' >", run_folder, "lgs.fofn"))
#             create_run_cfg(run_cfg_file=paste0(run_folder, "run.cfg"),
#                            ref_fasta=reference_file,
#                            workdir=paste0(subset_name,"/"))
#             system2(paste0("cp ", reference_file, " ",run_folder, reference_file ))
#             system2(paste0("cd ",run_folder,"; nextPolish run.cfg"))
#             system2(paste0("sed -i '1s/.*/>",subset_name,"'/ ", run_folder, subset_name, "/genome.nextpolish.fasta "))
#             system2(paste0("cat ",run_folder, subset_name, "/genome.nextpolish.fasta >>", save_file))
#             system2(paste0("rm -r ", run_folder, subset_name))
#             system2(paste0("cd ",run_folder,";rm lgs.fofn  run.cfg"))
#             
#         }
#     }
# }
```

Alternative to fill deletions faster:
```{r, eval=F}


### OPTION
t0 <- proc.time()
a <- str_locate_all(reads_aligned, "-+")
for(i in seq(a)){
    if(nrow(a[[i]])==0){next()}
    pos_to_replace <- IRanges(a[[i]][,"start"], a[[i]][,"end"])
    pos_to_average_start <- IRanges(a[[i]][,"start"]-1,a[[i]][,"start"]-1)
    pos_to_average_end <- IRanges(a[[i]][,"end"]+1,a[[i]][,"end"]+1)

    if(start(pos_to_average_start)[1]==0){
            start(pos_to_average_start)[1] <- end(pos_to_average_start)[1] <- start(pos_to_average_end)[1] 
    }
    n <- length(pos_to_average_start)
    pos_max <- width(reads_aligned)[i]
    if(start(pos_to_average_end)[n]>pos_max){
        start(pos_to_average_end)[n] <-  start(pos_to_average_start)[n]      
        end(pos_to_average_end)[n] <- start(pos_to_average_start)[n]   
    }
    
    before <- extractAt(scores_aligned[[i]], at = pos_to_average_start)
    before <- unlist(as(PhredQuality(before),"NumericList"))
    after <- extractAt(scores_aligned[[i]], at = pos_to_average_end)
    after <- unlist(as(PhredQuality(after),"NumericList"))
    both <- cbind(before,after)
    if(gaps_weights=="mean"){
        replace_qscore <- apply(both,1, mean)
    }else if(gaps_weights=="minimum"){
        replace_qscore <- apply(both,1, min)
    }
    replace_qscore <- lapply(replace_qscore, function(x){as(x,"PhredQuality")})
    replace_qscore_v <- vapply(seq(replace_qscore), function(x){paste0(rep(replace_qscore[[x]], times=width(pos_to_replace)[x]), collapse="")})
    # print(scores_aligned[[i]])
    scores_aligned[i] <- replaceAt(scores_aligned[i],at = pos_to_replace, value = replace_qscore_v)
    # print(scores_aligned[[i]])
}
t1 <- proc.time()

t1-t0
  


### ORIGINAL
if(verbose){message("Assign values to deletions\n")}
if(verbose){pb=utils::txtProgressBar(min =0,max=length(reads_aligned),style = 3)}
t2 <- proc.time()

for(i in seq(length(reads_aligned))){
    # if(verbose){utils::setTxtProgressBar(pb,i)}
    ## REPLACE DELETION SCORES
    del_positions <- BiocGenerics::start(Biostrings::vmatchPattern("-",reads_aligned[i])[[1]])
    ndel <- length(del_positions)
    
    if(ndel==0){next()}
    nr <- width(reads_aligned)[i]
    
    #If they are at the beginning, I replace by the first Qscore
    if(del_positions[1]==1){
        n_del_start=1
        if(length(del_positions)==1){
            n_del_start <- 1
        }else{
            condition = del_positions[n_del_start+1]==del_positions[n_del_start]+1
            while(condition){
                n_del_start <- n_del_start+1
                condition1 = is.na(del_positions[n_del_start+1])
                if(condition1){
                    condition=FALSE
                }else{
                    condition = del_positions[n_del_start+1]==del_positions[n_del_start]+1
                }
            }
        }
        scores_aligned[i] <- Biostrings::replaceAt(scores_aligned[i],at=IRanges(1,n_del_start),as.character(subseq(scores_aligned[i], 1, n_del_start)))
    }else{n_del_start=NA}
    
    #If they are at the end, I replace by the last Qscore
    if(del_positions[ndel]==nr){
        #detect which are the values on the end of the string
        Breaks <- c(0, which(diff(del_positions) != 1), length(del_positions))
        Breaks_ini <- Breaks[length(Breaks)-1]+1
        Breaks_end <- Breaks[length(Breaks)]
        del_pos_ini <- del_positions[Breaks_ini]
        del_pos_end <- del_positions[Breaks_end]
        replacement <- paste0(rep(as.character(subseq(scores_aligned[i], del_pos_ini-1, del_pos_ini-1)), Breaks_end-Breaks_ini+1),collapse="")
        scores_aligned[i] <-Biostrings::replaceAt(scores_aligned[i],
                                                  at=IRanges(del_pos_ini,del_pos_end),
                                                  replacement)
        n_del_end <- length(Breaks_ini:Breaks_end)
    }else{n_del_end=NA}
    
    if(!is.na(n_del_end)){del_positions <- del_positions[-seq(ndel-n_del_end+1,ndel)]}
    if(!is.na(n_del_start)){del_positions <- del_positions[-seq(n_del_start)]}
    if(length(del_positions)==0){next()}
    ndel <- length(del_positions)
    for(j in seq_along(del_positions)){
        jmin <- j
        while(j < ndel && del_positions[j+1]==del_positions[j]+1){ j <- j+1 }
        jmax<-j
        
        q_upstr <- subseq(scores_aligned[i], del_positions[jmin]-1, del_positions[jmin]-1)
        q_dwstr <- subseq(scores_aligned[i], del_positions[jmax]+1, del_positions[jmax]+1)
        
        if(gaps_weights=="mean"){
            prob_upstr <- as(q_upstr,"NumericList")[[1]]
            prob_dwstr  <- as(q_dwstr,"NumericList")[[1]]
            prob <- mean(prob_dwstr,prob_upstr)
        }else if(gaps_weights=="minimum"){
            prob_upstr <- as(q_upstr,"NumericList")[[1]]
            prob_dwstr  <- as(q_dwstr,"NumericList")[[1]]
            prob <- min(prob_dwstr,prob_upstr)
        }
        
        qscore_new <- as(prob,"PhredQuality")
        qscore_new_vec <- paste0(rep(as.character(qscore_new), jmax-jmin+1), collapse="")
        scores_aligned[i] <- Biostrings::replaceAt(scores_aligned[i],at=IRanges(del_positions[jmin],del_positions[jmax]),qscore_new_vec)
    }
    
}

t3 <- proc.time()

t3-t2



```



# Session info {.unnumbered}

```{r sessionInfo, echo=FALSE}
sessionInfo()
```
