figure_2b.Rmd

---
title: "Fig.2B. Influence of repeats on mtDNA deletions distributon"
output:
  workflowr::wflow_html:
    toc: true
    toc_float: yes
    theme: journal
    highlight: textmate
    code_folding: hide
    df_print: paged
---

```{r setup, echo=FALSE, include=FALSE}
knitr::opts_chunk$set(
  autodep        = TRUE,
  cache          = FALSE,
  cache.lazy     = FALSE,
  dev            = c("png", "pdf"),
  echo           = TRUE,
  error          = FALSE,
  fig.align      = "center",
  fig.width      = 8,
  fig.asp        = 0.618,
  message        = FALSE,
  warning        = FALSE
)
# Load tidyverse infrastructure packages
suppressPackageStartupMessages({
  library(tidyverse)
  library(here)
  library(ggpubr)
})

# Set paths
src_dir   <- here('code')
data_dir  <- here('2_Derived')
plots_dir <- here("3_Results")

# set seed
reseed <- 42
set.seed(seed = reseed)
```

## 1: READ VICTOR's FILE

```{r load}
Rep <- read.table(
    here(data_dir, "Homo_sapiens_triangles.txt"),
    sep = "\t",
    header = FALSE
)
names(Rep) <- c(
    "DelStart",
    "DelEnd",
    "DelLength",
    "MasterRepeat",
    "AlternRepeats1",
    "AlternRepeats2",
    "AlternRepeats3",
    "AlternRepeats4"
)
```

MasterRepeat is a repeat with arms located close to given deletion
breakpoints AlternativeRepeats1-4 are repeats where: first arm == arm1
of Master; first arm == arm2 of Master; second arm == arm1 of Master;
second arm == arm2 of Master; mb_del == realised deletions (exist in
MitoBreak); non_del - non realised deletions == don't exist in MitoBreak

## 2: filter out deletions within major arc:

-   OH: 110-441
-   OL: 5721-5781

```{r filt-maj-arc}
nrow(Rep)
Rep <- Rep[Rep$DelStart > 5781 &
    Rep$DelStart < 16569 &
    Rep$DelEnd > 5781 &
    Rep$DelEnd < 16569, ]
nrow(Rep)
```

## 3: concatenate all alternative deletions

```{r concat-alt-del}
Rep$AllAlternRepeats <- paste(
    Rep$AlternRepeats1,
    Rep$AlternRepeats2,
    Rep$AlternRepeats3,
    Rep$AlternRepeats4,
    sep = ","
)
Rep$AllAlternRepeats <- gsub("\\,\\[\\]\\,", ",", Rep$AllAlternRepeats)
Rep$AllAlternRepeats <- gsub("^\\[", "", Rep$AllAlternRepeats)
Rep$AllAlternRepeats <- gsub("\\]$", "", Rep$AllAlternRepeats)
```

## 4: prepare dataset for analysis of realized and nonrealized deletions line by line

```{r prepare}
Rep$CenterOfRealizedRepeats <- 0
Rep$CenterOfNonRealizedRepeats <- 0

Rep$LengthOfRealizedRepeats <- 0
Rep$LengthOfNonRealizedRepeats <- 0

Rep$StartOfRealizedRepeats <- 0
Rep$StartOfNonRealizedRepeats <- 0

Rep$EndOfRealizedRepeats <- 0
Rep$EndOfNonRealizedRepeats <- 0

for (i in (1:nrow(Rep))) {
    # i = 1
    # format of data to get a dataset of master and all alternative repeats for each deletion (for each line of the dataset)
    temp <- Rep[i, ]
    AltRep <- unlist(strsplit(temp$AllAlternRepeats, "\\]\\,\\["))
    AltRep <- AltRep[AltRep != ""]
    AltRep <- data.frame(AltRep)
    names(AltRep) <- c("WholeLine")
    if (nrow(AltRep) > 0) {
        AltRep$RepeatType <- "alternative"
        AltRep$WholeLine <- as.character(AltRep$WholeLine)

        MasterRepeat <- as.character(temp$MasterRepeat)
        MasterRepeat <- gsub("^\\[", "", MasterRepeat)
        MasterRepeat <- gsub("\\]$", "", MasterRepeat)

        MasterRepeat <- paste(MasterRepeat, "mb_del", sep = " ")
        MasterRepeat <- data.frame(MasterRepeat)
        names(MasterRepeat) <- c("WholeLine")
        MasterRepeat$RepeatType <- "master"

        AllRep <- rbind(MasterRepeat, AltRep)
        ReturnFifth <- function(x) {
            unlist(strsplit(x, " "))[5]
        }
        AllRep$RealisedRepeat <- apply(as.matrix(AllRep$WholeLine), 1, FUN = ReturnFifth)
        ReturnFirst <- function(x) {
            as.numeric(unlist(strsplit(x, " "))[1])
        }
        AllRep$RepStart <- apply(as.matrix(AllRep$WholeLine), 1, FUN = ReturnFirst)
        ReturnSecond <- function(x) {
            as.numeric(unlist(strsplit(x, " "))[2])
        }
        AllRep$RepEnd <- apply(as.matrix(AllRep$WholeLine), 1, FUN = ReturnSecond)
        AllRep <- AllRep[AllRep$RepStart > 5781 &
            AllRep$RepStart < 16569 &
            AllRep$RepEnd > 5781 &
            AllRep$RepEnd < 16569, ]
        if (nrow(AllRep) > 0) {
            AllRep$Center <- (AllRep$RepEnd - AllRep$RepStart) / 2 + AllRep$RepStart
            Rep$CenterOfRealizedRepeats[i] <- mean(AllRep[AllRep$RealisedRepeat == "mb_del", ]$Center)
            Rep$CenterOfNonRealizedRepeats[i] <- mean(AllRep[AllRep$RealisedRepeat == "non_del", ]$Center)

            AllRep$Length <- AllRep$RepEnd - AllRep$RepStart
            Rep$LengthOfRealizedRepeats[i] <- mean(AllRep[AllRep$RealisedRepeat == "mb_del", ]$Length)
            Rep$LengthOfNonRealizedRepeats[i] <- mean(AllRep[AllRep$RealisedRepeat == "non_del", ]$Length)

            Rep$StartOfNonRealizedRepeats[i] <- mean(AllRep[AllRep$RealisedRepeat == "non_del", ]$RepStart)
            Rep$StartOfRealizedRepeats[i] <- mean(AllRep[AllRep$RealisedRepeat == "mb_del", ]$RepStart)

            Rep$EndOfNonRealizedRepeats[i] <- mean(AllRep[AllRep$RealisedRepeat == "non_del", ]$RepEnd)
            Rep$EndOfRealizedRepeats[i] <- mean(AllRep[AllRep$RealisedRepeat == "mb_del", ]$RepEnd)

            if (i == 1) {
                FinalAllRep <- AllRep
            }
            if (i > 1) {
                FinalAllRep <- rbind(FinalAllRep, AllRep)
            }
        }
    }
}
nrow(Rep) # 703

Rep <- Rep[Rep$CenterOfRealizedRepeats > 0, ]
nrow(Rep) # 643
Rep <- Rep[Rep$CenterOfNonRealizedRepeats > 0, ] # it means there is at least one extra (third) arm to analyze
nrow(Rep) # 643
Rep <- Rep[!is.na(Rep$LengthOfRealizedRepeat), ] # 618
nrow(Rep)
```

## 5: Copare position of repeats with realized deletions vs rest repeats

```{r plt-boxplot-position, fig.width=4, fig.asp=1.618}
## center is a bit higher in realized repeats
wilcox.test(Rep$CenterOfRealizedRepeats,
    Rep$CenterOfNonRealizedRepeats,
    paired = TRUE
) # significant
t.test(Rep$CenterOfRealizedRepeats,
    Rep$CenterOfNonRealizedRepeats,
    paired = TRUE
) # significant
summary(Rep$CenterOfRealizedRepeats)
summary(Rep$CenterOfNonRealizedRepeats)
boxplot(Rep$CenterOfRealizedRepeats,
    Rep$CenterOfNonRealizedRepeats,
    notch = TRUE,
    names = c(
        "CenterOfRealizedRepeats",
        "CenterOfNonRealizedRepeats"
    )
)
```

```{r plt-violin-position}
RepReal <- Rep %>%
    select(
        Realized = CenterOfRealizedRepeats,
        NonRealized = CenterOfNonRealizedRepeats
    ) %>%
    gather(., Realized, NonRealized,
        key = "Deletion",
        value = "CenterOfRepeats"
    )

pltViolRepCenterRelease <- ggstatsplot::ggbetweenstats(
    data = RepReal,
    x = Deletion,
    y = CenterOfRepeats,
    notch = TRUE,
    # show notched box plot
    mean.ci = TRUE,
    # whether to display confidence interval for means
    k = 5,
    # number of decimal places for statistical results
    # outlier.tagging = TRUE, # whether outliers need to be tagged
    # outlier.label = ContactZone, # variable to be used for the outlier tag
    xlab = "Realisation of deletion",
    # label for the x-axis variable
    ylab = "Center of Repeats",
    # label for the y-axis variable
    title = "The effect of repeats' position on deletion realisation",
    # title text for the plot
    ggtheme = ggthemes::theme_fivethirtyeight(),
    # choosing a different theme
    ggstatsplot.layer = FALSE,
    # turn off `ggstatsplot` theme layer
    package = "wesanderson",
    # package from which color palette is to be taken
    palette = "Royal1",
    # choosing a different color palette
    messages = TRUE
)
# Note: Shapiro-Wilk Normality Test for Center of Repeats: p-value = 0.003
# Note: Bartlett's test for homogeneity of variances for factor Realisation of deletion: p-value = < 0.001
cowplot::save_plot(
    plot = pltViolRepCenterRelease,
    base_height = 8,
    base_asp = 1.618,
    file = normalizePath(
        file.path(plots_dir, "violin_rep_center_release.pdf")
    )
)

pltViolRepCenterRelease
```

## 6: Copare length of repeats with realized deletions vs rest repeats

```{r plt-boxplot-length, fig.width=4, fig.asp=1.618}
## deletion is longer in realized repeats
wilcox.test(
    Rep$LengthOfRealizedRepeats,
    Rep$LengthOfNonRealizedRepeats,
    paired = TRUE
) # significant
t.test(
    Rep$LengthOfRealizedRepeats,
    Rep$LengthOfNonRealizedRepeats,
    paired = TRUE
) # significant
summary(Rep$LengthOfRealizedRepeats) # 5643  mean
summary(Rep$LengthOfNonRealizedRepeats) # 3450 Median 3645 mean
summary(Rep$LengthOfRealizedRepeats - Rep$LengthOfNonRealizedRepeats) # 1998 mean
boxplot(Rep$LengthOfRealizedRepeats,
    Rep$LengthOfNonRealizedRepeats,
    notch = TRUE,
    names = c(
        "LengthOfRealizedRepeats",
        "LengthOfNonRealizedRepeats"
    )
)
```

```{r plt-violin-length}
RepReal <-
    Rep %>%
    select(
        Realized = LengthOfRealizedRepeats,
        NonRealized = LengthOfNonRealizedRepeats
    ) %>%
    gather(., Realized, NonRealized,
        key = "Deletion",
        value = "DistanceBetweenRepeats"
    ) %>%
    full_join(., RepReal)

pltViolRepFlankLengthRelease <-
    ggstatsplot::ggbetweenstats(
        data = RepReal,
        x = Deletion,
        y = DistanceBetweenRepeats,
        notch = TRUE,
        # show notched box plot
        mean.ci = TRUE,
        # whether to display confidence interval for means
        k = 5,
        # number of decimal places for statistical results
        # outlier.tagging = TRUE, # whether outliers need to be tagged
        # outlier.label = ContactZone, # variable to be used for the outlier tag
        xlab = "Realisation of deletion",
        # label for the x-axis variable
        ylab = "Distance between Repeats",
        # label for the y-axis variable
        title = "The effect of repeats' distance on deletion realisation",
        # title text for the plot
        ggtheme = ggthemes::theme_fivethirtyeight(),
        # choosing a different theme
        ggstatsplot.layer = FALSE,
        # turn off `ggstatsplot` theme layer
        package = "wesanderson",
        # package from which color palette is to be taken
        palette = "Royal1",
        # choosing a different color palette
        messages = TRUE
    )
# Note: Bartlett's test for homogeneity of variances for factor Realisation of deletion: p-value = < 0.001
cowplot::save_plot(
    plot = pltViolRepFlankLengthRelease,
    base_height = 8,
    base_asp = 1.618,
    file = normalizePath(
        file.path(plots_dir, "violin_rep_flanked_length_release.pdf")
    )
)

pltViolRepFlankLengthRelease
```

## 7: Start of repeats comparison for realized deletions vs rest repeats

```{r test-starts}
## start and end in realized versus non-realized repeats
wilcox.test(Rep$StartOfRealizedRepeats,
    Rep$StartOfNonRealizedRepeats,
    paired = TRUE
) # significant
t.test(Rep$StartOfRealizedRepeats,
    Rep$StartOfNonRealizedRepeats,
    paired = TRUE
) # significant
summary(Rep$StartOfRealizedRepeats) # mean 8678: 25%: 7401 - 75% 9764:
quantile(Rep$StartOfRealizedRepeats, 0.1) # 6465
quantile(Rep$StartOfRealizedRepeats, 0.9) # 10954

summary(Rep$StartOfNonRealizedRepeats) # mean 9394
summary(Rep$StartOfNonRealizedRepeats - Rep$StartOfRealizedRepeats) # mean 716 => non realised start later (~700 bp): 9394 - 8678

```

## 8: End of repeats comparison for realized deletions vs rest repeats

```{r test-ends}
wilcox.test(Rep$EndOfRealizedRepeats, Rep$EndOfNonRealizedRepeats, paired = TRUE) # significant
t.test(Rep$EndOfRealizedRepeats, Rep$EndOfNonRealizedRepeats, paired = TRUE) # significant
summary(Rep$EndOfRealizedRepeats) # 14321 mean
summary(Rep$EndOfNonRealizedRepeats) # 13039 mean
summary(Rep$EndOfNonRealizedRepeats - Rep$EndOfRealizedRepeats) # mean -1281=> non realised end earlier (~1300 bp): 13039 - 14321
quantile(Rep$EndOfRealizedRepeats, 0.1) # 13286
quantile(Rep$EndOfRealizedRepeats, 0.9) # 15863
```

## 9: Visualize it in terms of X(start) and Y(end)

```{r plt-simple-repeats-cols-by-realization}
plot(
    FinalAllRep[FinalAllRep$RealisedRepeat == "non_del", ]$RepStart,
    FinalAllRep[FinalAllRep$RealisedRepeat == "non_del", ]$RepEnd,
    pch = 16,
    col = "grey",
    xlim = c(5781, 16569),
    ylim = c(16569, 5781),
    xlab = "",
    ylab = ""
)
par(new = TRUE)
plot(
    FinalAllRep[FinalAllRep$RealisedRepeat == "mb_del", ]$RepStart,
    FinalAllRep[FinalAllRep$RealisedRepeat == "mb_del", ]$RepEnd,
    pch = 16,
    col = "red",
    xlim = c(5781, 16569),
    ylim = c(16569, 5781),
    xlab = "Start",
    ylab = "End"
)

```

```{r plt-margin-repeats-cols-by-realization}
# Visualise difference between density distributions of realized vs non-realized repeats.
sp <- ggplot(
  FinalAllRep,
  aes(
    RepStart,
    RepEnd
  )
) +
  aes(colour = RealisedRepeat) +
  geom_point() +
  scale_y_reverse() +
  theme_minimal(17) +
  xlab("5’") +
  ylab("3’") +
  scale_x_continuous(breaks = c(15000, 12000, 9000, 6000), position = "bottom") +
  scale_color_manual(
    values = c("red", "grey"),
    labels = c(
      "Realized repeats",
      "Non-realized repeats"
    )
  ) +
  theme(
    legend.position = "bottom",
    legend.title = element_blank(),
    legend.box = "horizontal"
  )

gg_realised <-
  ggExtra::ggMarginal(
    sp,
    type = "density",
    margins = "both",
    size = 3,
    groupColour = TRUE,
    groupFill = TRUE
  )

cowplot::save_plot(here(plots_dir, "real-vs-nonreal-repeats.svg"), gg_realised, base_asp = 0.868, base_height = 7)
```

```{r test-repeats-position-by-realisation}
wilcox.test(FinalAllRep[FinalAllRep$RealisedRepeat == "mb_del", ]$RepStart, FinalAllRep[FinalAllRep$RealisedRepeat == "non_del", ]$RepStart)
wilcox.test(FinalAllRep[FinalAllRep$RealisedRepeat == "mb_del", ]$RepEnd, FinalAllRep[FinalAllRep$RealisedRepeat == "non_del", ]$RepEnd)

FinalAllRep$RealisedRepeat <- as.factor(FinalAllRep$RealisedRepeat)
summary(FinalAllRep$RealisedRepeat)
```

## 10: Enrichment of repeats realised as a deletion in the contact zone

```{r}
StartA <- 6000 # 7000
StartB <- 9000 # 10000
EndA <- 13000
EndB <- 16000
R.In <- nrow(Rep[Rep$StartOfRealizedRepeats >= StartA &
    Rep$StartOfRealizedRepeats <= StartB &
    Rep$EndOfRealizedRepeats >= EndA &
    Rep$EndOfRealizedRepeats <= EndB, ])
R.Out <- nrow(Rep[(Rep$StartOfRealizedRepeats < StartA |
    Rep$StartOfRealizedRepeats > StartB) |
    (Rep$EndOfRealizedRepeats < EndA |
        Rep$EndOfRealizedRepeats > EndB), ])
R.In + R.Out # 324 + 294 = 618
NR.In <- nrow(Rep[Rep$StartOfNonRealizedRepeats >= StartA &
    Rep$StartOfNonRealizedRepeats <= StartB &
    Rep$EndOfNonRealizedRepeats >= EndA &
    Rep$EndOfNonRealizedRepeats <= EndB, ])
NR.Out <- nrow(Rep[(Rep$StartOfNonRealizedRepeats < StartA |
    Rep$StartOfNonRealizedRepeats > StartB) |
    (Rep$EndOfNonRealizedRepeats < EndA |
        Rep$EndOfNonRealizedRepeats > EndB), ])
NR.In + NR.Out # 79 + 598 = 618
X <- cbind(c(R.In, R.Out), c(NR.In, NR.Out))
fisher.test(X) # odds 7.5, p-value < 2.2e-16
```

Visualize it as mosaic plot

```{r contact-zone-deletion-realisation-of-repeat, fig.width=4}
X <- data.frame(X)
names(X) <- c("R", "NR")
row.names(X) <- c("IN", "OUT")
mosaicplot(X, main = "")
```

## Figure 2B:

```{r, fig.width=9, fig.asp=0.868}
gg_realised
```

Visualise difference between density distributions of realized vs
non-realized repeats.