Amplified_Prediction_script.Rmd

---
title: "Amplified Prediction Script"
output: html_notebook
---

#by Gandhar Mahadeshwar

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

```{r}
library(dplyr,warn.conflicts = FALSE)
library(foreign)
```

File intake.
```{r}
profile <- read.table("~/Desktop/23S Results/1mg2mn_ecoli_23s/1mg2mn_ecoli_23s_ecoli_23s_april23_2021_profile.txt", header = TRUE)
unt_mutation_counts <- read.table("~/Desktop/23S Results/1mg2mn_ecoli_23s/1mg2mn_ecoli_23s_Untreated_ecoli_23s_april23_2021_mutation_counts.txt", header = TRUE)
mod_mutation_counts <- read.table("~/Desktop/23S Results/1mg2mn_ecoli_23s/1mg2mn_ecoli_23s_Modified_ecoli_23s_april23_2021_mutation_counts.txt", header = TRUE)
```

We make two new columns for raw "correct" counts (unmutated).
```{r}
tmp_first <- profile %>% mutate(
  Modified_correct = Modified_effective_depth - Modified_mutations, 
  Untreated_correct = Untreated_effective_depth - Untreated_mutations
)
```

Select a few columns, rename others.
```{r}
tmp_filtered <- tmp_first %>% select(-ends_with("mapped_depth")) %>% select(-ends_with("read_depth")) %>% 
  select(-starts_with("Denatured")) %>% select(-ends_with("profile")) %>% select(-ends_with("err")) %>%
  dplyr::rename(
    Mn_mutations = Modified_mutations,
    Mn_effective_depth = Modified_effective_depth,
    Mn_rate = Modified_rate,
    Mn_correct = Modified_correct
    )
```


```{r}
mut_difference_table <- tmp_filtered %>% mutate(
  mut_difference = pmax((Mn_rate - Untreated_rate),0)) %>% select(c(1,11))

#write.csv(mut_difference_table,"../ForGandhar/SARS-naturalmods_data/sarscov2_mutdifftable.csv")

# plot(tmp_filtered$Mn_rate, tmp_filtered$Untreated_rate, xlim = c(0,0.08), ylim = c(0,0.08))
# 
# barplot(mut_difference_table$mut_difference, names.arg = mut_difference_table$Nucleotide)
```


Now, we look at the marginal mut rates file, calculate rates of certain marginal rates occurring. Also more column renaming.
```{r}
ratesfun <- function(x,y) {
  x / y
}

countedmut_Mn_filtered <- mod_mutation_counts %>% 
  mutate_all( ~ ratesfun(.,effective_depth)) %>%
  select(-ends_with("depth")) %>%
  dplyr::rename(
    'Mn_A-' = 'A.',
    'Mn_T-' = 'T.',
    'Mn_G-' = 'G.',
    'Mn_C-' = 'C.',
    'Mn_-A' = 'X.A',
    'Mn_-T' = 'X.T',
    'Mn_-G' = 'X.G',
    'Mn_-N' = 'X.N',
    'Mn_-C' = 'X.C',
    'Mn_AT' = 'AT',
    'Mn_AG' = 'AG',
    'Mn_AC' = 'AC',
    'Mn_TA' = 'TA',
    'Mn_TG' = 'TG',
    'Mn_TC' = 'TC',
    'Mn_GA' = 'GA',
    'Mn_GT' = 'GT',
    'Mn_GC' = 'GC',
    'Mn_CA' = 'CA',
    'Mn_CT' = 'CT',
    'Mn_CG' = 'CG',
    'Mn_multinuc_deletion' = 'multinuc_deletion',
    'Mn_multinuc_insertion' = 'multinuc_insertion',
    'Mn_multinuc_mismatch' = 'multinuc_mismatch',
    'Mn_complex_deletion' = 'complex_deletion',
    'Mn_complex_insertion' = 'complex_insertion'
  )

countedmut_Mn_filtered
```

Do the same for the untreated marginal mutation rates.
```{r}
countedmut_untreated_filtered <- unt_mutation_counts %>% 
  mutate_all( ~ ratesfun(.,effective_depth)) %>%
  select(-ends_with("depth")) %>%
  dplyr::rename(
    'A-' = 'A.',
    'T-' = 'T.',
    'G-' = 'G.',
    'C-' = 'C.',
    '-A' = 'X.A',
    '-T' = 'X.T',
    '-G' = 'X.G',
    '-N' = 'X.N',
    '-C' = 'X.C'
  )

countedmut_untreated_filtered
```

Bind three datasets together. Then, make the cumulative dataset (cumdf) robust by filtering.
```{r}
cumdf_pre <- cbind(tmp_filtered, countedmut_Mn_filtered, countedmut_untreated_filtered)
cumdf <- cumdf_pre %>% filter(Mn_effective_depth > 0) %>% filter(Untreated_effective_depth > 0)
```

```{r}
cumdf$`Mn_rate` <- as.numeric(cumdf$`Mn_rate`)

cumdf$`Mn_T-` <- cumdf$`Mn_T-` / cumdf$`Mn_rate`
cumdf$`Mn_A-` <- cumdf$`Mn_A-` / cumdf$`Mn_rate`
cumdf$`Mn_G-` <- cumdf$`Mn_G-` / cumdf$`Mn_rate`
cumdf$`Mn_C-` <- cumdf$`Mn_C-` / cumdf$`Mn_rate`
cumdf$`Mn_-A` <- cumdf$`Mn_-A` / cumdf$`Mn_rate`
cumdf$`Mn_-G` <- cumdf$`Mn_-G` / cumdf$`Mn_rate`
cumdf$`Mn_-C` <- cumdf$`Mn_-C` / cumdf$`Mn_rate`
cumdf$`Mn_-T` <- cumdf$`Mn_-T` / cumdf$`Mn_rate`
cumdf$`Mn_TA` <- cumdf$`Mn_TA` / cumdf$`Mn_rate`
cumdf$`Mn_TG` <- cumdf$`Mn_TG` / cumdf$`Mn_rate`
cumdf$`Mn_TC` <- cumdf$`Mn_TC` / cumdf$`Mn_rate`
cumdf$`Mn_GA` <- cumdf$`Mn_GA` / cumdf$`Mn_rate`
cumdf$`Mn_GT` <- cumdf$`Mn_GT` / cumdf$`Mn_rate`
cumdf$`Mn_GC` <- cumdf$`Mn_GC` / cumdf$`Mn_rate`
cumdf$`Mn_CA` <- cumdf$`Mn_CA` / cumdf$`Mn_rate`
cumdf$`Mn_CT` <- cumdf$`Mn_CT` / cumdf$`Mn_rate`
cumdf$`Mn_CG` <- cumdf$`Mn_CG` / cumdf$`Mn_rate`
cumdf$`Mn_AC` <- cumdf$`Mn_AC` / cumdf$`Mn_rate`
cumdf$`Mn_AG` <- cumdf$`Mn_AG` / cumdf$`Mn_rate`
cumdf$`Mn_AT` <- cumdf$`Mn_AT` / cumdf$`Mn_rate`

cumdf$`Untreated_rate` <- as.numeric(cumdf$`Untreated_rate`)

cumdf$`T-` <- cumdf$`T-` / cumdf$`Untreated_rate`
cumdf$`A-` <- cumdf$`A-` / cumdf$`Untreated_rate`
cumdf$`G-` <- cumdf$`G-` / cumdf$`Untreated_rate`
cumdf$`C-` <- cumdf$`C-` / cumdf$`Untreated_rate`
cumdf$`-A` <- cumdf$`-A` / cumdf$`Untreated_rate`
cumdf$`-G` <- cumdf$`-G` / cumdf$`Untreated_rate`
cumdf$`-C` <- cumdf$`-C` / cumdf$`Untreated_rate`
cumdf$`-T` <- cumdf$`-T` / cumdf$`Untreated_rate`
cumdf$`TA` <- cumdf$`TA` / cumdf$`Untreated_rate`
cumdf$`TG` <- cumdf$`TG` / cumdf$`Untreated_rate`
cumdf$`TC` <- cumdf$`TC` / cumdf$`Untreated_rate`
cumdf$`GA` <- cumdf$`GA` / cumdf$`Untreated_rate`
cumdf$`GT` <- cumdf$`GT` / cumdf$`Untreated_rate`
cumdf$`GC` <- cumdf$`GC` / cumdf$`Untreated_rate`
cumdf$`CA` <- cumdf$`CA` / cumdf$`Untreated_rate`
cumdf$`CT` <- cumdf$`CT` / cumdf$`Untreated_rate`
cumdf$`CG` <- cumdf$`CG` / cumdf$`Untreated_rate`
cumdf$`AC` <- cumdf$`AC` / cumdf$`Untreated_rate`
cumdf$`AG` <- cumdf$`AG` / cumdf$`Untreated_rate`
cumdf$`AT` <- cumdf$`AT` / cumdf$`Untreated_rate`

is.nan.data.frame <- function(x)
do.call(cbind, lapply(x, is.nan))

cumdf[is.nan(cumdf)] <- 0
cumdf

#cumdf <- cumdf[-c(which(cumdf$`TC` > 5)), ]
#cumdf_uracils_testing.mut <- cumdf_uracils_testing.mut[-c(which(cumdf_uracils_testing.mut$`TA` > 5)), ]
#cumdf_uracils_testing.mut <- cumdf_uracils_testing.mut[-c(which(cumdf_uracils_testing.mut$`TG` > 5)), ]
```

Separate into four datasets by nucleotide type.
```{r}
#Separate into four datasets by nucleotide type.
tmpa <- cumdf %>% filter(Sequence == 'A')
tmpu <- cumdf %>% filter(Sequence == 'U')
tmpc <- cumdf %>% filter(Sequence == 'C')
tmpg <- cumdf %>% filter(Sequence == 'G')

#adenines:total Mn2+ mutation rate above 0.5 were classified as a “m1A”
wc_cutoff_A <- which(tmpa$Mn_rate > 0.5)

#uridines:total Mn2+ mutation rate above 0.5 were classified as a “m3U or m1acp3Y”
wc_cutoff_U <- which(tmpu$Mn_rate > 0.5)

#guanosines: total Mn2+ mutation rate above 0.05 and total Mg2+ mutation rate above 0.0015 were classified as “m7G”
wc_cutoff_G <- which(tmpg$Mn_rate > 0.05 && tmpg$Untreated_rate > 0.0015)
```

Uracils final prediction set
```{r}
cumdf_uracils <- cumdf %>% filter(Sequence == 'U') %>% select(-('Mn_-T')) %>% select(-starts_with('Mn_A')) %>% 
  select(-starts_with('Mn_C')) %>% select(-starts_with('Mn_G')) %>% select(-('-T')) %>%
  select(-starts_with('A')) %>% select(-starts_with('C')) %>% select(-starts_with('G')) %>%
  select(-('Mn_-N')) %>% select(-('-N'))

cumdf_uracils_nuc <- cumdf_uracils$Nucleotide

cumdf_uracils_fin <- cumdf_uracils%>% select(-ends_with('depth')) %>% select(-ends_with('correct')) %>% 
  select(-ends_with('deletion')) %>% select(-ends_with('mismatch')) %>% select(-ends_with('insertion')) %>% 
  select(-ends_with('Nucleotide')) %>% select(-ends_with('Sequence')) %>% select(-ends_with('mutations')) %>% 
  select(-('T-')) %>% select(-('-A')) %>% select(-('-G')) %>% select(-('-C'))

modsq_uracils <- rep("?", length(cumdf_uracils_fin$Mn_rate))

sars2_u <- cbind(cumdf_uracils_fin, modsq_uracils) %>% rename('Modifications'='modsq_uracils')
sars2_u$Modifications <- as.factor(sars2_u$Modifications)
write.arff(sars2_u, file = "predict_uracils_amp.arff")

tx <- readLines("predict_uracils_amp.arff")
for (i in 1:length(tx)) {
  if (substr(tx[i],1,24)== "@attribute Modifications") {
    tx_line15 <- gsub(pattern = "'\\?'", replace = "none,pU,Um", x = tx[i])
    tx[i] <- tx_line15
    break
  }
}
tx3 <- gsub(pattern = "'", replace = "", x = tx)

writeLines(tx3, "predict_uracils_amp.arff")
```

Now, for cytosines
```{r}
cumdf_cytosines <- cumdf %>% filter(Sequence == 'C') %>% select(-('Mn_-C')) %>% select(-starts_with('Mn_A')) %>% 
  select(-starts_with('Mn_T')) %>% select(-starts_with('Mn_G')) %>% select(-('-C')) %>%
  select(-starts_with('A')) %>% select(-starts_with('T')) %>% select(-starts_with('G')) %>%
  select(-('Mn_-N')) %>% select(-('-N'))

cumdf_cytosines_nuc <- cumdf_cytosines$Nucleotide

cumdf_cytosines_fin <- cumdf_cytosines%>% select(-ends_with('depth')) %>% select(-ends_with('correct')) %>% 
  select(-ends_with('deletion')) %>% select(-ends_with('mismatch')) %>% select(-ends_with('insertion')) %>% 
  select(-ends_with('Nucleotide')) %>% select(-ends_with('Sequence')) %>% select(-ends_with('mutations')) %>% 
  select(-('C-')) %>% select(-('-A')) %>% select(-('-G')) %>% select(-('-T'))

modsq_cytosines <- rep("?", length(cumdf_cytosines_fin$Mn_rate))

sars2_c <- cbind(cumdf_cytosines_fin, modsq_cytosines) %>% rename('Modifications'='modsq_cytosines')
sars2_c$Modifications <- as.factor(sars2_c$Modifications)
write.arff(sars2_c, file = "predict_cytosines_amp.arff")

tx <- readLines("predict_cytosines_amp.arff")
for (i in 1:length(tx)) {
  if (substr(tx[i],1,24)== "@attribute Modifications") {
    tx_line15 <- gsub(pattern = "'\\?'", replace = "Cm,none", x = tx[i])
    tx[i] <- tx_line15
    break
  }
}
tx3 <- gsub(pattern = "'", replace = "", x = tx)

writeLines(tx3, "predict_cytosines_amp.arff")
```

Now, for guanines
```{r}
cumdf_guanine <- cumdf %>% filter(Sequence == 'G') %>% select(-('Mn_-G')) %>% select(-starts_with('Mn_A')) %>% 
  select(-starts_with('Mn_T')) %>% select(-starts_with('Mn_C')) %>% select(-('-G')) %>%
  select(-starts_with('A')) %>% select(-starts_with('T')) %>% select(-starts_with('C')) %>%
  select(-('Mn_-N')) %>% select(-('-N'))

cumdf_guanine_nuc <- cumdf_guanine$Nucleotide

cumdf_guanine_fin <- cumdf_guanine %>% select(-ends_with('depth')) %>% select(-ends_with('correct')) %>% 
  select(-ends_with('deletion')) %>% select(-ends_with('mismatch')) %>% select(-ends_with('insertion')) %>% 
  select(-ends_with('Nucleotide')) %>% select(-ends_with('Sequence')) %>% select(-ends_with('mutations')) %>% 
  select(-('G-')) %>% select(-('-A')) %>% select(-('-C')) %>% select(-('-T'))

modsq_guanine <- rep("?", length(cumdf_guanine_fin$Mn_rate))

sars2_g <- cbind(cumdf_guanine_fin, modsq_guanine) %>% rename('Modifications'='modsq_guanine')

sars2_g$Modifications <- as.factor(sars2_g$Modifications)
write.arff(sars2_g, file = "predict_guanine_amp.arff")

tx <- readLines("predict_guanine_amp.arff")
for (i in 1:length(tx)) {
  if (substr(tx[i],1,24)== "@attribute Modifications") {
    tx_line15 <- gsub(pattern = "'\\?'", replace = "Gm,none", x = tx[i])
    tx[i] <- tx_line15
    break
  }
}
tx3 <- gsub(pattern = "'", replace = "", x = tx)

writeLines(tx3, "predict_guanine_amp.arff")
```

Now, for adenines
```{r}
cumdf_adenine <- cumdf %>% filter(Sequence == 'A') %>% select(-('Mn_-A')) %>% select(-starts_with('Mn_G')) %>% 
  select(-starts_with('Mn_T')) %>% select(-starts_with('Mn_C')) %>% select(-('-A')) %>%
  select(-starts_with('G')) %>% select(-starts_with('T')) %>% select(-starts_with('C')) %>%
  select(-('Mn_-N')) %>% select(-('-N'))

cumdf_adenine_nuc <- cumdf_adenine$Nucleotide

cumdf_adenine_fin <- cumdf_adenine %>% select(-ends_with('depth')) %>% select(-ends_with('correct')) %>% 
  select(-ends_with('deletion')) %>% select(-ends_with('mismatch')) %>% select(-ends_with('insertion')) %>% 
  select(-ends_with('Nucleotide')) %>% select(-ends_with('Sequence')) %>% select(-ends_with('mutations')) %>% 
  select(-('A-')) %>% select(-('-G')) %>% select(-('-C')) %>% select(-('-T'))

modsq_adenine <- rep("?", length(cumdf_adenine_fin$Mn_rate))

sars2_a <- cbind(cumdf_adenine_fin, modsq_adenine) %>% rename('Modifications'='modsq_adenine')
sars2_a$Modifications <- as.factor(sars2_a$Modifications)

write.arff(sars2_a, file = "predict_adenines_amp.arff")

tx <- readLines("predict_adenines_amp.arff")
for (i in 1:length(tx)) {
  if (substr(tx[i],1,24)== "@attribute Modifications") {
    tx_line15 <- gsub(pattern = "'\\?'", replace = "Am,none", x = tx[i])
    tx[i] <- tx_line15
    break
  }
}
tx3 <- gsub(pattern = "'", replace = "", x = tx)

writeLines(tx3, "predict_adenines_amp.arff")
```