---
title: "eval_duplex_reads"
author: "bhaas"
date: '2024-10-13'
output: github_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
library(tidyverse)
```

The SGNEx ONT reads were each aligned to themselves by their reverse-complement. Best local alignments were obtained by Smith-Waterman alignment as per: 
https://github.com/TrinityCTAT/CTAT-LR-fusion/blob/devel/util/misc/find_duplex_reads.py

limited to reads at least 1kb in length, and for which the 100 bases at each termini were first removed before aligning to avoid terminal polyA or other sequences that would tend to show up as being soft-clipped in a genome alignment.

The fraction of all reads with at least a frac of total such reads aliging to its rev-comp at 0.1, 0.2, ..., 0.8 of the read length were assessed and are plotted below.


```{r}

data = read.csv("data/summary.tsv", header=T, sep="\t", check.names=F)

data %>% head()
```


```{r}

data$sample = str_replace(data$sample, "-EV_", "_")


```


```{r}

data$seqtype = sapply(data$sample, function(x) { str_split(x, "_")[[1]][3]})

data %>% head()


```


```{r}

data = data %>% gather(key=frac_rc_align, value=frac_of_reads, `0.1`, `0.2`, `0.3`, `0.4`, `0.5`, `0.6`, `0.7`, `0.8`)

data %>% head()
```

```{r}

data$sample = str_replace(data$sample, "SGNex_", "")
data$sample = str_replace(data$sample, ".duplex_read_report.tsv.gz", "")

```


```{r}

duplex_frac_plot = data %>% ggplot(aes(x=sample, y=frac_of_reads, fill=frac_rc_align)) + 
    theme_bw() +
    geom_bar(stat='identity', position='dodge') +
    facet_wrap(~seqtype, scale='free_x') +
    theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
    ggtitle("Fraction of reads self-aligning as duplex type")

duplex_frac_plot

```

```{r}

ggsave(duplex_frac_plot, file="SGNEx_duplex_frac_plot.svg", width=8, height=8)


```





