---
title: "Intership - tRNA studies"
author: "Makar Dorohuntsev"
format: html
editor: visual
embed-resources: true
---

## tRNA,TLSs and aaRS

<details>

<summary>Key theoretical aspects</summary>

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

### tRNA

In prokaryotes, tRNA genes are often organized in operons-clusters transcribed as a single RNA and later processed into individual tRNAs. Unlike eukaryotes, prokaryotic tRNA genes lack introns and require minimal processing (e.g., cleavage by RNase P, 3’ trimming, and CCA addition). Archaea show hybrid features, with some tRNA genes containing introns, resembling eukaryotic systems. tRNA abundance is tightly linked to cellular growth rates and environmental conditions.

In addition there is a wobble pairing present in mRNA-tRNA interactions, which allows for flexibility in codon-anti-codon pairing and usage of inosine at the third position of the anti-codon, enabling a single tRNA to recognize multiple codons. [Easy explanation](https://www.youtube.com/watch?v=ZUseoA_UW_4)

<details>

<summary>Wobble pairing influence on tRNA diversity</summary>

![](images/clipboard-3444307326.png)

[Wobble pairing analysis](https://www.frontiersin.org/journals/bioengineering-and-biotechnology/articles/10.3389/fbioe.2020.01031/full)

</details>

<details>

<summary>Wobble pairing evolution</summary>

![](images/clipboard-2889972363.png)

[Possible Co-evolution Theory of the Genetic Code](https://www.researchgate.net/publication/298409606_Coevolution_Theory_of_the_Genetic_Code_at_Age_Forty_Pathway_to_Translation_and_Synthetic_Life)

</details>

1.  tRNA is very varied between each other and have some key identity nucleotides and structures

<details>

<summary>tRNA types</summary>

<img src="tRNA types" style="width: 100%; max-width: 600px;"/>

[Page 4](https://cnrs.hal.science/hal-04023056v1)

</details>

2.  tRNA have different conformations and 3D structures depending on organism and it's type

<details>

<summary>tRNA structures</summary>

<img src="tRNA structures" style="width: 100%; max-width: 600px;"/>

[Page 15](https://cnrs.hal.science/hal-04023056v1)

</details>

3.  tRNA has complex life cycle in cell: It can be modified in nucleus and then used in mitochondrion or cytoplasm

<details>

<summary>tRNA cycle</summary>

<img src="tRNA cycle" style="width: 100%; max-width: 600px;"/>

[Page 4](https://www.tandfonline.com/doi/epdf/10.1080/15476286.2020.1809197?needAccess=true)

</details>

4.  tRNA has specific biogenesis and suggested to be modified after transcription

<details>

<summary>tRNA biogenesis</summary>

<img src="tRNA biogenesis" style="width: 100%; max-width: 600px;"/>

[Fig 1](https://www.cell.com/trends/biochemical-sciences/abstract/S0968-0004(20)30127-4?_returnURL=https%3A%2F%2Flinkinghub.elsevier.com%2Fretrieve%2Fpii%2FS0968000420301274%3Fshowall%3Dtrue)

</details>

5.  tRNA are recognized and verified during translation

<details>

<summary>tRNA in translation</summary>

<img src="tRNA recognition" style="width: 100%; max-width: 600px;"/>

[Page 8](https://www.tandfonline.com/doi/epdf/10.1080/15476286.2020.1809197?needAccess=true)

</details>

6.  Some tRNA nucleotides interact with distant ones during aaRS verification

<details>

<summary>Nucleotide interactions</summary>

<img src="Geom interactions" style="width: 100%; max-width: 600px;"/>

[Page 13](https://cnrs.hal.science/hal-04023056v1)

</details>

7.  Regarding alternative codons genesis, there could be arguments around statistical chance of mutations

<details>

<summary>Mutation risk vs alternative codons</summary>

Most of orange-red codons are susceptible to mutations which accrues between alternative codons, it could be explained by wobble mechanism which could lead to assignment of wrong tRNA anti-codon to mRNA codon. After inspection of standard codon table and alternative codon substitutions from NCBI database I can see clear tendency for third position to local change codon assignment to neighboring amino acid.

![](images/clipboard-3111457440.png)

DNA codon table

![](images/clipboard-1992191471.png)

[Alternative codons](https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi)

![](images/clipboard-2531825621.png)

[Mutation risk analysis](https://www.researchgate.net/publication/51200357_Codon_usage_in_vertebrates_is_associated_with_a_low_risk_of_acquiring_nonsense_mutations)

</details>

<details>

<summary>Useful databases</summary>

[Genetic Code in Taxonomy Tree](https://www.ncbi.nlm.nih.gov/Taxonomy/CommonTree/wwwcmt.cgi) - for group-specific analysis

[tRNA database](https://tpsic.igcz.poznan.pl/info/start/) - for secondary structure analysis and additional validation of search programs

</details>

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

### aaRS

Most information regarding structure and homology of aaRS can be found here: [AARS online](https://urzyme.github.io/)

1.  aaRS has a complex structure and is divided into two classes

<details>

<summary>aaRS classes and assigned amino acids</summary>

<img src="Amino classes" style="width: 100%; max-width: 600px;"/>

[Classes](https://en.wikipedia.org/wiki/Aminoacyl_tRNA_synthetase#cite_note-5)

</details>

2.  There are known anti-determinant elements in tRNA - aaRS interactions

<details>

<summary>Anti-determinant element examples</summary>

<img src="Antideterminant elements" style="width: 100%; max-width: 600px;"/>

[Page 13](https://cnrs.hal.science/hal-04023056v1)

</details>

3.  Mi-aaRS are coded in nuclear genome and are used in mitochondria, they have different structure and function compared to cytoplasmic aaRS

<details>

<summary>mi-aaRS pathway</summary>

![](images/clipboard-560633311.png)

[mi-aaRS in human disease context](https://www.researchgate.net/publication/318434152_Mitochondrial_Aminoacyl-tRNA_Synthetases_in_Human_Disease)

</details>

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

### Phage tRNA

1.  Phage encode their own tRNA genes to use host aaRS, more about this can be found in [Molecular biology of bacteriophage T4](https://archive.org/details/molecularbiology0000unse_i3c9/page/108/mode/1up?q=tRNA)

2.  Phage tRNA genes are often found in clusters. Their main role is to adapt host translation machinery to phage needs, like creating a codon usage bias to help with translation of A+T rich phage genome. Phage also cleave host tRNA genes to create same A+T bias:

    1.  tRNA-Leu (CUG codon)

    2.  [tRNA-Lys](https://www.researchgate.net/publication/221919247_RloC_A_Translation-Disabling_tRNase_Implicated_in_Phage_Exclusion_During_Recovery_from_DNA_Damage) with explained mechanism of cleavage

3.  Phages follow their host genetic code with noted exceptions [Phage code diversity](https://www.researchgate.net/publication/354165737_Stop_codon_recoding_is_widespread_in_diverse_phage_lineages_and_has_the_potential_to_regulate_translation_of_late_stage_and_lytic_genes) There are 3 consensus types of phage genetic code thriving in alternatively coded host:

-   Using same code as host (with own tRNA adopted to host code)

-   Adopted to host code with toleration of specific differences (without own tRNA, fully relying on host tRNA)

-   Using own code with specific mechanisms accounting for translation ("recoded phages" most common in extreme environments)

### TLSs

1.  TLS (tRNA-like structures) have similar functions and primary used by some viruses in hosts

<details>

<summary>TLS description</summary>

TLS have different structures and functions but mostly defined by type of amino acid attached through interaction with aaRS of host. They are proven to use anti-codon triplets and some identity elements of tRNA

[PubMed data](https://pubmed.ncbi.nlm.nih.gov/?term=tRNA-like)

</details>

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

### Additional information

1.  tRNA can be used by viruses to transport specific fragments

<details>

<summary>Viral Manipulation of Host tRNAs</summary>

<img src="tRNA virus1" style="width: 100%; max-width: 600px;"/>

[Page 6](https://www.cell.com/trends/biochemical-sciences/abstract/S0968-0004(20)30127-4?_returnURL=https%3A%2F%2Flinkinghub.elsevier.com%2Fretrieve%2Fpii%2FS0968000420301274%3Fshowall%3Dtrue)

</details>

<details>

<summary>tRNAs Are Used as Reverse Transcription Primers by Retroviruses</summary>

<img src="tRNA virus2" style="width: 100%; max-width: 600px;"/>

[Page 7](https://www.cell.com/trends/biochemical-sciences/abstract/S0968-0004(20)30127-4?_returnURL=https%3A%2F%2Flinkinghub.elsevier.com%2Fretrieve%2Fpii%2FS0968000420301274%3Fshowall%3Dtrue)

</details>

2.  There are known alternative genetic codes without dedicated STOP codons

[Source](https://pmc.ncbi.nlm.nih.gov/articles/PMC4967479/)

</details>

## Research notebook

<details>

<summary>Methodological backtrack</summary>

#### Research is divided in separate elements. Each element is described below:

<details>

<summary>Genomes used in analysis</summary>

Download links for 4 genomes (10.06.2025):

-   [Lancefieldella parvula DSM 20469](https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000024225.1/)
-   [Escherichia coli str. K-12 substr. MG1655](https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000005845.2/)
-   [Mesoplasma florum L1](https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000008305.1/)
-   [Kineococcus radiotolerans SRS30216 = ATCC BAA-149](https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000017305.1/)

</details>

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

<details>

<summary>tRNAscan-SE search</summary>

Use [tRNAscan](https://github.com/UCSC-LoweLab/tRNAscan-SE) to search for tRNA genes in downloaded genomes (13.06.2025):

``` bash
tRNAscan-SE 2.0.12 (Nov 2022)

tRNAscan-SE -B --brief Genomes/Ecoli.fna -o Outputs/Ecoli
tRNAscan-SE -B --brief Genomes/Kineococcus.fna -o Outputs/Kineococcus
tRNAscan-SE -B --brief Genomes/Lancefieldella.fna -o Outputs/Lancefieldella
tRNAscan-SE -B --brief Genomes/Mesoplasma.fna -o Outputs/Mesoplasma
```

After verification of completeness of tRNA genes all aminoacids appeared at least once in each genome except for SeC in Lancefieldella and Mesoplasma.

1.  Most amino acids comes in single coding DNA codon

    ![](images/clipboard-1357051920.png)

<!-- -->

2.  After verification of anticodon hits from tRNAscan-se resuts with codon usage database I can see numerous anticodons for one or more amino acid codons

![](images/clipboard-3706893255.png)

[Codon Usage Database](https://www.kazusa.or.jp/codon/)

</details>

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

<details>

<summary>ARAGORN search</summary>

As alternative approach I used [ARAGORN](https://academic.oup.com/nar/article/32/1/11/1194008) to search for tRNA genes in downloaded genomes (21.06.2025):

[Download using conda](https://anaconda.org/bioconda/aragorn)

``` bash
ARAGORN v1.2.41   Dean Laslett

aragorn -gcstd -w -t -o ARAGORN_output/Ecoli.txt Ecoli.fna
aragorn -gcstd -w -t -o ARAGORN_output/Kineococcus.txt Kineococcus.fna 
aragorn -gcstd -w -t -o ARAGORN_output/Lancefieldella.txt Lancefieldella.fna
aragorn -gcstd -w -t -o ARAGORN_output/Mesoplasma.txt Mesoplasma.fna
```

ARAGORN output differs from tRNAscan-SE by omitting some tRNA genes corresponding to amino acids present in tRNAscan output (like fMet or Ile2) or finding additional tRNA genes for amino acids absent in tRNAscan output. Additionally ARAGORN acquire more hits for some tRNA genes (like Arg in Kineococcus). tRNA sequences are longed in ARAGORN due to CAA ends which tRNAscan removes.

1.  Number of tRNA gene hits for each corresponding amino acid from ARAGORN output:

    ![](images/clipboard-1409456313.png)

2.  Analogically to tRNAscan made heatmaps using output from ARAGORN:

![](images/clipboard-579054460.png)

</details>

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

<details>

<summary>aaRS search</summary>

aaRS gene completeness check in genomes using MMSEQS2 (19.06.2025):

I used 24 sequences of Ecoli aaRS to search Ecoli and Kineococcus genomes:

[Proteins used for query](https://www.uniprot.org/uniprotkb?facets=model_organism%3A83333&fields=accession%2Cprotein_name%2Cgene_names%2Cdate_created%2Cgene_synonym%2Corganism_id&query=%28keyword%3AKW-0030%29+Ecoli&view=table)

<details>

<summary>Self verification</summary>

For sequence corelation check I used given query as search database and as search query:

Data for table was acquired using grep count on mmseq2 results:

``` bash
mmseqs createdb ../../Ecoli_aaRS.fasta aaRS_DB
mmseqs createdb ../../Ecoli_aaRS.fasta genome_DB
mmseqs search aaRS_DB genome_DB result_DB tmp --min-seq-id 0.3 -e 1e-50
mmseqs convertalis aaRS_DB genome_DB result_DB results.m8
cut -f1 Query:Prot_query.m8 | sed 's/^.*|//' | sed 's/_ECOLI.*//' | sort | uniq -c | sort -nr
```

``` bash
mmseqs search aaRS_DB genome_DB result_DB tmp --min-seq-id 0.3 -e 1e-9
mmseqs search aaRS_DB genome_DB result_DB tmp --min-seq-id 0.3 -e 1e-50
```

```{r,echo=FALSE, message=FALSE, warning=FALSE}
{raw_text1 <- "
3 Glutamyl_Q
3 Glutamine
3 Glutamate
2 Lysine_2
2 Lysine_1
2 Isoleucine
1 Valine
1 Tyrosine
1 Tryptophan
1 Threonine
1 Serine
1 Proline
1 Phenylalanine_B
1 Phenylalanine_A
1 Methionine
1 Leucine
1 Histidine
1 Glycine_B
1 Glycine_A
1 Cysteine
1 Aspartate
1 Asparagine
1 Arginine
1 Alanine
"
raw_text2 <- "
2 Lysine_2
2 Lysine_1
1 Valine
1 Tyrosine
1 Tryptophan
1 Threonine
1 Serine
1 Proline
1 Phenylalanine_B
1 Phenylalanine_A
1 Methionine
1 Leucine
1 Isoleucine
1 Histidine
1 Glycine_B
1 Glycine_A
1 Glutamyl_Q
1 Glutamine
1 Glutamate
1 Cysteine
1 Aspartate
1 Asparagine
1 Arginine
1 Alanine"
}
lines1 <- strsplit(trimws(raw_text1), "\n")[[1]]
parts1 <- strsplit(lines1, " ", fixed = TRUE)
lines2 <- strsplit(trimws(raw_text2), "\n")[[1]]
parts2 <- strsplit(lines2, " ", fixed = TRUE)

Ecoli_check <- data.frame(
  AminoAcid = sapply(parts1, function(x) x[2]),
  Count_e5 = as.numeric(sapply(parts1, function(x) x[1])),
  Count_e50 = as.numeric(sapply(parts2, function(x) x[1]))
)
knitr::kable(Ecoli_check, caption = "Self querry", format = "html", escape = FALSE)
```

Given the results, I can point some aaRS which has multiple hits, which means close sequence homology between them and other aaRS. To account for this I used e threshold of 1e-50 to reduce number of hits to only highly homologous sequences.

</details>

<details>

<summary>Ecoli output</summary>

As database for MMseqs2 I used CDS.faa acquired from prodigal prediction of [Ecoli.fna genome](https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000005845.2/)

Search command used:

``` bash
mmseqs search aaRS_DB genome_DB result_DB tmp --min-seq-id 0.3 -e 1e-50
```

```{r echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
library(dplyr)
library(kableExtra)
library(readr)
Ecoli_Ecoli <- read_delim("~/Praktyki_Bioinfa/aaRS_search/genomes_for_search/Ecoli:Ecoli/results.m8", 
    delim = "\t", escape_double = FALSE, 
    col_names = FALSE, trim_ws = TRUE)
Ecoli_Ecoli$X11 <- formatC(Ecoli_Ecoli$X11, digits = 6, format = "e")
Ecoli_Ecoli <- Ecoli_Ecoli %>%
  arrange(X1)
knitr::kable(Ecoli_Ecoli, caption = "Ecoli:Ecoli results", format = "html", escape = FALSE)
```

</details>

<details>

<summary>Kineococcus output</summary>

As database for MMseqs2 I used CDS.faa acquired from prodigal prediction of [Kineococcus.fna genome](https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000017305.1/)

Search command used:

``` bash
mmseqs search aaRS_DB genome_DB result_DB tmp --min-seq-id 0.3 -e 1e-50
```

```{r echo=FALSE, message=FALSE, warning=FALSE}
Ecoli_Kineococcus <- read_delim("~/Praktyki_Bioinfa/aaRS_search/genomes_for_search/Ecoli:Kineococcus/results.m8", 
    delim = "\t", escape_double = FALSE, 
    col_names = FALSE, trim_ws = TRUE)
Ecoli_Kineococcus$X11 <- formatC(Ecoli_Kineococcus$X11, digits = 6, format = "e")
Ecoli_Kineococcus <- Ecoli_Kineococcus %>%
  arrange(X1)
knitr::kable(Ecoli_Kineococcus, caption = "Ecoli:Kineococcus results", format = "html", escape = FALSE)
```

</details>

<details>

<summary>Lancefiedella output</summary>

As database for MMseqs2 I used CDS.faa acquired from prodigal prediction of [Lancefieldella.fna genome](https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000024225.1/)

Search command used:

``` bash
mmseqs search aaRS_DB genome_DB result_DB tmp --min-seq-id 0.3 -e 1e-50
```

```{r echo=FALSE, message=FALSE, warning=FALSE}
Ecoli_Lancefieldella <- read_delim("~/Praktyki_Bioinfa/aaRS_search/genomes_for_search/Ecoli:Lancefiedella/results.m8", 
    delim = "\t", escape_double = FALSE, 
    col_names = FALSE, trim_ws = TRUE)
Ecoli_Lancefieldella$X11 <- formatC(Ecoli_Lancefieldella$X11, digits = 6, format = "e")
Ecoli_Lancefieldella <- Ecoli_Lancefieldella %>%
  arrange(X1)
knitr::kable(Ecoli_Lancefieldella, caption = "Ecoli:Lancefiedella results", format = "html", escape = FALSE)
```

</details>

<details>

<summary>Mesoplasma output</summary>

As database for MMseqs2 I used CDS.faa acquired from prodigal prediction of [Mesoplasma.fna genome](https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000008305.1/)

Search command used:

``` bash
mmseqs search aaRS_DB genome_DB result_DB tmp --min-seq-id 0.3 -e 1e-50
```

```{r echo=FALSE, message=FALSE, warning=FALSE}
library(knitr)
Ecoli_Mesoplasma <- read_delim("~/Praktyki_Bioinfa/aaRS_search/genomes_for_search/Ecoli:Mesoplasma/results.m8", 
    delim = "\t", escape_double = FALSE, 
    col_names = FALSE, trim_ws = TRUE)
Ecoli_Mesoplasma$X11 <- formatC(Ecoli_Mesoplasma$X11, digits = 6, format = "e")
Ecoli_Mesoplasma <- Ecoli_Mesoplasma %>%
  arrange(X1)
knitr::kable(Ecoli_Mesoplasma, caption = "Ecoli:Mesoplasma results", format = "html", escape = FALSE)
```

</details>

In conclusion to aaRS search I can ascertain that most of the hits are not intersecting with highly homologous sequences from self-check, this

```{r echo=FALSE, message=FALSE, warning=FALSE}
library(dplyr)
library(tidyr)
library(knitr)

ecoli_aars <- data.frame(AminoAcid = unique(Ecoli_Ecoli$X1), Ecoli = "✓")
kineo_aars <- data.frame(AminoAcid = unique(Ecoli_Kineococcus$X1), Kineococcus = "✓")
lancefieldella_aars <- data.frame(AminoAcid = unique(Ecoli_Lancefieldella$X1), Lancefieldella = "✓")
mesoplasma_aars <- data.frame(AminoAcid = unique(Ecoli_Mesoplasma$X1), Mesoplasma = "✓")
combined <- full_join(ecoli_aars, kineo_aars, by = "AminoAcid") %>%
  mutate(
    Ecoli = ifelse(is.na(Ecoli), "✗", Ecoli),
    Kineococcus = ifelse(is.na(Kineococcus), "✗", Kineococcus)
  ) %>%
  arrange(Ecoli == "✗", Kineococcus == "✗")
combined <- full_join(combined, lancefieldella_aars, by = "AminoAcid") %>%
  mutate(
    Lancefieldella = ifelse(is.na(Lancefieldella), "✗", Lancefieldella)
  ) %>%
  arrange(Lancefieldella == "✗") %>%
  full_join(mesoplasma_aars, by = "AminoAcid") %>%
  mutate(
    Mesoplasma = ifelse(is.na(Mesoplasma), "✗", Mesoplasma)
  ) %>%
  arrange(Mesoplasma == "✗")

knitr::kable(
  combined,
  caption = "Presence of aaRS in E. coli and Kineococcus by amino acid",
  format = "html",
  escape = FALSE
)

```

Phylogenetic trees for validation of aaRS search hits:

<details>

<summary>Phylogenetic trees</summary>

1\. [Tree made using](https://itol.embl.de/tree/821652210732821751311211) ![](images/clipboard-4277154786.png)

2.  [Tree made using (16s based tree)](https://tygs.dsmz.de/)

    ![](images/clipboard-967686318.png){width="692" height="224"}

3.  [Tree made using (whole genome based tree)](https://tygs.dsmz.de/)

    ![](images/clipboard-2634125862.png){width="677" height="262"}

Additionally phylogenetic tree was created to validate and predict number of hits.

</details>

After self-check threshold for e-value was set to 1e-50. This potentially could lead to loss of some aaRS hits for more evolutionary distant organisms but in my opinion is necessary to reduce false positive hits.

Here is heatmap for E-val threshold check:

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height=10,fig.width=16}
library(ggplot2)
library(stringr)
files <- list.files(path = "~/Praktyki_Bioinfa/aaRS_search/genomes_for_search/Evale_scale ", pattern = "^E\\d+\\.txt$", full.names = TRUE)

# Przetwarzanie danych
df_all <- lapply(files, function(file) {
  evalue <- str_extract(basename(file), "E\\d+")
  data <- read.table(file, header = FALSE, col.names = c("count", "aa"))
  data$evalue <- evalue
  return(data)
}) %>% bind_rows()

df_heatmap <- df_all %>%
  group_by(evalue, aa) %>%
  summarise(count = sum(count), .groups = "drop") %>%
  mutate(evalue_numeric = as.numeric(str_remove(evalue, "E")))

ggplot(df_heatmap, aes(
  x = aa,
  y = reorder(evalue, evalue_numeric),
  fill = count
)) +
  geom_tile(color = "white", alpha = 0.6) +
  scale_fill_gradient(low = "blue", high = "red") +
  theme_minimal(base_size = 12) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1),
    panel.grid = element_blank()
  ) +
  labs(
    x = "aaRS domain",
    y = "E-value threshold",
    fill = "Hits"
  )
```

To compare similarity between aaRS I used Ecoli query and created two heatmaps first with raw evolutionary distances and second with scaled % similarity

Raw distances:

```{r echo=FALSE,  message=FALSE,output= FALSE ,warning=FALSE, fig.height= 16,fig.width= 16}
library(Biostrings)
library(stringr)
library(ComplexHeatmap)
library(DECIPHER)
library(circlize)
library(RColorBrewer)
# 1. Load sequences (replace with your actual file)
aars <- readAAStringSet("~/Praktyki_Bioinfa/aaRS_search/Ecoli_aaRS.fasta")

# 2. Extract simplified names (handles all your cases)
extract_aa_name <- function(full_name) {
  # Extract the part between | and | (e.g., "Alanine")
  aa_part <- str_extract(full_name, "\\|([^\\|]+)\\|") %>% 
    str_remove_all("\\|")
  
  # Handle special cases with suffixes
  case_when(
    str_detect(full_name, "Glutamyl_Q") ~ "Glutamyl_Q",
    str_detect(full_name, "Lysine_1") ~ "Lysine_1",
    str_detect(full_name, "Lysine_2") ~ "Lysine_2",
    str_detect(full_name, "Glycine_A") ~ "Glycine_A",
    str_detect(full_name, "Glycine_B") ~ "Glycine_B",
    str_detect(full_name, "Phenylalanine_A") ~ "Phenylalanine_A",
    str_detect(full_name, "Phenylalanine_B") ~ "Phenylalanine_B",
    aa_part == "Alanine" ~ "Alanine",
    aa_part == "Glutamate" ~ "Glutamate",
    # Add more mappings as needed
    TRUE ~ aa_part
  )
}

names(aars) <- sapply(names(aars), extract_aa_name)

# 3. Align sequences
aligned <- AlignSeqs(aars, processors = NULL)

# 4. Calculate similarity matrix (0-100%)
dist_mat <- DistanceMatrix(aligned, type = "matrix", correction = "none")
sim_mat <- round(100 * (1 - dist_mat/max(dist_mat)))
diag(sim_mat) <- 100  # Set self-comparisons to 100%

# 5. Create sensitive color gradient focused on biological range
color_gradient <- colorRampPalette(brewer.pal(9, "YlOrRd"))(100)
custom_breaks <- seq(floor(min(sim_mat)/10)*10, 100, length.out = 100)

ht_raw <- Heatmap(
  dist_mat,
  name = "Distance", 
  col = colorRamp2(
    breaks = c(min(dist_mat), median(dist_mat), max(dist_mat)),
    colors = c("blue", "white", "red")
  ),
  rect_gp = gpar(col = "white", lwd = 0.5),
  cluster_rows = TRUE,
  cluster_columns = TRUE,
  clustering_distance_rows = function(x) as.dist(x),
  clustering_distance_columns = function(x) as.dist(x),
  cell_fun = function(j, i, x, y, width, height, fill) {
    grid.text(
      sprintf("%.2f", dist_mat[i, j]),
      x, y, 
      gp = gpar(fontsize = 7,
                col = ifelse(dist_mat[i, j] > median(dist_mat), 
                             "white", "black"))
    )
  },
  row_names_gp = gpar(fontsize = 9),
  column_names_gp = gpar(fontsize = 9, rot = 45),
  heatmap_legend_param = list(
    at = round(c(min(dist_mat), median(dist_mat), max(dist_mat)), 2),
    title_position = "leftcenter-rot"
  )
)
```

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 16}
ht_raw
```

Scaled similarity:

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 16}
ht_scaled <- Heatmap(
  sim_mat,
  name = "Similarity (%)",
  col = colorRamp2(
    breaks = c(0, 50, 100),
    colors = c("red", "white", "blue")  # Reversed scale for similarity
  ),
  rect_gp = gpar(col = "white", lwd = 0.5),
  cluster_rows = TRUE,
  cluster_columns = TRUE,
  clustering_distance_rows = function(x) as.dist(x),
  clustering_distance_columns = function(x) as.dist(x),
  cell_fun = function(j, i, x, y, width, height, fill) {
    grid.text(
      sprintf("%.0f%%", sim_mat[i, j]),  # Whole percentages
      x, y, 
      gp = gpar(fontsize = 7,
                col = ifelse(sim_mat[i, j] < 50, 
                             "white", "black"))
    )
  },
  row_names_gp = gpar(fontsize = 9),
  column_names_gp = gpar(fontsize = 9, rot = 45),
  heatmap_legend_param = list(
    at = c(0, 50, 100),
    labels = c("0%", "50%", "100%"),
    title_position = "leftcenter-rot"
  )
)
ht_scaled
```

Selected points of similarity between aaRS (\>10% similarity):

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 12,fig.width= 12}
library(tibble)
library(tidyverse)
similarity_df <- as.data.frame(sim_mat) %>%
  rownames_to_column("Synthase1") %>%
  pivot_longer(cols = -Synthase1, names_to = "Synthase2", values_to = "Similarity") %>%
  filter(Synthase1 != Synthase2)

high_similarity_pairs <- similarity_df %>%
  filter(Similarity > 10) %>%
  mutate(
    Pair = map2_chr(Synthase1, Synthase2, ~paste(sort(c(.x, .y)), collapse = " - "))
  ) %>%
  distinct(Pair, .keep_all = TRUE) %>%
  arrange(desc(Similarity)) %>%
  dplyr::select(Synthase1, Synthase2, Similarity)

high_similarity_pairs <- high_similarity_pairs[-1,]

ggplot(high_similarity_pairs, aes(x = Synthase1, y = Synthase2, fill = Similarity)) +
  geom_tile(color = "white", linewidth = 0.5) +
  geom_text(aes(label = round(Similarity, 1)), 
            color = "white",
            size = 4,
            fontface = "bold") +
  scale_fill_gradientn(
    colors = c("#6D9EC1", "#E46726", "#B2182B"),
    name = "Similarity score",
    limits = c(min(high_similarity_pairs$Similarity, na.rm = TRUE), 
               max(high_similarity_pairs$Similarity, na.rm = TRUE))) +
      coord_fixed() +
      theme_minimal(base_size = 12) +
      labs(
        title = "High Similarity Pairs of Aminoacyl-tRNA Synthetases",
        subtitle = "Pairs with similarity score > 10",
        x = "Aminoacyl-tRNA Synthetase",
        y = "Aminoacyl-tRNA Synthetase"
      ) +
      theme(
        axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
        axis.text.y = element_text(hjust = 1),
        panel.grid = element_blank(),
        plot.title = element_text(face = "bold", hjust = 0.5)
      )

```

Coverage to identity plot:

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 6,fig.width= 6}
library(tidyverse)
library(ggrepel)
library(patchwork)

# Wczytanie danych z MMseqs2
mmseq_data <- read_tsv("~/Praktyki_Bioinfa/aaRS_search/genomes_for_search/Cov_Id/results.m8",
                       col_names = c("query", "target", "pident", "alnlen", "mismatches", "gaps", "qstart", "qend","tstart","tend",
                                     "evalue", "bitscore")) %>%
  mutate(
    # Poprawne obliczenie długości z alignmentu
    qlen = pmax(qend, alnlen),  # Długość query
    tlen = pmax(tend, alnlen),  # Długość targetu
    # Obliczenie pokrycia
    qcov = 100 * alnlen / qlen,
    tcov = 100 * alnlen / tlen,
    # Poprawna klasyfikacja
    status = case_when(
      pident >= 0.90 & qcov >= 90 & tcov >= 90 ~ "High",
      pident >= 0.70 & (qcov >= 70 | tcov >= 70) ~ "Medium", 
      TRUE ~ "Low"
    )
  ) 

# 2. Poprawiona wizualizacja
ggplot(mmseq_data, aes(x = pident, y = qcov, color = status)) +
  geom_point(alpha = 0.5) +
  geom_text_repel(
    aes(label = ifelse(query!=target,paste0(query, " to ",target),query)),
    size = 2.5,
    max.overlaps = 25
  ) +
  scale_color_manual(values = c("High" = "#1a9850", "Medium" = "#fdae61", "Low" = "#d73027")) +
  labs(
    title = "Cov to Identity in aaRS mmseq2 hits with eval > 1e-50",
    subtitle = "Status: High (ident ≥90%, cov ≥90%) | Medium (ident ≥70%, cov ≥70%)",
    x = "Identity (%)",
    y = "Coverage (%)"
  ) +
  theme_minimal()
```

</details>

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

<details>

<summary>Codetta search</summary>

I used [Codetta](https://github.com/kshulgina/codetta?tab=readme-ov-file) to predict genetic code of genomes and searched for any inconsistent results like STOP codons used in CDS sequences (22.06.2025). The hits were taken from output of alignment of profile HMM database to genomes nucl. sequences.

The analysis consists of three steps:

1.  Aligning the input nucleotide sequence to a database of profile Hidden Markov models (HMMs) of proteins (such as the Pfam database)
2.  Collating the resulting alignments into a single output file
3.  Inferring the genetic code from the alignment output file

For my analysis I used output of the second step, where I could search for individual alignments in HMM predictions for triplet in whole genome.

<details>

<summary>Ecoli search</summary>

``` bash
Ecoli predicted genetic code:
          FFLLSSSSYY??CC?WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
```

```{r echo=FALSE, message=FALSE, warning=FALSE}
codetta_hits_Ecoli <- read.csv("~/Praktyki_Bioinfa/data/codetta_hits/codetta_hits_Ecoli.txt", header = FALSE,sep = ",")
codetta_hits_Ecoli$V8 <- formatC(codetta_hits_Ecoli$V8, digits = 6, format = "e")
colnames(codetta_hits_Ecoli) <- c("C1","C2", "DNA_codon","C4","C5","C6", "Gene", "E-value","C9")
knitr::kable(
  codetta_hits_Ecoli,
  caption = "Ecoli hits for STOP codons",
  format = "html",
  escape = FALSE
)
```

</details>

<details>

<summary>Kineococcus search</summary>

``` bash
Kineococcus predicted genetic code:
          FF?LSSSSYY??CC?WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
```

```{r echo=FALSE, message=FALSE, warning=FALSE}
codetta_hits_Kineococcus <- read.csv("~/Praktyki_Bioinfa/data/codetta_hits/codetta_hits_Kineococcus.txt", header = FALSE,sep = ",")
codetta_hits_Kineococcus$V8 <- formatC(codetta_hits_Kineococcus$V8, digits = 6, format = "e")
colnames(codetta_hits_Kineococcus) <- c("C1","C2", "DNA_codon","C4","C5","C6", "Gene", "E-value","C9")
knitr::kable(
  codetta_hits_Kineococcus,
  caption = "Kineococcus hits for STOP codons",
  format = "html",
  escape = FALSE
)
```

</details>

<details>

<summary>Lancefieldella search</summary>

``` bash
Lancefieldella predicted genetic code:
          FFLLSSSSYY??CC?WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
```

``` bash
Lancefieldella hits for STOP codons:
10,15.5029800,TAA,15,5,29800,DUF6467,1.7e-13,33
```

```{r echo=FALSE, message=FALSE, warning=FALSE}
codetta_hits_Lancefieldella <- data.frame(C1 = 10,C2 = 15.5029800,DNA_codon = "TAA",C4 = 15,C5 = 5,C6 = 29800,Gene = "DUF6467",`E-value` = 1.7e-13,C9 = 33)

knitr::kable(
  codetta_hits_Lancefieldella,
  caption = "Mesoplasma hits for STOP codons",
  format = "html",
  escape = FALSE
)
```

</details>

<details>

<summary>Mesoplasma search</summary>

``` bash
Mesoplasma predicted genetic code with alternative codon:
          FFLLSSSSYY??CCWWLLLLPPPPHHQQRRR?IIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
Standard:              |*|
```

```{r echo=FALSE, message=FALSE, warning=FALSE}
codetta_hits_Mesoplasma <- read.csv("~/Praktyki_Bioinfa/data/codetta_hits/codetta_hits_Mesoplasma.txt", header = FALSE,sep = ",")
codetta_hits_Mesoplasma$V8 <- formatC(codetta_hits_Mesoplasma$V8, digits = 6, format = "e")
colnames(codetta_hits_Mesoplasma) <- c("C1","C2", "DNA_codon","C4","C5","C6", "Gene", "E-value","C9")
knitr::kable(
  codetta_hits_Mesoplasma,
  caption = "Mesoplasma hits for STOP codons",
  format = "html",
  escape = FALSE
)
```

In Mesoplasma were found 625 hits corresponding to substitution of TGA codon from STOP to Trp, which represent [alternative genetic code № 4](https://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi?chapter=cgencodes#SG4)

</details>

<details>

<summary>Codetta - tRNA genes cross-validation</summary>

In this part I want to cross-validate tRNA genes found with genetic code inferring from third step of Codetta output. I will use all tRNA genes from tRNAscan-SE and ARAGORN output and compare them with Codetta output to find nucleotide triplets which appeared in genome but have no corresponding tRNA gene. (28.06.2025)

Function:

1.  Takes tRNA gene sequences from tRNAscan-SE and ARAGORN as single list in + strand

2.  Translate it to RNA sequences for comparison with Codetta output

3.  Create list of unique triplets covered by tRNA genes additionally [accounting for wobble pairing](https://www.frontiersin.org/journals/bioengineering-and-biotechnology/articles/10.3389/fbioe.2020.01031/full)

4.  Map list on 4x4x4 matrix with codons and amino acids giving specific list of triplets which are not covered by tRNA genes in any way but were still counted in Codetta

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 16}
library(ggplot2)
library(gridExtra)
library(glue)

tRNA_codons_Ecoli <- c("GGC","TGC","ACG","CCG","CCT","TCT","GTT","GTC","GCA","CTG","TTG","TTC","CCC","GCC","TCC","GTG","GAT","CAT","CAA","CAG","GAG","TAA","TAG","TTT","GAA","CGG","GGG","TGG","TCA","CGA","GCT","GGA","TGA","CGT","GGT","TGT","CCA","GTA","NNN","GAC","TAC","TCG")
tRNA_codons_Kineococcus <- c("CGC","GGC","TGC","ACG","CCG","CCT","TCT","GTT","GTC","GCA","CTG","TTG","CTC","TTC","CCC","GCC","TCC","GTG","GAT","CAT","CAA","CAG","GAG","TAA","TAG","CTT","TTT","GAA","CGG","GGG","TGG","TCA","CGA","GCT","GGA","TGA","CGT","GGT","TGT","CCA","GTA","CAC","GAC","TAC","GCG","TCG")
tRNA_codons_Lancefieldella <- c("CGC","GGC","TGC","ACG","CCG","CCT","TCT","GTT","GTC","GCA","CTG","TTG","CTC","TTC","CCC","GCC","TCC","GTG","GAT","CAT","CAA","CAG","GAG","TAA","TAG","CTT","TTT","GAA","CGG","GGG","TGG","CGA","GCT","GGA","TGA","CGT","GGT","TGT","CCA","GTA","NNN","CAC","GAC","TAC")
tRNA_codons_Mesoplasma <- c("TGC","ACG","TCT","GTT","GTC","GCA","TTG","TTC","TCC","GTG","GAT","CAT","CAA","TAA","TAG","TTT","GAA","TGG","GCT","TGA","TCA","AGT","TGT","CCA","GTA","TAC")

#place for alternative code matrix
{
basic_code <- c(
  "UUU" = "F", "UUC" = "F", "UUA" = "L", "UUG" = "L",
  "UCU" = "S", "UCC" = "S", "UCA" = "S", "UCG" = "S",
  "UAU" = "Y", "UAC" = "Y", "UAA" = "*", "UAG" = "*",
  "UGU" = "C", "UGC" = "C", "UGA" = "*", "UGG" = "W",
  "CUU" = "L", "CUC" = "L", "CUA" = "L", "CUG" = "L",
  "CCU" = "P", "CCC" = "P", "CCA" = "P", "CCG" = "P",
  "CAU" = "H", "CAC" = "H", "CAA" = "Q", "CAG" = "Q",
  "CGU" = "R", "CGC" = "R", "CGA" = "R", "CGG" = "R",
  "AUU" = "I", "AUC" = "I", "AUA" = "I", "AUG" = "M",
  "ACU" = "T", "ACC" = "T", "ACA" = "T", "ACG" = "T",
  "AAU" = "N", "AAC" = "N", "AAA" = "K", "AAG" = "K",
  "AGU" = "S", "AGC" = "S", "AGA" = "R", "AGG" = "R",
  "GUU" = "V", "GUC" = "V", "GUA" = "V", "GUG" = "V",
  "GCU" = "A", "GCC" = "A", "GCA" = "A", "GCG" = "A",
  "GAU" = "D", "GAC" = "D", "GAA" = "E", "GAG" = "E",
  "GGU" = "G", "GGC" = "G", "GGA" = "G", "GGG" = "G"
)

code_4 <- c(
    "UUU" = "F", "UUC" = "F", "UUA" = "L", "UUG" = "L",
    "UCU" = "S", "UCC" = "S", "UCA" = "S", "UCG" = "S",
    "UAU" = "Y", "UAC" = "Y", "UAA" = "*", "UAG" = "*",
    "UGU" = "C", "UGC" = "C", "UGA" = "W", "UGG" = "W",
    "CUU" = "L", "CUC" = "L", "CUA" = "L", "CUG" = "L",
    "CCU" = "P", "CCC" = "P", "CCA" = "P", "CCG" = "P",
    "CAU" = "H", "CAC" = "H", "CAA" = "Q", "CAG" = "Q",
    "CGU" = "R", "CGC" = "R", "CGA" = "R", "CGG" = "R",
    "AUU" = "I", "AUC" = "I", "AUA" = "I", "AUG" = "M",
    "ACU" = "T", "ACC" = "T", "ACA" = "T", "ACG" = "T",
    "AAU" = "N", "AAC" = "N", "AAA" = "K", "AAG" = "K",
    "AGU" = "S", "AGC" = "S", "AGA" = "R", "AGG" = "R",
    "GUU" = "V", "GUC" = "V", "GUA" = "V", "GUG" = "V",
    "GCU" = "A", "GCC" = "A", "GCA" = "A", "GCG" = "A",
    "GAU" = "D", "GAC" = "D", "GAA" = "E", "GAG" = "E",
    "GGU" = "G", "GGC" = "G", "GGA" = "G", "GGG" = "G"
  )
}

map_tRNA_to_mRNA <- function(antykodony_DNA, species = "Organism", genetic_code) {
  tRNA_RNA <- toupper(chartr("T", "U", antykodony_DNA))
  tRNA_RNA <- tRNA_RNA[nchar(tRNA_RNA) == 3 & grepl("^[ACGU]+$", tRNA_RNA)]
  

  wobble_map <- list(
    "G" = c("U","C"),
    "C" = c("G","A"),
    "A" = c("U","C","A"),
    "U" = c("A","G","U","C"),
    "I" = c("A","U","C"),
    "*" = c("A","U","C", "G")
  )
  
  bazy <- c("U", "C", "A", "G")
  kodony_mRNA <- as.vector(outer(outer(bazy, bazy, paste0), bazy, paste0))
  
  
  rna_comp <- function(x) chartr("ACGU", "UGCA", x)
  
  
  rozpoznawalne <- unique(unlist(lapply(tRNA_RNA, function(anty) {
    a1 <- substr(anty, 1, 1)  # 5'
    a2 <- substr(anty, 2, 2)
    a3 <- substr(anty, 3, 3)  # 3'
    Filter(function(kodon) {
      k1 <- substr(kodon, 1, 1)  # 5'
      k2 <- substr(kodon, 2, 2)
      k3 <- substr(kodon, 3, 3)  # 3'
      k3 %in% wobble_map[[a1]] &&
        k2 == rna_comp(a2) &&
        k1 == rna_comp(a3)
    }, kodony_mRNA)
  })))
  
  
  if(genetic_code[51] == "*") {
    rozpoznawalne <- rozpoznawalne[!rozpoznawalne %in% c("UGA")]
  }
  if(genetic_code[35] == "*") {
    rozpoznawalne <- rozpoznawalne[!rozpoznawalne %in% c("UAA")]
  }
  if(genetic_code[36] == "*") {
    rozpoznawalne <- rozpoznawalne[!rozpoznawalne %in% c("UAG")]
  }
  
  
  df <- data.frame(
    codon = kodony_mRNA,
    b1 = substr(kodony_mRNA, 1, 1),
    b2 = substr(kodony_mRNA, 2, 2),
    b3 = substr(kodony_mRNA, 3, 3),
    aa = genetic_code[kodony_mRNA],
    status = ifelse(kodony_mRNA %in% rozpoznawalne, "Found", "Missing")
  ) %>%
    mutate(
      b1 = factor(b1, levels = rev(bazy)),
      b2 = factor(b2, levels = bazy),
      b3 = factor(b3, levels = bazy),
      label = paste0(codon, "\n", aa)
    )
  
  
  p <- ggplot(df, aes(x = b3, y = b1, fill = status)) +
    geom_tile(color = "black", width = 0.95, height = 0.95) +
    geom_text(aes(label = label), size = 3.2, fontface = "bold", lineheight = 0.8) +
    facet_wrap(~b2, nrow = 2, labeller = label_both) +
    scale_fill_manual(values = c("Found" = "lightblue", "Missing" = "pink")) +
    theme_minimal(base_size = 14) +
    theme(
      panel.grid = element_blank(),
      strip.text = element_text(size = 14),
      axis.title = element_text(size = 14),
      plot.title = element_text(hjust = 0.5)
    ) +
    labs(
      title = glue("Map of uncodable amino acids in {species}"),
      x = "Third nucleotide", y = "First nucleotide", fill = "Status"
    )
  
  return(list(
    wykres = p,
    rozpoznawalne_kodony = sort(rozpoznawalne),
    brakujace_kodony = sort(setdiff(kodony_mRNA, rozpoznawalne)),
    tRNA_RNA = tRNA_RNA
  ))
}

w1 <- map_tRNA_to_mRNA(tRNA_codons_Ecoli, species = "Ecoli",basic_code)
w2 <- map_tRNA_to_mRNA(tRNA_codons_Kineococcus, species = "Kineococcus",basic_code)
w3 <- map_tRNA_to_mRNA(tRNA_codons_Lancefieldella, species = "Lancefieldella",basic_code)
w4 <- map_tRNA_to_mRNA(tRNA_codons_Mesoplasma, species = "Mesoplasma",code_4)

grid.arrange(
  w1$wykres, w2$wykres, w3$wykres, w4$wykres,
  ncol = 2, top = "tRNA anticodon mapping to mRNA codons"
)
```

1.  STOP codons are missing with exception for UGA, which is present due to wobble pairing but was excluded in case of classic genetic code and genetic code without STOP codon substitution

2.  Complete codon sets in Ecoli, Kineococcus and Lancefieldella

3.  Missing tRNA gene accounting for CGG codon in Mesoplasma (It can be explained as statistically insignificant codon as codetta noted only 3 hit alignments of this codon)

In context of alternative genomes or rare amino acids, there is no significant difference in results because function maps tRNA anti-codons without specific amino acid. 4x4x4 matrix is part of overall biochemistry and does not change context of tRNA anti-codon mapping. In organisms living in extreme conditions internal mechanisms often minimize wobble pairing to stabilize genome and translation process. This could account for less number of covered codons than in classic wobble pairing matrix, but should not

</details>

</details>

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

<details>

<summary>ACPM pipeline for big-data analysis</summary>

I created bash pipeline using steps from previous analysis (7-8.07.2025):

1.  tRNAscan-SE and ARAGORN codons outputs collapsed into single file

2.  Codetta analysis with output of inference of genetic code

3.  R script with function for wobble accounting and plot generation

Output of this pipeline is a plot containing ACPM (Active Codon Profiled Matrix) of given genome and txt file with unused codons and type of genetic code if found in NCBI

Pipeline require specific directory structure and conda environment with tRNAscan, ARAGORN and Codetta

``` bash
conda create -n trna_pipeline_env python=3.10

conda activate trna_pipeline_env

conda config --add channels defaults
conda config --add channels bioconda
conda config --add channels conda-forge

conda install trnascan-se aragorn
```

In addition tRNAscan, ARAGORN and [Codetta](https://github.com/kshulgina/codetta?tab=readme-ov-file) can require specific packages which should be downloaded during setup.

Codetta folder should be downloaded in directory "codetta" and genome .fna's or .fasta'a should be downloaded in directory "genomes".

Overall look of full pipeline directory should look like this:

Pipeline/

-   pipeline.sh

-   codon_plot.R

-   codetta/

-   genomes/

-   temp/

-   results/

    -   genome_name_codon_summary.txt

    -   genome_name_codon_histogram.png

    -   ...

<details>

<summary>bash pipeline</summary>

```{bash eval=FALSE}
#!/bin/bash

# Base path (script location)
BASE_DIR=$(dirname "$0")
GENOMES_DIR="$BASE_DIR/genomes"
TEMP_DIR="$BASE_DIR/temp"
RESULTS_DIR="$BASE_DIR/results"
CODETTA_DIR="$BASE_DIR/codetta"

# Create directories if they don't exist
mkdir -p "$TEMP_DIR"
mkdir -p "$RESULTS_DIR"

# Iterate over all .fna files
for genome_path in "$GENOMES_DIR"/*.fna; do
    # Filename without path and extension
    genome_file=$(basename "$genome_path")
    genome_name="${genome_file%.fna}"

    echo "Processing: $genome_name"

    # Temporary files
    trna_out="$TEMP_DIR/${genome_name}_tRNAscan.out"
    aragorn_out="$TEMP_DIR/${genome_name}_aragorn.out"
    codon_summary="$TEMP_DIR/${genome_name}_codons_summary.txt"
    codetta_file="$TEMP_DIR/${genome_name}_codetta_results.txt"

    # Run tRNAscan-SE
    tRNAscan-SE -B --brief -o "$trna_out" "$genome_path"

    # Run ARAGORN
    aragorn -gcstd -w -t -o "$aragorn_out" "$genome_path"

    # Run Codetta
    python3.10 "$CODETTA_DIR/codetta.py" -e 1e-50 --results_summary "$codetta_file" "$genome_path"

    # Extract codons from tRNAscan-SE
    if [[ -f "$trna_out" ]]; then
        awk '$6 ~ /^[ATGC]{3}$/ { print $6 }' "$trna_out" > "$TEMP_DIR/${genome_name}_codons_tRNAscan.txt"
    else
        echo "Missing file $trna_out"
    fi

    # Extract codons from ARAGORN
    if [[ -f "$aragorn_out" ]]; then
        grep -v "^#" "$aragorn_out" | \
        awk '{ gsub(/\(|\)/, "", $5); codon=toupper($5); if (codon ~ /^[ATGC]{3}$/) print codon }' \
        > "$TEMP_DIR/${genome_name}_codons_aragorn.txt"
    else
        echo "Missing file $aragorn_out"
    fi

    # Combine and count unique codons
    cat "$TEMP_DIR/${genome_name}_codons_tRNAscan.txt" "$TEMP_DIR/${genome_name}_codons_aragorn.txt" \
        | sort | uniq -c | awk '{print $2","$1}' > "$codon_summary"

    # Run R script for visualization and analysis
    histogram_file="$RESULTS_DIR/${genome_name}_codon_histogram.png"
    unused_codons="$RESULTS_DIR/${genome_name}_unused_codons.txt"

    if [[ -f "$codon_summary" && -f "$codetta_file" ]]; then
        Rscript codon_plot.R "$codon_summary" "$unused_codons" "$histogram_file" "$genome_name" "$codetta_file"
    else
        echo "Missing file $codon_summary or $codetta_file for visualization."
    fi

    # Remove all temporary files for this genome
    rm -f "$TEMP_DIR/${genome_name}"_*
    find "$GENOMES_DIR" -type f -name "${genome_name}*" ! -name "${genome_name}.fna" ! -name "${genome_name}.fasta" -delete


    echo "Saved results for $genome_name in results/ directory"
done
```

</details>

<details>

<summary>R script</summary>

```{r eval=FALSE}
library(ggplot2)
library(glue)
# Define standard genetic code
basic_code <- c(
  "UUU" = "F", "UUC" = "F", "UUA" = "L", "UUG" = "L",
  "UCU" = "S", "UCC" = "S", "UCA" = "S", "UCG" = "S",
  "UAU" = "Y", "UAC" = "Y", "UAA" = "*", "UAG" = "*",
  "UGU" = "C", "UGC" = "C", "UGA" = "*", "UGG" = "W",
  "CUU" = "L", "CUC" = "L", "CUA" = "L", "CUG" = "L",
  "CCU" = "P", "CCC" = "P", "CCA" = "P", "CCG" = "P",
  "CAU" = "H", "CAC" = "H", "CAA" = "Q", "CAG" = "Q",
  "CGU" = "R", "CGC" = "R", "CGA" = "R", "CGG" = "R",
  "AUU" = "I", "AUC" = "I", "AUA" = "I", "AUG" = "M",
  "ACU" = "T", "ACC" = "T", "ACA" = "T", "ACG" = "T",
  "AAU" = "N", "AAC" = "N", "AAA" = "K", "AAG" = "K",
  "AGU" = "S", "AGC" = "S", "AGA" = "R", "AGG" = "R",
  "GUU" = "V", "GUC" = "V", "GUA" = "V", "GUG" = "V",
  "GCU" = "A", "GCC" = "A", "GCA" = "A", "GCG" = "A",
  "GAU" = "D", "GAC" = "D", "GAA" = "E", "GAG" = "E",
  "GGU" = "G", "GGC" = "G", "GGA" = "G", "GGG" = "G"
)

# Function to create genome-specific genetic code from Codetta CSV output
create_genome_specific_code <- function(codetta_file, base_code) {
  # Read Codetta output as CSV
  codetta_data <- read.csv(codetta_file, header = TRUE, stringsAsFactors = FALSE)
  
  # Extract inferred_gencode (assuming one row per file)
  genetic_code_str <- codetta_data$inferred_gencode[1]
  
  # Define codon order as per Codetta (UUU, UUC, ..., GGG)
  codons <- c("UUU", "UUC", "UUA", "UUG", "UCU", "UCC", "UCA", "UCG",
              "UAU", "UAC", "UAA", "UAG", "UGU", "UGC", "UGA", "UGG",
              "CUU", "CUC", "CUA", "CUG", "CCU", "CCC", "CCA", "CCG",
              "CAU", "CAC", "CAA", "CAG", "CGU", "CGC", "CGA", "CGG",
              "AUU", "AUC", "AUA", "AUG", "ACU", "ACC", "ACA", "ACG",
              "AAU", "AAC", "AAA", "AAG", "AGU", "AGC", "AGA", "AGG",
              "GUU", "GUC", "GUA", "GUG", "GCU", "GCC", "GCA", "GCG",
              "GAU", "GAC", "GAA", "GAG", "GGU", "GGC", "GGA", "GGG")
  
  # Split genetic code string into individual amino acids
  inferred_aa <- strsplit(genetic_code_str, "")[[1]]
  
  # Create new genetic code matrix
  genome_code <- base_code
  for (i in 1:length(codons)) {
    inferred <- inferred_aa[i]
    standard <- base_code[codons[i]]
    # Update only if inferred is not "?" and different from standard
    if (inferred != "?" && inferred != standard) {
      genome_code[codons[i]] <- inferred
    }
  }
  
  return(genome_code)
}

# Function to map tRNA anticodons to mRNA codons
map_tRNA_to_mRNA <- function(antykodony_DNA, species = "Organism", genetic_code) {
  tRNA_RNA <- toupper(chartr("T", "U", antykodony_DNA))
  tRNA_RNA <- tRNA_RNA[nchar(tRNA_RNA) == 3 & grepl("^[ACGU]+$", tRNA_RNA)]
  
  wobble_map <- list(
    "G" = c("U","C"),
    "C" = c("G","A"),
    "A" = c("U","C","A"),
    "U" = c("A","G","U","C"),
    "I" = c("A","U","C"),
    "*" = c("A","U","C", "G")
  )
  
  bazy <- c("U", "C", "A", "G")
  kodony_mRNA <- as.vector(outer(outer(bazy, bazy, paste0), bazy, paste0))
  
  rna_comp <- function(x) chartr("ACGU", "UGCA", x)
  
  rozpoznawalne <- unique(unlist(lapply(tRNA_RNA, function(anty) {
    a1 <- substr(anty, 1, 1)  # 5'
    a2 <- substr(anty, 2, 2)
    a3 <- substr(anty, 3, 3)  # 3'
    Filter(function(kodon) {
      k1 <- substr(kodon, 1, 1)  # 5'
      k2 <- substr(kodon, 2, 2)
      k3 <- substr(kodon, 3, 3)  # 3'
      k3 %in% wobble_map[[a1]] &&
        k2 == rna_comp(a2) &&
        k1 == rna_comp(a3)
    }, kodony_mRNA)
  })))
  
  if(genetic_code[51] == "*") {
    rozpoznawalne <- rozpoznawalne[!rozpoznawalne %in% c("UGA")]
  }
  if(genetic_code[35] == "*") {
    rozpoznawalne <- rozpoznawalne[!rozpoznawalne %in% c("UAA")]
  }
  if(genetic_code[36] == "*") {
    rozpoznawalne <- rozpoznawalne[!rozpoznawalne %in% c("UAG")]
  }
  
  df <- data.frame(
    codon = kodony_mRNA,
    b1 = substr(kodony_mRNA, 1, 1),
    b2 = substr(kodony_mRNA, 2, 2),
    b3 = substr(kodony_mRNA, 3, 3),
    aa = genetic_code[kodony_mRNA],
    status = ifelse(kodony_mRNA %in% rozpoznawalne, "Found", "Missing")
  )
  df$b1 <- factor(df$b1, levels = rev(bazy))
  df$b2 <- factor(df$b2, levels = bazy)
  df$b3 <- factor(df$b3, levels = bazy)
  df$label <- paste0(df$codon, "\n", df$aa)
  
  p <- ggplot(df, aes(x = b3, y = b1, fill = status)) +
    geom_tile(color = "black", width = 0.95, height = 0.95) +
    geom_text(aes(label = label), size = 3.2, fontface = "bold", lineheight = 0.8) +
    facet_wrap(~b2, nrow = 2, labeller = label_both) +
    scale_fill_manual(values = c("Found" = "lightblue", "Missing" = "pink")) +
    theme_minimal(base_size = 14) +
    theme(
      panel.grid = element_blank(),
      strip.text = element_text(size = 14),
      axis.title = element_text(size = 14),
      plot.title = element_text(hjust = 0.5)
    ) +
    labs(
      title = glue("Map of uncodable amino acids in {species}"),
      x = "Second nucleotide", y = "First nucleotide", fill = "Status"
    )
  
  return(list(
    wykres = p,
    rozpoznawalne_kodony = sort(rozpoznawalne),
    brakujace_kodony = sort(setdiff(kodony_mRNA, rozpoznawalne)),
    tRNA_RNA = tRNA_RNA
  ))
}

# Read command line arguments
args <- commandArgs(trailingOnly = TRUE)
codon_summary_file <- args[1]
unused_codons_file <- args[2]
histogram_file <- args[3]
species_name <- args[4]
codetta_file <- args[5]

# Read codon summary data
codon_data <- read.csv(codon_summary_file, header = FALSE, col.names = c("Codon", "Count"))

# Map DNA codons to RNA
codon_data$Codon <- toupper(chartr("T", "U", codon_data$Codon))

# Create genome-specific genetic code from Codetta output
genome_specific_code <- create_genome_specific_code(codetta_file, basic_code)

# Use genome-specific genetic code in analysis
result <- map_tRNA_to_mRNA(codon_data$Codon, species = species_name, genetic_code = genome_specific_code)

# Save plot
ggsave(histogram_file, plot = result$wykres, width = 10, height = 8, dpi = 300)
# Save unused codons
writeLines(c(result$brakujace_kodony,"","Genetic code: ",paste(genome_specific_code,collapse = "")),  unused_codons_file)
```

</details>

In addition here is short bash code to cound specific genetic codes for big data analysis (inside results folder):

```{bash eval=FALSE}
grep -h -A1 "Genetic code:" *.txt | grep -v "Genetic code:" | grep -v "^--$" | sort | uniq -c | sort -nr
```

</details>

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

<details>

<summary>Phage in tRNA research</summary>

Inspired by [Molecular biology of bacteriophage T4](https://archive.org/details/molecularbiology0000unse_i3c9/page/108/mode/1up?q=tRNA) I want to explore phage tRNA genes and possible logic hidden behind genomic data.

<details>

<summary>Preliminary ACPM analysis</summary>

[Escherichia phage T4](https://www.ncbi.nlm.nih.gov/datasets/taxonomy/2681598/) contains 8 tRNA genes, main theoretical purpose of which is to adopt to host and be more efficient in translation. (14-15.07.2025)

Pipeline analysis showed standard genetic code and not translation capable ACPM:

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 16}
library(tidyverse)
df <- read_csv("~/Praktyki_Bioinfa/Pipeline/Phage_results/tidy_codon_matrix.csv")

bazy <- c("U", "C", "A", "G")
kodony_mRNA <- as.vector(outer(outer(bazy, bazy, paste0), bazy, paste0))

phage_codons <- df %>%
  filter(Genome == "all_codons")
rozpoznawalne <- phage_codons$Codon

{
  base_code <- c(
    "UUU" = "F", "UUC" = "F", "UUA" = "L", "UUG" = "L",
    "UCU" = "S", "UCC" = "S", "UCA" = "S", "UCG" = "S",
    "UAU" = "Y", "UAC" = "Y", "UAA" = "*", "UAG" = "*",
    "UGU" = "C", "UGC" = "C", "UGA" = "*", "UGG" = "W",
    "CUU" = "L", "CUC" = "L", "CUA" = "L", "CUG" = "L",
    "CCU" = "P", "CCC" = "P", "CCA" = "P", "CCG" = "P",
    "CAU" = "H", "CAC" = "H", "CAA" = "Q", "CAG" = "Q",
    "CGU" = "R", "CGC" = "R", "CGA" = "R", "CGG" = "R",
    "AUU" = "I", "AUC" = "I", "AUA" = "I", "AUG" = "M",
    "ACU" = "T", "ACC" = "T", "ACA" = "T", "ACG" = "T",
    "AAU" = "N", "AAC" = "N", "AAA" = "K", "AAG" = "K",
    "AGU" = "S", "AGC" = "S", "AGA" = "R", "AGG" = "R",
    "GUU" = "V", "GUC" = "V", "GUA" = "V", "GUG" = "V",
    "GCU" = "A", "GCC" = "A", "GCA" = "A", "GCG" = "A",
    "GAU" = "D", "GAC" = "D", "GAA" = "E", "GAG" = "E",
    "GGU" = "G", "GGC" = "G", "GGA" = "G", "GGG" = "G"
  )
  
  genetic_code_str <- "FFLLSSSSYY**CC*WLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG"
  # Define codon order as per Codetta (UUU, UUC, ..., GGG)
  codons <- c("UUU", "UUC", "UUA", "UUG", "UCU", "UCC", "UCA", "UCG",
              "UAU", "UAC", "UAA", "UAG", "UGU", "UGC", "UGA", "UGG",
              "CUU", "CUC", "CUA", "CUG", "CCU", "CCC", "CCA", "CCG",
              "CAU", "CAC", "CAA", "CAG", "CGU", "CGC", "CGA", "CGG",
              "AUU", "AUC", "AUA", "AUG", "ACU", "ACC", "ACA", "ACG",
              "AAU", "AAC", "AAA", "AAG", "AGU", "AGC", "AGA", "AGG",
              "GUU", "GUC", "GUA", "GUG", "GCU", "GCC", "GCA", "GCG",
              "GAU", "GAC", "GAA", "GAG", "GGU", "GGC", "GGA", "GGG")
  
  # Split genetic code string into individual amino acids
  inferred_aa <- strsplit(genetic_code_str, "")[[1]]
  
  # Create new genetic code matrix
  genome_code <- base_code
  for (i in 1:length(codons)) {
    inferred <- inferred_aa[i]
    standard <- base_code[codons[i]]
    # Update only if inferred is not "?" and different from standard
    if (inferred != "?" && inferred != standard) {
      genome_code[codons[i]] <- inferred
    }
  }
  
}
genetic_code <- genome_code
species <- "Tequatrovirus taxon"
df <- data.frame(
  codon = kodony_mRNA,
  b1 = substr(kodony_mRNA, 1, 1),
  b2 = substr(kodony_mRNA, 2, 2),
  b3 = substr(kodony_mRNA, 3, 3),
  aa = genetic_code[kodony_mRNA],
  status = ifelse(kodony_mRNA %in% rozpoznawalne, "Found", "Missing")
)
df$b1 <- factor(df$b1, levels = rev(bazy))
df$b2 <- factor(df$b2, levels = bazy)
df$b3 <- factor(df$b3, levels = bazy)
df$label <- paste0(df$codon, "\n", df$aa)

ggplot(df, aes(x = b3, y = b1, fill = status)) +
  geom_tile(color = "black", width = 0.95, height = 0.95) +
  geom_text(aes(label = label), size = 6, fontface = "bold", lineheight = 0.8) +
  facet_wrap(~b2, nrow = 2, labeller = label_both) +
  scale_fill_manual(values = c("Found" = "lightblue", "Missing" = "pink")) +
  theme_minimal(base_size = 14) +
  theme(
    panel.grid = element_blank(),
    strip.text = element_text(size = 14),
    axis.title = element_text(size = 14),
    plot.title = element_text(hjust = 0.5)
  ) +
  labs(
    title = glue("Map of uncodable amino acids in {species}"),
    x = "Third nucleotide", y = "First nucleotide", fill = "Status"
  )
```

After preliminary singe genome analysis I downloaded [Tequatrovirus](https://www.ncbi.nlm.nih.gov/datasets/genome/?taxon=10663) taxon (347 phage genomes) from NCBI and created pipeline to perform statistical check of codon frequency and ACPM for all genomes in this taxon. This pipeline resulted in long format table with information about ACPM for every genome (every possible codon covered by present tRNA genes) (16.07.2025-18.07.2025)

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 16}
phage_codons <- read_csv("~/Praktyki_Bioinfa/Pipeline/Phage_results/tidy_codon_matrix.csv")

knitr::kable(
  phage_codons %>%
    group_by(Genome) %>%
    summarize(Codons_found = sum(Presence == 1)) %>%
    filter(Codons_found == 0),
  caption = "Genomes with no covered codons",
  format = "html",
  escape = FALSE
)

knitr::kable(
  phage_codons %>%
    group_by(Genome) %>%
    summarize(Codons_found = sum(Presence == 1)) %>%
    filter(Codons_found == 36 & Genome != "all_codons"),
  caption = "Genomes with max covered codons",
  format = "html",
  escape = FALSE
)

df <- data.frame(
  codon = kodony_mRNA,
  b1 = substr(kodony_mRNA, 1, 1),
  b2 = substr(kodony_mRNA, 2, 2),
  b3 = substr(kodony_mRNA, 3, 3),
  aa = genetic_code[kodony_mRNA],
  status = ifelse(kodony_mRNA %in% rozpoznawalne, "Found", "Missing")
)
df$b1 <- factor(df$b1, levels = rev(bazy))
df$b2 <- factor(df$b2, levels = bazy)
df$b3 <- factor(df$b3, levels = bazy)

counts <- phage_codons %>% 
  group_by(Codon) %>%
  summarize(Count = sum(Presence == 1))
  
df <- df %>%
  left_join(counts, by = c("codon" = "Codon")) %>%
  mutate(Count = coalesce(Count, 0))

df$label <- paste0(df$codon, "\n", df$Count)

ggplot(df, aes(x = b3, y = b1, fill = status)) +
  geom_tile(color = "black", width = 0.95, height = 0.95) +
  geom_text(aes(label = label), size = 6, fontface = "bold", lineheight = 0.8) +
  facet_wrap(~b2, nrow = 2, labeller = label_both) +
  scale_fill_manual(values = c("Found" = "lightblue", "Missing" = "pink")) +
  theme_minimal(base_size = 14) +
  theme(
    panel.grid = element_blank(),
    strip.text = element_text(size = 14),
    axis.title = element_text(size = 14),
    plot.title = element_text(hjust = 0.5)
  ) +
  labs(
    title = glue("Count map of amino acid hits in {species}"),
    x = "Third nucleotide", y = "First nucleotide", fill = "Status"
  )

```

</details>

<details>

<summary>tRNA gene homology in host</summary>

I extracted tRNA genes from [Ecoli genome](https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_000005845.2/) using ARAGORN which on average gives more hits in phages than tRNAscan-SE and output sequences are longer because of CCA tale that ARAGORN keeps (It should not change similarity score as it affects every sequence. Using Ecoli tRNA genes as query and target I pulled mmseq2 hits and plotted average identity between tRNA genes (Only a portion of all genes got mmseq2 hits, other sequences are too distinct) (21.07.2025)

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 16}
Ecoli_Ecoli_tRNA <- read_delim("~/Praktyki_Bioinfa/tRNA_Homology/Ecoli:Ecoli/results.m8", 
                      delim = "\t", escape_double = FALSE, 
                      col_names = c("query", "target", "pident", "alnlen", "mismatches", "gaps", "qstart", "qend","tstart","tend","evalue", "bitscore"), locale = locale(), 
                      trim_ws = TRUE)

Ecoli_Ecoli_tRNA$qAA <- paste(str_split_fixed(Ecoli_Ecoli_tRNA$query, "_", 6)[,4],toupper(str_split_fixed(Ecoli_Ecoli_tRNA$query, "_", 6)[,5]),sep = "_")
Ecoli_Ecoli_tRNA$tAA <- paste(str_split_fixed(Ecoli_Ecoli_tRNA$target, "_", 6)[,4],toupper(str_split_fixed(Ecoli_Ecoli_tRNA$target, "_", 6)[,5]),sep = "_")
Ecoli_Ecoli_tRNA$query <- sapply(str_split(Ecoli_Ecoli_tRNA$query, "_"), function(x) paste(head(x, 3), collapse = "_"))
Ecoli_Ecoli_tRNA$target <- sapply(str_split(Ecoli_Ecoli_tRNA$target, "_"), function(x) paste(head(x, 3), collapse = "_"))

{
Ecoli_Ecoli_tRNA_plot <- Ecoli_Ecoli_tRNA %>%
  group_by(qAA, tAA) %>%
  summarise(pident = mean(pident, na.rm = TRUE)) %>%
  ungroup()

ggplot(Ecoli_Ecoli_tRNA_plot, aes(x = qAA, y = tAA, fill = pident)) +
  geom_tile(color = "white", linewidth = 0.5) +
  geom_text(aes(label = round(pident, 2)),
            color = "white", size = 2, fontface = "bold") +
  scale_fill_gradientn(
    colors = c("#6D9EC1", "#E46726", "#B2182B"),
    name = "Avg pident",
    limits = c(min(Ecoli_Ecoli_tRNA_plot$pident, na.rm = TRUE), 
               max(Ecoli_Ecoli_tRNA_plot$pident, na.rm = TRUE))) +
  coord_fixed() +
  theme_minimal(base_size = 12) +
  labs(
    title = "Average sequences similarity between E. coli tRNA genes",
    subtitle = "Grouped by AA pairs and averaged",
    x = "E. coli tRNA (qAA)",
    y = "E. coli tRNA (tAA)"
  ) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
    axis.text.y = element_text(hjust = 1),
    panel.grid = element_blank(),
    plot.title = element_text(face = "bold", hjust = 0.5)
  )
}
```

</details>

<details>

<summary>tRNA gene homology host to virus</summary>

Using previous Ecoli tRNAs as query I pulled all tRNA genes from [Tequatrovirus](https://www.ncbi.nlm.nih.gov/datasets/genome/?taxon=10663) taxon (347 genomes) with accession number annotation to one target database and performed same mmseq2 search. Hits from this search have on average lover E-value (e-12 instead of e-48) which still can be considered as significant. (21.07.2025)

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 10,fig.width= 16}
Ecoli_Phage_tRNA <- read_delim("~/Praktyki_Bioinfa/tRNA_Homology/Ecoli:Tequatrovirus/results.m8", 
                      delim = "\t", escape_double = FALSE, 
                      col_names = c("query", "target", "pident", "alnlen", "mismatches", "gaps", "qstart", "qend","tstart","tend","evalue", "bitscore"), locale = locale(), 
                      trim_ws = TRUE)

Ecoli_Phage_tRNA$qAA <- paste(str_split_fixed(Ecoli_Phage_tRNA$query, "_", 6)[,4],toupper(str_split_fixed(Ecoli_Phage_tRNA$query, "_", 6)[,5]),sep = "_")
Ecoli_Phage_tRNA$tAA <- sapply(str_split(Ecoli_Phage_tRNA$target, "_"), function(x) paste(tail(head(x, -1), 2), collapse = "_"))
Ecoli_Phage_tRNA$query <- sapply(str_split(Ecoli_Phage_tRNA$query, "_"), function(x) paste(head(x, 3), collapse = "_"))
Ecoli_Phage_tRNA$target <- sapply(str_split(Ecoli_Phage_tRNA$target, "_"), function(x) paste(head(x, 2), collapse = "_"))

{
Ecoli_Phage_tRNA_plot <- Ecoli_Phage_tRNA %>%
  group_by(qAA, tAA) %>%
  summarise(pident = mean(pident, na.rm = TRUE)) %>%
  ungroup()

ggplot(Ecoli_Phage_tRNA_plot, aes(x = qAA, y = tAA, fill = pident)) +
  geom_tile(color = "white", linewidth = 0.5) +
  geom_text(aes(label = round(pident, 2)),
            color = "white", size = 7, fontface = "bold") +
  scale_fill_gradientn(
    colors = c("#6D9EC1", "#E46726", "#B2182B"),
    name = "Avg pident",
    limits = c(min(Ecoli_Phage_tRNA_plot$pident, na.rm = TRUE), 
               max(Ecoli_Phage_tRNA_plot$pident, na.rm = TRUE))) +
  coord_fixed() +
  theme_minimal(base_size = 12) +
  labs(
    title = "Average sequence similarity between phage and E. coli tRNA genes",
    subtitle = "Grouped by AA pairs and averaged",
    x = "E. coli tRNA (qAA)",
    y = "Phage tRNA (tAA)"
  ) +
  theme(
    axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
    axis.text.y = element_text(hjust = 1),
    panel.grid = element_blank(),
    plot.title = element_text(face = "bold", hjust = 0.5)
  )
}
```

In addition I performed two searches for E.coli non-parasitic phage taxons [Fromanvirus](https://www.ncbi.nlm.nih.gov/datasets/genome/?taxon=186764) and [Herelleviridae](https://www.ncbi.nlm.nih.gov/datasets/genome/?taxon=2560065) with respectively 332 and 549 viral genomes. After ARAGORN search Fromanvirus taxon showed little number of tRNA genes compared to Herelleviridae 277 hits and 3018 hits respectively. (06.08.2025)

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 8,fig.width= 16}
{
Ecoli_Tequatrovirus_tRNA <- read_delim("~/Praktyki_Bioinfa/tRNA_Homology/Ecoli:Tequatrovirus/results.m8", 
                               delim = "\t", escape_double = FALSE, 
                               col_names = c("query", "target", "pident", "alnlen", "mismatches", "gaps", "qstart", "qend","tstart","tend","evalue", "bitscore"), locale = locale(), 
                               trim_ws = TRUE)

Ecoli_Herelleviridae_tRNA <- read_delim("~/Praktyki_Bioinfa/tRNA_Homology/Ecoli:Herelleviridae/results.m8", 
                                       delim = "\t", escape_double = FALSE, 
                                       col_names = c("query", "target", "pident", "alnlen", "mismatches", "gaps", "qstart", "qend","tstart","tend","evalue", "bitscore"), locale = locale(), 
                                       trim_ws = TRUE)

Ecoli_Fromanvirus_tRNA <- read_delim("~/Praktyki_Bioinfa/tRNA_Homology/Ecoli:Fromanvirus/results.m8", 
                                        delim = "\t", escape_double = FALSE, 
                                        col_names = c("query", "target", "pident", "alnlen", "mismatches", "gaps", "qstart", "qend","tstart","tend","evalue", "bitscore"), locale = locale(), 
                                        trim_ws = TRUE)
}

#Query - bacteria, target - phage
{
Ecoli_Tequatrovirus_tRNA$qAA <- sapply(str_split(Ecoli_Tequatrovirus_tRNA$query, "_"), function(x) paste(toupper(tail(head(x, -2), 2)), collapse = "_"))
Ecoli_Tequatrovirus_tRNA$tAA <- sapply(str_split(Ecoli_Tequatrovirus_tRNA$target, "_"), function(x) paste(toupper(tail(head(x, -1), 2)), collapse = "_"))
Ecoli_Tequatrovirus_tRNA$query <- sapply(str_split(Ecoli_Tequatrovirus_tRNA$query, "_"), function(x) paste(head(x, 2), collapse = "_"))
Ecoli_Tequatrovirus_tRNA$target <- sapply(str_split(Ecoli_Tequatrovirus_tRNA$target, "_"), function(x) paste(head(x, 3), collapse = "_"))
  
Ecoli_Herelleviridae_tRNA$qAA <- sapply(str_split(Ecoli_Herelleviridae_tRNA$query, "_"), function(x) paste(toupper(tail(head(x, -2), 2)), collapse = "_"))
Ecoli_Herelleviridae_tRNA$tAA <- sapply(str_split(Ecoli_Herelleviridae_tRNA$target, "_"), function(x) paste(toupper(tail(head(x, -1), 2)), collapse = "_"))
Ecoli_Herelleviridae_tRNA$query <- sapply(str_split(Ecoli_Herelleviridae_tRNA$query, "_"), function(x) paste(head(x, 3), collapse = "_"))
Ecoli_Herelleviridae_tRNA$target <- sapply(str_split(Ecoli_Herelleviridae_tRNA$target, "_"), function(x) paste(head(x, 2), collapse = "_"))

Ecoli_Fromanvirus_tRNA$qAA <- sapply(str_split(Ecoli_Fromanvirus_tRNA$query, "_"), function(x) paste(toupper(tail(head(x, -2), 2)), collapse = "_"))
Ecoli_Fromanvirus_tRNA$tAA <- sapply(str_split(Ecoli_Fromanvirus_tRNA$target, "_"), function(x) paste(toupper(tail(head(x, -1), 2)), collapse = "_"))
Ecoli_Fromanvirus_tRNA$query <- sapply(str_split(Ecoli_Fromanvirus_tRNA$query, "_"), function(x) paste(head(x, 3), collapse = "_"))
Ecoli_Fromanvirus_tRNA$target <- sapply(str_split(Ecoli_Fromanvirus_tRNA$target, "_"), function(x) paste(head(x, 2), collapse = "_"))
}

{
Thiothrix_Tequatrovirus_tRNA <- read_delim("~/Praktyki_Bioinfa/tRNA_Homology/Thiothrix subterranea:Tequatrovirus/results.m8", 
                                         delim = "\t", escape_double = FALSE, 
                                         col_names = c("query", "target", "pident", "alnlen", "mismatches", "gaps", "qstart", "qend","tstart","tend","evalue", "bitscore"), locale = locale(), 
                                         trim_ws = TRUE)
Thiothrix_Herelleviridae_tRNA <- read_delim("~/Praktyki_Bioinfa/tRNA_Homology/Thiothrix subterranea:Herelleviridae/results.m8", 
                                         delim = "\t", escape_double = FALSE, 
                                         col_names = c("query", "target", "pident", "alnlen", "mismatches", "gaps", "qstart", "qend","tstart","tend","evalue", "bitscore"), locale = locale(), 
                                         trim_ws = TRUE)
Thiothrix_Fromanvirus_tRNA <- read_delim("~/Praktyki_Bioinfa/tRNA_Homology/Thiothrix subterranea:Fromanvirus/results.m8", 
                                         delim = "\t", escape_double = FALSE, 
                                         col_names = c("query", "target", "pident", "alnlen", "mismatches", "gaps", "qstart", "qend","tstart","tend","evalue", "bitscore"), locale = locale(), 
                                         trim_ws = TRUE)
}

{
  Thiothrix_Tequatrovirus_tRNA$qAA <- sapply(str_split(Thiothrix_Tequatrovirus_tRNA$query, "_"), function(x) paste(toupper(tail(head(x, -1), 2)), collapse = "_"))
  Thiothrix_Tequatrovirus_tRNA$tAA <- sapply(str_split(Thiothrix_Tequatrovirus_tRNA$target, "_"), function(x) paste(toupper(tail(head(x, -1), 2)), collapse = "_"))
  Thiothrix_Tequatrovirus_tRNA$query <- sapply(str_split(Thiothrix_Tequatrovirus_tRNA$query, "_"), function(x) paste(head(x, 3), collapse = "_"))
  Thiothrix_Tequatrovirus_tRNA$target <- sapply(str_split(Thiothrix_Tequatrovirus_tRNA$target, "_"), function(x) paste(head(x, 2), collapse = "_"))

  Thiothrix_Herelleviridae_tRNA$qAA <- sapply(str_split(Thiothrix_Herelleviridae_tRNA$query, "_"), function(x) paste(toupper(tail(head(x, -1), 2)), collapse = "_"))
  Thiothrix_Herelleviridae_tRNA$tAA <- sapply(str_split(Thiothrix_Herelleviridae_tRNA$target, "_"), function(x) paste(toupper(tail(head(x, -1), 2)), collapse = "_"))
  Thiothrix_Herelleviridae_tRNA$query <- sapply(str_split(Thiothrix_Herelleviridae_tRNA$query, "_"), function(x) paste(head(x, 3), collapse = "_"))
  Thiothrix_Herelleviridae_tRNA$target <- sapply(str_split(Thiothrix_Herelleviridae_tRNA$target, "_"), function(x) paste(head(x, 2), collapse = "_"))

  Thiothrix_Fromanvirus_tRNA$qAA <- sapply(str_split(Thiothrix_Fromanvirus_tRNA$query, "_"), function(x) paste(toupper(tail(head(x, -1), 2)), collapse = "_"))
  Thiothrix_Fromanvirus_tRNA$tAA <- sapply(str_split(Thiothrix_Fromanvirus_tRNA$target, "_"), function(x) paste(toupper(tail(head(x, -1), 2)), collapse = "_"))
  Thiothrix_Fromanvirus_tRNA$query <- sapply(str_split(Thiothrix_Fromanvirus_tRNA$query, "_"), function(x) paste(head(x, 3), collapse = "_"))
  Thiothrix_Fromanvirus_tRNA$target <- sapply(str_split(Thiothrix_Fromanvirus_tRNA$target, "_"), function(x) paste(head(x, 2), collapse = "_"))

}

library(gridExtra)
wykres_tRNA <- function(Ecoli_Phage_tRNA) {
  name <- deparse(substitute(Ecoli_Phage_tRNA))
  queryN <- sapply(str_split(name, "_"), function(x) paste(head(x, 1), collapse = "_"))
  targetN <- sapply(str_split(name, "_"), function(x) paste(head(tail(x,-1), 1), collapse = "_"))
  Ecoli_Phage_tRNA_plot <- Ecoli_Phage_tRNA %>%
    group_by(qAA, tAA) %>%
    summarise(pident = mean(pident, na.rm = TRUE)) %>%
    filter(pident > 0.85) %>%
    ungroup()
  
  ggplot(Ecoli_Phage_tRNA_plot, aes(x = qAA, y = tAA, fill = pident)) +
    geom_tile(color = "white", linewidth = 0.5) +
    geom_text(aes(label = round(pident, 2)),
              color = "white", size = 7, fontface = "bold") +
    scale_fill_gradientn(
      colors = c("#6D9EC1", "#E46726", "#B2182B"),
      name = "Avg pident",
      limits = c(min(Ecoli_Phage_tRNA_plot$pident, na.rm = TRUE), 
                 max(Ecoli_Phage_tRNA_plot$pident, na.rm = TRUE))) +
    coord_fixed() +
    theme_minimal(base_size = 12) +
    labs(
      title = glue("Average sequence similarity between {queryN} and {targetN}"),
      subtitle = "Grouped by AA pairs and averaged",
      x = glue("{queryN} tRNA (qAA)"),
      y = glue("{targetN} tRNA (tAA)")
    ) +
    theme(
      axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
      axis.text.y = element_text(hjust = 1),
      panel.grid = element_blank(),
      plot.title = element_text(face = "bold", hjust = 0.5)
    )
}

grid.arrange(
  wykres_tRNA(Ecoli_Tequatrovirus_tRNA),
  wykres_tRNA(Ecoli_Herelleviridae_tRNA),
  wykres_tRNA(Ecoli_Fromanvirus_tRNA),
  ncol = 3
  )

grid.arrange(
  wykres_tRNA(Thiothrix_Tequatrovirus_tRNA),
  wykres_tRNA(Thiothrix_Herelleviridae_tRNA),
  wykres_tRNA(Thiothrix_Fromanvirus_tRNA),
  ncol = 3
)
```

</details>

<details>

<summary>Phage tRNA bias check</summary>

In this part I will check how tRNA genes are distributed in specific phage taxons and in bigger class. I extracted tRNA genes using ARAGORN from [Tequatrovirus](https://www.ncbi.nlm.nih.gov/datasets/genome/?taxon=10663)(genus of 347 genomes), [Fromanvirus](https://www.ncbi.nlm.nih.gov/datasets/genome/?taxon=186764)(genus 333 genomes) and [Herelleviridae](https://www.ncbi.nlm.nih.gov/datasets/genome/?taxon=2560065)(family of 550 genomes). Additionally I extracted tRNA genes from whole [Caudoviricetes](https://www.ncbi.nlm.nih.gov/datasets/genome/?taxon=2731619&annotated_only=true&refseq_annotation=true&genbank_annotation=true&typical_only=true&exclude_mags=true&exclude_multi_isolates=true) taxon (from 07.2025, class of 14834 genomes).

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 20}
Caudoviricetes_tRNA_stat <- read_table("~/Praktyki_Bioinfa/Virus_tRNA/Virus_short_stat.txt", 
                                       col_names = c("Counts","Type"))
Herelleviridae_tRNA_stat <- read_table("~/Praktyki_Bioinfa/tRNA_Homology/Herelleviridae_stat.txt", 
                                                 col_names = c("Counts","Type"))
Fromanvirus_tRNA_stat <- read_table("~/Praktyki_Bioinfa/tRNA_Homology/Fromanvirus_stat.txt", 
                                         col_names = c("Counts","Type"))
Tequatrovirus_tRNA_stat <- read_table("~/Praktyki_Bioinfa/tRNA_Homology/Tequatrovirus_stat.txt", 
                                                 col_names = c("Counts","Type"))
{
# Delete unknown types
Caudoviricetes_tRNA_stat <- Caudoviricetes_tRNA_stat %>% 
  filter(Type != "unknown_NNN")
Herelleviridae_tRNA_stat <- Herelleviridae_tRNA_stat %>%
  filter(Type != "unknown_NNN")
Fromanvirus_tRNA_stat <- Fromanvirus_tRNA_stat %>%
  filter(Type != "unknown_NNN")
Tequatrovirus_tRNA_stat <- Tequatrovirus_tRNA_stat %>%
  filter(Type != "unknown_NNN")

Caudoviricetes_tRNA_stat <- Caudoviricetes_tRNA_stat %>%
  mutate(Normalized_Counts = Counts / 14834)
Herelleviridae_tRNA_stat <- Herelleviridae_tRNA_stat %>%
  mutate(Normalized_Counts = Counts / 550)
Fromanvirus_tRNA_stat <- Fromanvirus_tRNA_stat %>%
  mutate(Normalized_Counts = Counts / 333)
Tequatrovirus_tRNA_stat <- Tequatrovirus_tRNA_stat %>%
  mutate(Normalized_Counts = Counts / 347)
}
{
#tRNA codons from Herelleviridae type
Herelleviridae_codons <- Herelleviridae_tRNA_stat %>%
  pull(Type) %>% 
  str_split("_", simplify = TRUE) %>%
  .[, 2]
#tRNA codons from Fromanvirus type
Fromanvirus_codons <- Fromanvirus_tRNA_stat %>%
  pull(Type) %>% 
  str_split("_", simplify = TRUE) %>%
  .[, 2]
#tRNA codons from Tequatrovirus type
Tequatrovirus_codons <- Tequatrovirus_tRNA_stat %>%
  pull(Type) %>% 
  str_split("_", simplify = TRUE) %>%
  .[, 2]
}

H2 <- map_tRNA_to_mRNA(
  antykodony_DNA = Herelleviridae_codons,
  species = "Herelleviridae family",
  genetic_code = basic_code
)

F2 <- map_tRNA_to_mRNA(
  antykodony_DNA = Fromanvirus_codons,
  species = "Fromanvirus genus",
  genetic_code = basic_code
)

T2 <- map_tRNA_to_mRNA(
  antykodony_DNA = Tequatrovirus_codons,
  species = "Tequatrovirus genus",
  genetic_code = basic_code
)

grid.arrange(
  H2$wykres,
  F2$wykres,
  T2$wykres,
  ncol = 3
)

```

Important note: Distribution is normalized by overall number of genomes including there with no tRNA genes. This means that high presence of specific codon can be caused by overall popularity or by very high counts in rare genomes

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 20}
grid.arrange(
ggplot(Caudoviricetes_tRNA_stat, aes(x = reorder(Type,-Normalized_Counts), y = Normalized_Counts)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  labs(title = "Distribution of tRNA Types in Caudoviricetes",
       x = "tRNA Type",
       y = "Normalized Counts") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)),
ggplot(Herelleviridae_tRNA_stat, aes(x = reorder(Type,-Normalized_Counts), y = Normalized_Counts)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  labs(title = "Distribution of tRNA Types in Herelleviridae",
       x = "tRNA Type",
       y = "Normalized Counts") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)),
ggplot(Fromanvirus_tRNA_stat, aes(x = reorder(Type,-Normalized_Counts), y = Normalized_Counts)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  labs(title = "Distribution of tRNA Types in Fromanvirus",
       x = "tRNA Type",
       y = "Normalized Counts") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)),
ggplot(Tequatrovirus_tRNA_stat, aes(x = reorder(Type,-Normalized_Counts), y = Normalized_Counts)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  labs(title = "Distribution of tRNA Types in Tequatrovirus",
       x = "tRNA Type",
       y = "Normalized Counts") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)),
ncol = 2
)
```

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 20}
combined_df <- bind_rows(
  Caudoviricetes_tRNA_stat %>% mutate(Virus = "Caudoviricetes"),
  Herelleviridae_tRNA_stat %>% mutate(Virus = "Herelleviridae"),
  Fromanvirus_tRNA_stat %>% mutate(Virus = "Fromanvirus"), 
  Tequatrovirus_tRNA_stat %>% mutate(Virus = "Tequatrovirus")
)

combined_df <- combined_df %>%
  mutate(Normalized_Counts = ifelse(is.na(Normalized_Counts), 0, Normalized_Counts))

# Utwórz pełną siatkę wszystkich kombinacji Type i Virus
all_combinations <- expand.grid(
  Type = unique(combined_df$Type),
  Virus = unique(combined_df$Virus)
)

# Połącz z oryginalnymi danymi, wypełniając brakujące wartości 0
final_df <- all_combinations %>%
  left_join(combined_df, by = c("Type", "Virus")) %>%
  mutate(
    Normalized_Counts = ifelse(is.na(Normalized_Counts), 0, Normalized_Counts)
  ) %>%
  select(Type, Virus, Normalized_Counts, everything())

# Jeśli chcesz tylko Type i Normalised_Count:
final_simple <- final_df %>%
  select(Type, Virus, Normalized_Counts)

ggplot(final_simple, aes(x = Virus, y = Normalized_Counts, fill = Type)) +
  geom_bar(stat = "identity", position = "fill") +
  scale_y_continuous(labels = scales::percent) +
  labs(
    title = "Proportion of tRNA Types by Virus Family",
    x = "Virus Family",
    y = "Proportion",
    fill = "tRNA Type"
  ) +
  theme_minimal() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold"))

```

In this part I will explore possible clusters of tRNA genes in multidimensional data of 14834 Caudoviricetes genomes. I will use t-SNE and UMAP to visualize this data and see if there are any clusters of tRNA genes in this class. (23.08.2025)

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 20}
  library(dplyr)
  library(stringr)
  library(tidyr)
  library(ggplot2)
  library(DT)
  library(Rtsne)
  library(umap)
  library(ggforce)
  library(gridExtra)
  library(cluster)
  library(factoextra)
{ # For silhouette score and PCA
  
  # Wczytanie danych
  Type_Per_Acc <- read_table("~/Praktyki_Bioinfa/Virus_tRNA/Type_Per_Acc.txt", 
                             col_names = c("Count", "Acc"))
  
  # Przetwarzanie kolumn Acc i Type
  Type_Per_Acc$Type <- Type_Per_Acc %>%
    pull(Acc) %>% 
    str_split("_", simplify = TRUE) %>%
    .[, 3:4] %>% 
    apply(1, function(x) paste(x, collapse = "_"))
  
  Type_Per_Acc$Acc <- Type_Per_Acc %>%
    pull(Acc) %>% 
    str_split("_", simplify = TRUE) %>%
    .[, 1:2] %>% 
    apply(1, function(x) paste(x, collapse = "_"))
  
  # Filtrowanie unknown
  Type_Per_Acc <- Type_Per_Acc %>% 
    filter(Type != "unknown_NNN")
  
  # Wydobycie aminokwasu
  df <- Type_Per_Acc %>% 
    mutate(Amino = sub("_.*", "", Type))
  
  # Tworzymy macierz Acc × Type (counts)
  mat <- tidyr::pivot_wider(df, 
                            id_cols = Acc,  
                            names_from = Type, 
                            values_from = Count, 
                            values_fill = list(Count = 0))
  
  # Obliczenie najczęstszego (dominant) Amino per Acc
  amino_info <- df %>%
    group_by(Acc, Amino) %>%
    summarise(total_count = sum(Count), .groups = "drop") %>%
    group_by(Acc) %>%
    slice_max(total_count, n = 1, with_ties = FALSE) %>%
    ungroup() %>%
    select(Acc, Amino)
  
  # Normalizacja CPM (counts per million)
  mat_numeric <- mat %>%
    select(-Acc) %>%
    mutate(across(everything(), as.numeric)) %>%
    mutate(across(everything(), ~ . / sum(.) * 1e6))  # CPM normalization
  
  # Macierz z rownames
  mat_matrix <- as.matrix(mat_numeric)
  rownames(mat_matrix) <- mat$Acc
  
  # --- PCA jako preprocessing ---
  pca_result <- prcomp(mat_matrix, scale. = TRUE)
  pca_data <- pca_result$x[, 1:50]  # Użycie 50 głównych składowych
  
  # --- Przygotowanie danych bez duplikatów ---
  duplicate_rows <- duplicated(pca_data)
  mat_nodup <- pca_data[!duplicate_rows, , drop = FALSE]
  acc_nodup <- rownames(mat_matrix)[!duplicate_rows]
  
  # Filtrowanie amino_info dla nodup Acc
  amino_info_nodup <- amino_info %>% 
    filter(Acc %in% acc_nodup)
  
  # --- t-SNE (parametry dla 14800 genomów) ---
  set.seed(123)
  tsne <- Rtsne(mat_nodup, perplexity = 50, theta = 0.5, max_iter = 2000, check_duplicates = FALSE)
  
  # Data frame dla t-SNE
  tsne_df <- data.frame(
    Dim1 = tsne$Y[, 1],
    Dim2 = tsne$Y[, 2],
    Acc = acc_nodup
  )
  
  # Dołączanie Amino (najczęstszy per Acc)
  tsne_df <- tsne_df %>% 
    left_join(amino_info_nodup, by = "Acc")
  
  # --- UMAP (parametry dla 14800 genomów) ---
  set.seed(123)
  umap_result <- umap(mat_nodup, n_neighbors = 30, min_dist = 0.1, metric = "euclidean", n_epochs = 500)
  
  # Data frame dla UMAP
  umap_df <- data.frame(
    Dim1 = umap_result$layout[, 1],
    Dim2 = umap_result$layout[, 2],
    Acc = acc_nodup
  )
  
  # Dołączanie Amino (najczęstszy per Acc)
  umap_df <- umap_df %>% 
    left_join(amino_info_nodup, by = "Acc")
  
  # --- Wizualizacja Aminokwasy ---
  autoplot_tsne <- ggplot(tsne_df, aes(x = Dim1, y = Dim2, color = Amino)) +
    geom_point(size = 2, alpha = 0.6) +
    labs(title = "t-SNE tRNA profiles - Most common AA by Individual",
         x = "t-SNE Dimension 1", y = "t-SNE Dimension 2") +
    theme_minimal() +
    theme(legend.position = "right")
  
  autoplot_umap <- ggplot(umap_df, aes(x = Dim1, y = Dim2, color = Amino)) +
    geom_point(size = 2, alpha = 0.6) +
    labs(title = "UMAP tRNA profiles - Most common AA by Individual",
         x = "UMAP Dimension 1", y = "UMAP Dimension 2") +
    theme_minimal() +
    theme(legend.position = "right")
  
  # --- Klasteryzacja (optymalizacja k) ---
  # Użycie silhouette score do wyboru k (do 10)
  set.seed(123)
  sil_tsne <- c()
  for(k in 2:10){
    kmeans_k <- kmeans(tsne$Y, centers = k, nstart = 50)
    ss <- silhouette(kmeans_k$cluster, dist(tsne$Y))
    sil_tsne[k-1] <- mean(ss[, 3])
  }
  optimal_k_tsne <- which.max(sil_tsne) + 1
  
  set.seed(123)
  sil_umap <- c()
  for(k in 2:10){
    kmeans_k <- kmeans(umap_result$layout, centers = k, nstart = 50)
    ss <- silhouette(kmeans_k$cluster, dist(umap_result$layout))
    sil_umap[k-1] <- mean(ss[, 3])
  }
  optimal_k_umap <- which.max(sil_umap) + 1
  
  # Klasteryzacja dla t-SNE
  set.seed(123)
  kmeans_tsne <- kmeans(tsne$Y, centers = optimal_k_tsne, nstart = 50)
  tsne_df$Cluster <- as.factor(kmeans_tsne$cluster)
  
  # Klasteryzacja dla UMAP
  set.seed(123)
  kmeans_umap <- kmeans(umap_result$layout, centers = optimal_k_umap, nstart = 50)
  umap_df$Cluster <- as.factor(kmeans_umap$cluster)
  
  # Wizualizacja z klasterami
  tsne_cluster_plot <- ggplot(tsne_df, aes(x = Dim1, y = Dim2, color = Cluster)) +
    geom_point(size = 2, alpha = 0.6) +
    geom_mark_hull(aes(fill = Cluster), alpha = 0.2, show.legend = FALSE) +
    labs(title = "t-SNE kmeans Clustering", x = "t-SNE Dimension 1", y = "t-SNE Dimension 2") +
    theme_minimal()
  
  umap_cluster_plot <- ggplot(umap_df, aes(x = Dim1, y = Dim2, color = Cluster)) +
    geom_point(size = 2, alpha = 0.6) +
    geom_mark_hull(aes(fill = Cluster), alpha = 0.2, show.legend = FALSE) +
    labs(title = "UMAP kmeans Clustering", x = "UMAP Dimension 1", y = "UMAP Dimension 2") +
    theme_minimal()
  
  # --- Porównanie Aminokwasy vs Klastery ---
  gridExtra::grid.arrange(
    autoplot_tsne,
    tsne_cluster_plot,
    autoplot_umap,
    umap_cluster_plot,
    ncol = 2, nrow = 2
  )}

{# Wczytanie danych
Type_Per_Acc <- read_table("~/Praktyki_Bioinfa/Virus_tRNA/Type_Per_Acc.txt", 
                           col_names = c("Count", "Acc"))

# Przetwarzanie kolumn Acc i Type (zachowujemy pełny Type jako Amino_Codon)
Type_Per_Acc$Type <- Type_Per_Acc %>%
  pull(Acc) %>% 
  str_split("_", simplify = TRUE) %>%
  .[, 3:4] %>% 
  apply(1, function(x) paste(x, collapse = "_"))  # Pełny Type (np. Met_CAT)

Type_Per_Acc$Acc <- Type_Per_Acc %>%
  pull(Acc) %>% 
  str_split("_", simplify = TRUE) %>%
  .[, 1:2] %>% 
  apply(1, function(x) paste(x, collapse = "_"))

# Filtrowanie unknown
Type_Per_Acc <- Type_Per_Acc %>% 
  filter(Type != "unknown_NNN")

# Wydobycie pełnego Type jako identyfikator (Amino_Codon)
df <- Type_Per_Acc %>% 
  mutate(Type_Full = Type)  # Zachowujemy pełny Type jako Amino_Codon

# Tworzymy macierz Type_Full × Acc (counts)
mat <- tidyr::pivot_wider(df, 
                          id_cols = Type_Full,  # Używamy pełnego Type (Amino_Codon)
                          names_from = Acc, 
                          values_from = Count, 
                          values_fill = list(Count = 0))

# Zachowujemy informację o Type_Full dla kolorowania
amino_info <- df %>% 
  select(Type_Full) %>% 
  distinct() %>%
  mutate(Amino = sub("_.*", "", Type_Full))  # Ekstrakcja tylko aminokwasu dla kolorowania

# Normalizacja CPM (counts per million)
mat_numeric <- mat %>%
  select(-Type_Full) %>%
  mutate(across(everything(), as.numeric)) %>%
  mutate(across(everything(), ~ . / sum(.) * 1e6))  # CPM normalization

# Macierz z rownames
mat_matrix <- as.matrix(mat_numeric)
rownames(mat_matrix) <- mat$Type_Full

# --- PCA jako preprocessing ---
pca_result <- prcomp(mat_matrix, scale. = TRUE)
pca_data <- pca_result$x[, 1:20]  # Użycie 20 głównych składowych (mniejszy zbiór)

# --- Przygotowanie danych bez duplikatów ---
duplicate_rows <- duplicated(pca_data)
mat_nodup <- pca_data[!duplicate_rows, , drop = FALSE]
type_nodup <- rownames(mat_matrix)[!duplicate_rows]

# Filtrowanie amino_info dla nodup Type
amino_info_nodup <- amino_info %>% 
  filter(Type_Full %in% type_nodup)

# --- t-SNE (parametry dla ~60-100 Type_Full) ---
set.seed(123)
tsne <- Rtsne(mat_nodup, perplexity = 20, theta = 0.5, max_iter = 1000, check_duplicates = FALSE)

# Data frame dla t-SNE
tsne_df <- data.frame(
  Dim1 = tsne$Y[, 1],
  Dim2 = tsne$Y[, 2],
  Type_Full = type_nodup
)

# Dołączanie Amino dla kolorowania
tsne_df <- tsne_df %>% 
  left_join(amino_info_nodup, by = "Type_Full")

# --- UMAP (parametry dla ~60-100 Type_Full) ---
set.seed(123)
umap_result <- umap(mat_nodup, n_neighbors = 15, min_dist = 0.05, metric = "cosine", n_epochs = 500)  # Użycie cosine dla sparse data

# Data frame dla UMAP
umap_df <- data.frame(
  Dim1 = umap_result$layout[, 1],
  Dim2 = umap_result$layout[, 2],
  Type_Full = type_nodup
)

# Dołączanie Amino dla kolorowania
umap_df <- umap_df %>% 
  left_join(amino_info_nodup, by = "Type_Full")

# --- Wizualizacja Aminokwasy ---
autoplot_tsne <- ggplot(tsne_df, aes(x = Dim1, y = Dim2, color = Amino)) +
  geom_point(size = 3, alpha = 0.7) +
  labs(title = "t-SNE tRNA profiles - tRNA Types with codons",
       x = "t-SNE Dimension 1", y = "t-SNE Dimension 2") +
  theme_minimal() +
  theme(legend.position = "right")

autoplot_umap <- ggplot(umap_df, aes(x = Dim1, y = Dim2, color = Amino)) +
  geom_point(size = 3, alpha = 0.7) +
  labs(title = "UMAP tRNA profiles - tRNA Types with codons",
       x = "UMAP Dimension 1", y = "UMAP Dimension 2") +
  theme_minimal() +
  theme(legend.position = "right")

# --- Klasteryzacja (optymalizacja k) ---
# Użycie silhouette score do wyboru k (do 10)
set.seed(123)
sil_tsne <- c()
for(k in 2:10){
  kmeans_k <- kmeans(tsne$Y, centers = k, nstart = 50)
  ss <- silhouette(kmeans_k$cluster, dist(tsne$Y))
  sil_tsne[k-1] <- mean(ss[, 3])
}
optimal_k_tsne <- which.max(sil_tsne) + 1

set.seed(123)
sil_umap <- c()
for(k in 2:10){
  kmeans_k <- kmeans(umap_result$layout, centers = k, nstart = 50)
  ss <- silhouette(kmeans_k$cluster, dist(umap_result$layout))
  sil_umap[k-1] <- mean(ss[, 3])
}
optimal_k_umap <- which.max(sil_umap) + 1

# Klasteryzacja dla t-SNE
set.seed(123)
kmeans_tsne <- kmeans(tsne$Y, centers = optimal_k_tsne, nstart = 50)
tsne_df$Cluster <- as.factor(kmeans_tsne$cluster)

# Klasteryzacja dla UMAP
set.seed(123)
kmeans_umap <- kmeans(umap_result$layout, centers = optimal_k_umap, nstart = 50)
umap_df$Cluster <- as.factor(kmeans_umap$cluster)

# Wizualizacja z klasterami
tsne_cluster_plot <- ggplot(tsne_df, aes(x = Dim1, y = Dim2, color = Cluster)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_mark_hull(aes(fill = Cluster), alpha = 0.2, show.legend = FALSE) +
  labs(title = "t-SNE kmeans Clusters", x = "t-SNE Dimension 1", y = "t-SNE Dimension 2") +
  theme_minimal()

umap_cluster_plot <- ggplot(umap_df, aes(x = Dim1, y = Dim2, color = Cluster)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_mark_hull(aes(fill = Cluster), alpha = 0.2, show.legend = FALSE) +
  labs(title = "UMAP kmeans Clusters", x = "UMAP Dimension 1", y = "UMAP Dimension 2") +
  theme_minimal()

# --- Porównanie Aminokwasy vs Klastery ---
gridExtra::grid.arrange(
  autoplot_tsne,
  tsne_cluster_plot,
  autoplot_umap,
  umap_cluster_plot,
  ncol = 2, nrow = 2
)

# --- Ekstrakcja pełnych typów tRNA dla każdego klastra ---
# Dla t-SNE
tsne_clusters_full <- tsne_df %>%
  group_by(Cluster) %>%
  summarise(
    Amino = toString(unique(Amino)),  # Łączenie unikalnych aminokwasów w jeden ciąg
    Type_Full = toString(unique(Type_Full))  # Łączenie unikalnych typów w jeden ciąg
  )

# Dla UMAP
umap_clusters_full <- umap_df %>%
  group_by(Cluster) %>%
  summarise(
    Amino = toString(unique(Amino)),  # Łączenie unikalnych aminokwasów w jeden ciąg
    Type_Full = toString(unique(Type_Full))  # Łączenie unikalnych typów w jeden ciąg
  )

tsne_dt <- datatable(
  tsne_clusters_full,
  rownames = FALSE,
  colnames = c("Cluster", "Amino Acids", "Types (AA_Codon)"),
  caption = htmltools::tags$caption(
    style = 'caption-side: top; text-align: left;',
    "t-SNE Clusters with tRNA Types"
  ),
  filter = "top",  # umożliwia wyszukiwanie po każdej kolumnie
  options = list(
    pageLength = 10,
    autoWidth = TRUE,
    searchHighlight = TRUE
  )
)

# --- UMAP ---
umap_dt <- datatable(
  umap_clusters_full,
  rownames = FALSE,
  colnames = c("Cluster", "Amino Acids", "Types (AA_Codon)"),
  caption = htmltools::tags$caption(
    style = 'caption-side: top; text-align: left;',
    "UMAP Clusters with tRNA Types"
  ),
  filter = "top",  # wyszukiwanie/filtrowanie
  options = list(
    pageLength = 10,
    autoWidth = TRUE,
    searchHighlight = TRUE
  )
)
}
```

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 20}
{
Type_Per_Acc <- read_table("~/Praktyki_Bioinfa/Virus_tRNA/Type_Per_Acc.txt", 
                           col_names = c("Count", "Acc"))

# Przetwarzanie kolumn Acc i Type
Type_Per_Acc$Type <- Type_Per_Acc %>%
  pull(Acc) %>% 
  str_split("_", simplify = TRUE) %>%
  .[, 3:4] %>% 
  apply(1, function(x) paste(x, collapse = "_"))

Type_Per_Acc$Acc <- Type_Per_Acc %>%
  pull(Acc) %>% 
  str_split("_", simplify = TRUE) %>%
  .[, 1:2] %>% 
  apply(1, function(x) paste(x, collapse = "_"))

# Filtrowanie unknown
Type_Per_Acc <- Type_Per_Acc %>% 
  filter(Type != "unknown_NNN")

# Obliczenie całkowitej liczby tRNA (Total tRNA Count) na genom
tRNA_total <- Type_Per_Acc %>%
  group_by(Acc) %>%
  summarise(Total_tRNA = sum(Count), .groups = "drop")

# Przygotowanie danych do klasteryzacji i wizualizacji
data_tRNA <- tRNA_total

# --- t-SNE na podstawie Total_tRNA ---
set.seed(123)
tsne_tRNA <- Rtsne(as.matrix(data_tRNA$Total_tRNA), perplexity = 50, theta = 0.5, max_iter = 2000, check_duplicates = FALSE)

# Data frame dla t-SNE
tsne_df <- data.frame(
  Dim1 = tsne_tRNA$Y[, 1],
  Dim2 = tsne_tRNA$Y[, 2],
  Acc = data_tRNA$Acc,
  Total_tRNA = data_tRNA$Total_tRNA
)

# --- UMAP na podstawie Total_tRNA ---
set.seed(123)
umap_tRNA <- umap(as.matrix(data_tRNA$Total_tRNA), n_neighbors = 30, min_dist = 0.1, metric = "euclidean", n_epochs = 500)

# Data frame dla UMAP
umap_df <- data.frame(
  Dim1 = umap_tRNA$layout[, 1],
  Dim2 = umap_tRNA$layout[, 2],
  Acc = data_tRNA$Acc,
  Total_tRNA = data_tRNA$Total_tRNA
)

# --- Wizualizacja Total tRNA ---
autoplot_tsne <- ggplot(tsne_df, aes(x = Dim1, y = Dim2, color = Total_tRNA)) +
  geom_point(size = 2, alpha = 0.6) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "t-SNE tRNA profiles - Total tRNA Count by Individual",
       x = "t-SNE Dimension 1", y = "t-SNE Dimension 2") +
  theme_minimal() +
  theme(legend.position = "right")

autoplot_umap <- ggplot(umap_df, aes(x = Dim1, y = Dim2, color = Total_tRNA)) +
  geom_point(size = 2, alpha = 0.6) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "UMAP tRNA profiles - Total tRNA Count by Individual",
       x = "UMAP Dimension 1", y = "UMAP Dimension 2") +
  theme_minimal() +
  theme(legend.position = "right")

# --- Klasteryzacja na podstawie Total_tRNA ---
# Użycie silhouette score do wyboru k (do 10)
set.seed(123)
sil_tsne <- c()
for(k in 2:10){
  kmeans_k <- kmeans(tsne_tRNA$Y, centers = k, nstart = 50)
  ss <- silhouette(kmeans_k$cluster, dist(tsne_tRNA$Y))
  sil_tsne[k-1] <- mean(ss[, 3])
}
optimal_k_tsne <- which.max(sil_tsne) + 1

set.seed(123)
sil_umap <- c()
for(k in 2:10){
  kmeans_k <- kmeans(umap_tRNA$layout, centers = k, nstart = 50)
  ss <- silhouette(kmeans_k$cluster, dist(umap_tRNA$layout))
  sil_umap[k-1] <- mean(ss[, 3])
}
optimal_k_umap <- which.max(sil_umap) + 1

# Klasteryzacja dla t-SNE
set.seed(123)
kmeans_tsne <- kmeans(tsne_tRNA$Y, centers = optimal_k_tsne, nstart = 50)
tsne_df$Cluster <- as.factor(kmeans_tsne$cluster)

# Klasteryzacja dla UMAP
set.seed(123)
kmeans_umap <- kmeans(umap_tRNA$layout, centers = optimal_k_umap, nstart = 50)
umap_df$Cluster <- as.factor(kmeans_umap$cluster)

# Wizualizacja z klasterami
tsne_cluster_plot <- ggplot(tsne_df, aes(x = Dim1, y = Dim2, color = Cluster)) +
  geom_point(size = 2, alpha = 0.6) +
  geom_mark_hull(aes(fill = Cluster), alpha = 0.2, show.legend = FALSE) +
  labs(title = "t-SNE kmeans Clustering", x = "t-SNE Dimension 1", y = "t-SNE Dimension 2") +
  theme_minimal()

umap_cluster_plot <- ggplot(umap_df, aes(x = Dim1, y = Dim2, color = Cluster)) +
  geom_point(size = 2, alpha = 0.6) +
  geom_mark_hull(aes(fill = Cluster), alpha = 0.2, show.legend = FALSE) +
  labs(title = "UMAP kmeans Clustering", x = "UMAP Dimension 1", y = "UMAP Dimension 2") +
  theme_minimal()

gridExtra::grid.arrange(
  autoplot_tsne,
  tsne_cluster_plot,
  autoplot_umap,
  umap_cluster_plot,
  ncol = 2, nrow = 2)
}
```

<details>

<summary>Tables for tRNA clusters</summary>

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 20}
tsne_dt
umap_dt
```

</details>

</details>

</details>

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

<details>

<summary>Exploring tRNA profile clusters</summary>

Inspired by phage clustered profiles from previous step, I will explore the tRNA profiles of: 1. Selected prokariotic taxons with alternative codon usage 2. UHGG database (4905 genomes of human gut microbiome)

<details>

<summary>tRNA profiles for chosen alternatively coded taxons</summary>

I chose to include three types of alternative genomes:

1.  [Mycoplasmoidaceae](https://www.ncbi.nlm.nih.gov/datasets/genome/?taxon=2790998&annotated_only=true&refseq_annotation=true&genbank_annotation=true) - type 4 (511 genomes)

2.  [Acholeplasmataceae](https://www.ncbi.nlm.nih.gov/datasets/genome/?taxon=2146&annotated_only=true&refseq_annotation=true&genbank_annotation=true) - type 11 (445 genomes)

3.  [Candidatus Altimarinota and Candidatus Absconditibacteriota](https://www.ncbi.nlm.nih.gov/datasets/genome/?taxon=363464,221235&annotated_only=true&refseq_annotation=true&genbank_annotation=true) - type 25 (216 genomes)

    I used ARAGORN to extract tRNA genes with specific parameter for alternative genomes (24.08.2025)

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 20}
df4 <- read_table("~/Praktyki_Bioinfa/tRNA_by_type/Type_Per_Acc_4_fixed.txt", 
                  col_names = c("Count", "Acc")) %>%
  mutate(Group = 4)

df11 <- read_table("~/Praktyki_Bioinfa/tRNA_by_type/Type_Per_Acc_11_fixed.txt", 
                   col_names = c("Count", "Acc")) %>%
  mutate(Group = 11)

df25 <- read_table("~/Praktyki_Bioinfa/tRNA_by_type/Type_Per_Acc_25_fixed.txt", 
                   col_names = c("Count", "Acc")) %>%
  mutate(Group = 25)

# Połączenie w jeden DF
Type_Per_Acc <- bind_rows(df4, df11, df25)

# Przetwarzanie kolumn Acc i Type
Type_Per_Acc$Type <- Type_Per_Acc %>%
  pull(Acc) %>% 
  str_split("_", simplify = TRUE) %>%
  .[, 3:4] %>% 
  apply(1, function(x) paste(x, collapse = "_"))

Type_Per_Acc$Acc <- Type_Per_Acc %>%
  pull(Acc) %>% 
  str_split("_", simplify = TRUE) %>%
  .[, 1:2] %>% 
  apply(1, function(x) paste(x, collapse = "_"))

# Filtrowanie unknown
Type_Per_Acc <- Type_Per_Acc %>% 
  filter(Type != "unknown_NNN")

# Wydobycie aminokwasu
df <- Type_Per_Acc %>% 
  mutate(Amino = sub("_.*", "", Type))

# Etap 1
{
# Tworzymy macierz Acc × Type (counts)
mat <- tidyr::pivot_wider(df, 
                          id_cols = Acc,
                          names_from = Type, 
                          values_from = Count, 
                          values_fill = list(Count = 0))

# Obliczenie całkowitej liczby tRNA (Total tRNA Count) na genom
tRNA_total <- df %>%
  group_by(Acc) %>%
  summarise(Total_tRNA = sum(Count), .groups = "drop")

# Obliczenie najczęstszego kodonu (Most Frequent Codon) na genom
codon_info <- df %>%
  group_by(Acc, Type) %>%
  summarise(total_count = sum(Count), .groups = "drop") %>%
  group_by(Acc) %>%
  slice_max(total_count, n = 1, with_ties = FALSE) %>%
  ungroup() %>%
  select(Acc, Type) %>%
  dplyr::rename(Most_Frequent_Codon = Type)

# Dołączanie Group do mat
mat <- mat %>%
  left_join(df %>% select(Acc, Group) %>% distinct(), by = "Acc")

# Normalizacja CPM (counts per million)
mat_numeric <- mat %>%
  select(-Acc, -Group) %>%
  mutate(across(everything(), as.numeric)) %>%
  mutate(across(everything(), ~ . / sum(.) * 1e6))  # CPM normalization

# Macierz z rownames
mat_matrix <- as.matrix(mat_numeric)
rownames(mat_matrix) <- mat$Acc

# --- PCA jako preprocessing - dostosowane do mniejszego zbioru ---
pca_result <- prcomp(mat_matrix, scale. = TRUE)
pca_data <- pca_result$x[, 1:30]  # Zmniejszono do 30 komponentów dla ~1300 genomów

# --- Przygotowanie danych bez duplikatów ---
duplicate_rows <- duplicated(pca_data)
mat_nodup <- pca_data[!duplicate_rows, , drop = FALSE]
acc_nodup <- rownames(mat_matrix)[!duplicate_rows]

# Filtrowanie metadanych dla nodup Acc
mat_nodup_df <- mat %>% filter(Acc %in% acc_nodup)
tRNA_total_nodup <- tRNA_total %>% filter(Acc %in% acc_nodup)
codon_info_nodup <- codon_info %>% filter(Acc %in% acc_nodup)

# --- t-SNE (parametry dostosowane do ~1300 genomów) ---
set.seed(123)
tsne <- Rtsne(mat_nodup, perplexity = 30, theta = 0.5, max_iter = 1000, check_duplicates = FALSE)  # Zmniejszono perplexity

# Data frame dla t-SNE
tsne_df <- data.frame(
  Dim1 = tsne$Y[, 1],
  Dim2 = tsne$Y[, 2],
  Acc = acc_nodup,
  Group = mat_nodup_df$Group,
  Total_tRNA = tRNA_total_nodup$Total_tRNA,
  Most_Frequent_Codon = codon_info_nodup$Most_Frequent_Codon
)

# --- UMAP (parametry dostosowane do ~1300 genomów) ---
set.seed(123)
umap_result <- umap(mat_nodup, n_neighbors = 20, min_dist = 0.1, metric = "euclidean", n_epochs = 300)  # Zmniejszono n_neighbors, zmniejszono n_epochs

# Data frame dla UMAP
umap_df <- data.frame(
  Dim1 = umap_result$layout[, 1],
  Dim2 = umap_result$layout[, 2],
  Acc = acc_nodup,
  Group = mat_nodup_df$Group,
  Total_tRNA = tRNA_total_nodup$Total_tRNA,
  Most_Frequent_Codon = codon_info_nodup$Most_Frequent_Codon
)

# --- Wizualizacja po Group (Etap 1) ---
autoplot_tsne_group <- ggplot(tsne_df, aes(x = Dim1, y = Dim2, color = as.factor(Group))) +
  geom_point(size = 2, alpha = 0.6) +
  labs(title = "t-SNE tRNA profiles - Alternative genomes types (4,11,25)",
       x = "t-SNE Dimension 1", y = "t-SNE Dimension 2") +
  theme_minimal() +
  theme(legend.position = "right")

autoplot_umap_group <- ggplot(umap_df, aes(x = Dim1, y = Dim2, color = as.factor(Group))) +
  geom_point(size = 2, alpha = 0.6) +
  labs(title = "UMAP tRNA profiles - Alternative genomes types (4,11,25)",
       x = "UMAP Dimension 1", y = "UMAP Dimension 2") +
  theme_minimal() +
  theme(legend.position = "right")

# --- Dodatkowy wykres grupujący po ogólnej liczbie genów tRNA ---
autoplot_tsne_total_tRNA <- ggplot(tsne_df, aes(x = Dim1, y = Dim2, color = Total_tRNA)) +
  geom_point(size = 2, alpha = 0.6) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "t-SNE tRNA profiles - by Total tRNA Count",
       x = "t-SNE Dimension 1", y = "t-SNE Dimension 2") +
  theme_minimal() +
  theme(legend.position = "right")

autoplot_umap_total_tRNA <- ggplot(umap_df, aes(x = Dim1, y = Dim2, color = Total_tRNA)) +
  geom_point(size = 2, alpha = 0.6) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "UMAP tRNA profiles - by Total tRNA Count",
       x = "UMAP Dimension 1", y = "UMAP Dimension 2") +
  theme_minimal() +
  theme(legend.position = "right")

# --- Klasteryzacja (na podstawie profilów tRNA, Etap 1) ---
set.seed(123)
sil_tsne <- c()
for(k in 2:4){
  kmeans_k <- kmeans(tsne$Y, centers = k, nstart = 50)
  ss <- silhouette(kmeans_k$cluster, dist(tsne$Y))
  sil_tsne[k-1] <- mean(ss[, 3])
}
optimal_k_tsne <- which.max(sil_tsne) + 1

set.seed(123)
sil_umap <- c()
for(k in 2:4){
  kmeans_k <- kmeans(umap_result$layout, centers = k, nstart = 50)
  ss <- silhouette(kmeans_k$cluster, dist(umap_result$layout))
  sil_umap[k-1] <- mean(ss[, 3])
}
optimal_k_umap <- which.max(sil_umap) + 1

# Klasteryzacja dla t-SNE
set.seed(123)
kmeans_tsne <- kmeans(tsne$Y, centers = optimal_k_tsne, nstart = 50)
tsne_df$Cluster <- as.factor(kmeans_tsne$cluster)

# Klasteryzacja dla UMAP
set.seed(123)
kmeans_umap <- kmeans(umap_result$layout, centers = optimal_k_umap, nstart = 50)
umap_df$Cluster <- as.factor(kmeans_umap$cluster)

# Wizualizacja z klasterami (Etap 1)
tsne_cluster_plot <- ggplot(tsne_df, aes(x = Dim1, y = Dim2, color = Cluster)) +
  geom_point(size = 2, alpha = 0.6) +
  geom_mark_hull(aes(fill = Cluster), alpha = 0.2, show.legend = FALSE) +
  labs(title = "t-SNE kmeans Clustering", x = "t-SNE Dimension 1", y = "t-SNE Dimension 2") +
  theme_minimal()

umap_cluster_plot <- ggplot(umap_df, aes(x = Dim1, y = Dim2, color = Cluster)) +
  geom_point(size = 2, alpha = 0.6) +
  geom_mark_hull(aes(fill = Cluster), alpha = 0.2, show.legend = FALSE) +
  labs(title = "UMAP kmeans Clustering", x = "UMAP Dimension 1", y = "UMAP Dimension 2") +
  theme_minimal()

gridExtra::grid.arrange(
  autoplot_tsne_group,
  autoplot_tsne_total_tRNA,
  tsne_cluster_plot,
  autoplot_umap_group,
  autoplot_umap_total_tRNA,
  umap_cluster_plot,
  ncol = 3, nrow = 2
)
}
```

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 20}
# Etap 2
{
# Tworzymy macierz Type × Acc (counts)
mat_type <- tidyr::pivot_wider(df, 
                               id_cols = Type, 
                               names_from = Acc, 
                               values_from = Count, 
                               values_fill = list(Count = 0))

# Zachowujemy informację o Amino dla kolorowania
amino_info_type <- df %>% 
  select(Type, Amino) %>% 
  distinct()

# Dołączanie Group do mat_type (agregacja po Acc)
group_info_type <- df %>%
  group_by(Type) %>%
  summarise(Group = toString(unique(Group))) %>%
  unnest(cols = c(Group))  # Rozdzielanie, jeśli są różne wartości

# Połączenie z macierzą typów (dla wizualizacji)
mat_type <- mat_type %>%
  left_join(group_info_type, by = "Type") %>%
  left_join(amino_info_type, by = "Type")

# Normalizacja CPM (counts per million)
mat_numeric_type <- mat_type %>%
  select(-Type, -Group, -Amino) %>%
  mutate(across(everything(), as.numeric)) %>%
  mutate(across(everything(), ~ . / sum(.) * 1e6))  # CPM normalization

# Macierz z rownames
mat_matrix_type <- as.matrix(mat_numeric_type)
rownames(mat_matrix_type) <- mat_type$Type

# --- PCA jako preprocessing ---
pca_result_type <- prcomp(mat_matrix_type, scale. = TRUE)
pca_data_type <- pca_result_type$x[, 1:15]  # Zmniejszono do 15 komponentów dla mniejszego zbioru

# --- Przygotowanie danych bez duplikatów ---
duplicate_rows_type <- duplicated(pca_data_type)
mat_nodup_type <- pca_data_type[!duplicate_rows_type, , drop = FALSE]
type_nodup <- rownames(mat_matrix_type)[!duplicate_rows_type]

# Filtrowanie metadanych dla nodup Type
mat_nodup_type_df <- mat_type %>% filter(Type %in% type_nodup)

# --- t-SNE dla typów ---
set.seed(123)
tsne_type <- Rtsne(mat_nodup_type, perplexity = 15, theta = 0.5, max_iter = 800, check_duplicates = FALSE)  # Zmniejszono perplexity i iteracje

# Data frame dla t-SNE
tsne_df_type <- data.frame(
  Dim1 = tsne_type$Y[, 1],
  Dim2 = tsne_type$Y[, 2],
  Type = type_nodup,
  Group = mat_nodup_type_df$Group,
  Amino = mat_nodup_type_df$Amino
)

# --- UMAP dla typów ---
set.seed(123)
umap_type <- umap(mat_nodup_type, n_neighbors = 10, min_dist = 0.05, metric = "cosine", n_epochs = 150)  # Zmniejszono n_neighbors i n_epochs

# Data frame dla UMAP
umap_df_type <- data.frame(
  Dim1 = umap_type$layout[, 1],
  Dim2 = umap_type$layout[, 2],
  Type = type_nodup,
  Group = mat_nodup_type_df$Group,
  Amino = mat_nodup_type_df$Amino
)

# --- Wizualizacja po Group (Etap 2) ---
autoplot_tsne_group_type <- ggplot(tsne_df_type, aes(x = Dim1, y = Dim2, color = Group)) +
  geom_point(size = 3, alpha = 0.7) +
  labs(title = "t-SNE tRNA profiles - Alternative genomes types and combinations(4,11,25)",
       x = "t-SNE Dimension 1", y = "t-SNE Dimension 2") +
  theme_minimal() +
  theme(legend.position = "right")

autoplot_umap_group_type <- ggplot(umap_df_type, aes(x = Dim1, y = Dim2, color = Group)) +
  geom_point(size = 3, alpha = 0.7) +
  labs(title = "UMAP tRNA profiles - Alternative genomes types and combinations(4,11,25)",
       x = "UMAP Dimension 1", y = "UMAP Dimension 2") +
  theme_minimal() +
  theme(legend.position = "right")

# --- Wizualizacja po Amino (Etap 2) ---
autoplot_tsne_amino_type <- ggplot(tsne_df_type, aes(x = Dim1, y = Dim2, color = Amino)) +
  geom_point(size = 3, alpha = 0.7) +
  labs(title = "t-SNE tRNA profiles - by Amino-codon types",
       x = "t-SNE Dimension 1", y = "t-SNE Dimension 2") +
  theme_minimal() +
  theme(legend.position = "right")

autoplot_umap_amino_type <- ggplot(umap_df_type, aes(x = Dim1, y = Dim2, color = Amino)) +
  geom_point(size = 3, alpha = 0.7) +
  labs(title = "UMAP tRNA profiles - by Amino-codon types",
       x = "UMAP Dimension 1", y = "UMAP Dimension 2") +
  theme_minimal() +
  theme(legend.position = "right")

# --- Klasteryzacja dla typów ---
set.seed(123)
sil_tsne_type <- c()
for(k in 2:4){
  kmeans_k <- kmeans(tsne_type$Y, centers = k, nstart = 50)
  ss <- silhouette(kmeans_k$cluster, dist(tsne_type$Y))
  sil_tsne_type[k-1] <- mean(ss[, 3])
}
optimal_k_tsne_type <- which.max(sil_tsne_type) + 1

set.seed(123)
sil_umap_type <- c()
for(k in 2:4){
  kmeans_k <- kmeans(umap_type$layout, centers = k, nstart = 50)
  ss <- silhouette(kmeans_k$cluster, dist(umap_type$layout))
  sil_umap_type[k-1] <- mean(ss[, 3])
}
optimal_k_umap_type <- which.max(sil_umap_type) + 1

# Klasteryzacja dla t-SNE
set.seed(123)
kmeans_tsne_type <- kmeans(tsne_type$Y, centers = optimal_k_tsne_type, nstart = 50)
tsne_df_type$Cluster <- as.factor(kmeans_tsne_type$cluster)

# Klasteryzacja dla UMAP
set.seed(123)
kmeans_umap_type <- kmeans(umap_type$layout, centers = optimal_k_umap_type, nstart = 50)
umap_df_type$Cluster <- as.factor(kmeans_umap_type$cluster)

# Wizualizacja z klasterami (Etap 2)
tsne_cluster_plot_type <- ggplot(tsne_df_type, aes(x = Dim1, y = Dim2, color = Cluster)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_mark_hull(aes(fill = Cluster), alpha = 0.2, show.legend = FALSE) +
  labs(title = "t-SNE kmeans Clustering", x = "t-SNE Dimension 1", y = "t-SNE Dimension 2") +
  theme_minimal()

umap_cluster_plot_type <- ggplot(umap_df_type, aes(x = Dim1, y = Dim2, color = Cluster)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_mark_hull(aes(fill = Cluster), alpha = 0.2, show.legend = FALSE) +
  labs(title = "UMAP kmeans Clustering", x = "UMAP Dimension 1", y = "UMAP Dimension 2") +
  theme_minimal()

gridExtra::grid.arrange(
  autoplot_tsne_group_type,
  autoplot_tsne_amino_type,
  tsne_cluster_plot_type,
  autoplot_umap_group_type,
  autoplot_umap_amino_type,
  umap_cluster_plot_type,
  ncol = 3, nrow = 2
)
}

# --- Ekstrakcja pełnych typów tRNA dla każdego klastra ---
# Dla t-SNE
tsne_clusters_full <- tsne_df_type %>%
  group_by(Cluster, Amino, Group, Type) %>% 
  summarise(.groups = "drop") %>%  # usunięcie duplikatów
  arrange(Cluster, Group, Amino)

# Dla UMAP
umap_clusters_full <- umap_df_type %>%
  group_by(Cluster, Amino, Group, Type) %>% 
  summarise(.groups = "drop") %>% 
  arrange(Cluster, Group, Amino)

# --- Tworzenie interaktywnych tabel DT ---
# Tabela dla t-SNE
tsne_DT <- datatable(
  tsne_clusters_full,
  colnames = c("Cluster", "Amino Acid", "Group", "Type"),
  caption = "t-SNE Clusters with tRNA Types",
  filter = "top",
  options = list(
    pageLength = 15,
    autoWidth = TRUE,
    scrollX = TRUE
  )
)

# Tabela dla UMAP
umap_DT <- datatable(
  umap_clusters_full,
  colnames = c("Cluster", "Amino Acid", "Group", "Type"),
  caption = "UMAP Clusters with tRNA Types",
  filter = "top",
  options = list(
    pageLength = 15,
    autoWidth = TRUE,
    scrollX = TRUE
  )
)


```

<details>

<summary>Tables with tRNA types per cluster and per genome type</summary>

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 16}
tsne_DT
umap_DT


tsne_clusters_full <- tsne_df_type %>%
  group_by(Cluster) %>%
  summarise(
    Amino = toString(unique(Amino)),  
    Group_Types = paste0(
      sapply(unique(Group), function(g) {
        types <- unique(Type[Group == g])
        paste0(g, ": (", toString(types), ")")
      }), collapse = ", "
    )   
  )

# Dla UMAP
umap_clusters_full <- umap_df_type %>%
  group_by(Cluster) %>%
  summarise(
    Amino = toString(unique(Amino)),  
    Group_Types = paste0(
      sapply(unique(Group), function(g) {
        types <- unique(Type[Group == g])
        paste0(g, ": (", toString(types), ")")
      }), collapse = ", "
    )
  )

tsne_table <- kable(
  tsne_clusters_full,
  format = "html",  
  col.names = c("Cluster", "Amino Acids", "Grouped Types (AA_Codon)"),
  caption = "t-SNE Clusters with tRNA Types"
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) %>%
  scroll_box(width = "100%", height = "400px")

umap_table <- kable(
  umap_clusters_full,
  format = "html",  
  col.names = c("Cluster", "Amino Acids", "Grouped Types (AA_Codon)"),
  caption = "UMAP Clusters with tRNA Types"
) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) %>%
  scroll_box(width = "100%", height = "400px")

tsne_table
```

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 16}
umap_table
```

</details>

</details>

<details>

<summary>UHGG database</summary>

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 20}
Type_Per_Acc <- read_table("~/Praktyki_Bioinfa/tRNA_by_type/Type_Per_Acc_UHGG.txt", 
                             col_names = c("Count", "Acc"))

Type_Per_Acc$Type <- Type_Per_Acc %>%
  pull(Acc) %>% 
  str_split("_", simplify = TRUE) %>%
  .[, 3:4] %>% 
  apply(1, function(x) paste(x, collapse = "_"))

Type_Per_Acc$Acc <- Type_Per_Acc %>%
  pull(Acc) %>% 
  str_split("_", simplify = TRUE) %>%
  .[, 1:2] %>% 
  apply(1, function(x) paste(x, collapse = "_"))

Type_Per_Acc <- Type_Per_Acc %>% 
  filter(Type != "unknown_NNN")

df <- Type_Per_Acc %>% 
  mutate(Amino = sub("_.*", "", Type))
{
# Etap 1: Analiza po osobnikach (genomach)
# Tworzymy macierz Acc × Type (counts)
mat <- tidyr::pivot_wider(df, 
                          id_cols = Acc,
                          names_from = Type, 
                          values_from = Count, 
                          values_fill = list(Count = 0))

# Obliczenie całkowitej liczby tRNA (Total tRNA Count) na genom
tRNA_total <- df %>%
  group_by(Acc) %>%
  summarise(Total_tRNA = sum(Count), .groups = "drop")

# Obliczenie najczęstszego kodonu (Most Frequent Codon) na genom
codon_info <- df %>%
  group_by(Acc, Type) %>%
  summarise(total_count = sum(Count), .groups = "drop") %>%
  group_by(Acc) %>%
  slice_max(total_count, n = 1, with_ties = FALSE) %>%
  ungroup() %>%
  select(Acc, Type) %>%
  dplyr::rename(Most_Frequent_Codon = Type)
# Normalizacja CPM (counts per million)
mat_numeric <- mat %>%
  select(-Acc) %>%
  mutate(across(everything(), as.numeric)) %>%
  mutate(across(everything(), ~ . / sum(.) * 1e6))  # CPM normalization

# Macierz z rownames
mat_matrix <- as.matrix(mat_numeric)
rownames(mat_matrix) <- mat$Acc

# --- PCA jako preprocessing ---
pca_result <- prcomp(mat_matrix, scale. = TRUE)
pca_data <- pca_result$x[, 1:50]  # Użycie 50 głównych składowych

# --- Przygotowanie danych bez duplikatów ---
duplicate_rows <- duplicated(pca_data)
mat_nodup <- pca_data[!duplicate_rows, , drop = FALSE]
acc_nodup <- rownames(mat_matrix)[!duplicate_rows]

# Filtrowanie metadanych dla nodup Acc
tRNA_total_nodup <- tRNA_total %>% filter(Acc %in% acc_nodup)
codon_info_nodup <- codon_info %>% filter(Acc %in% acc_nodup)

# --- t-SNE (parametry dla 14800 genomów) ---
set.seed(123)
tsne <- Rtsne(mat_nodup, perplexity = 50, theta = 0.5, max_iter = 2000, check_duplicates = FALSE)

# Data frame dla t-SNE
tsne_df <- data.frame(
  Dim1 = tsne$Y[, 1],
  Dim2 = tsne$Y[, 2],
  Acc = acc_nodup,
  Total_tRNA = tRNA_total_nodup$Total_tRNA,
  Most_Frequent_Codon = codon_info_nodup$Most_Frequent_Codon
)

# --- UMAP (parametry dla 14800 genomów) ---
set.seed(123)
umap_result <- umap(mat_nodup, n_neighbors = 30, min_dist = 0.1, metric = "euclidean", n_epochs = 500)

# Data frame dla UMAP
umap_df <- data.frame(
  Dim1 = umap_result$layout[, 1],
  Dim2 = umap_result$layout[, 2],
  Acc = acc_nodup,
  Total_tRNA = tRNA_total_nodup$Total_tRNA,
  Most_Frequent_Codon = codon_info_nodup$Most_Frequent_Codon
)

# --- Wizualizacja po Total tRNA Count (Etap 1) ---
autoplot_tsne_total <- ggplot(tsne_df, aes(x = Dim1, y = Dim2, color = Total_tRNA)) +
  geom_point(size = 2, alpha = 0.6) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "t-SNE - Total tRNA Count (Individual Genomes)",
       x = "t-SNE Dimension 1", y = "t-SNE Dimension 2") +
  theme_minimal() +
  theme(legend.position = "right")

autoplot_umap_total <- ggplot(umap_df, aes(x = Dim1, y = Dim2, color = Total_tRNA)) +
  geom_point(size = 2, alpha = 0.6) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "UMAP - Total tRNA Count (Individual Genomes)",
       x = "UMAP Dimension 1", y = "UMAP Dimension 2") +
  theme_minimal() +
  theme(legend.position = "right")

# --- Wizualizacja po Most Frequent Codon (Amino-Kodon) (Etap 1) ---
autoplot_tsne_amino_kodon <- ggplot(tsne_df, aes(x = Dim1, y = Dim2, color = Most_Frequent_Codon)) +
  geom_point(size = 2, alpha = 0.6) +
  labs(title = "t-SNE - Most Frequent Codon (Individual Genomes)",
       x = "t-SNE Dimension 1", y = "t-SNE Dimension 2") +
  theme_minimal() +
  theme(legend.position = "right",
        legend.text = element_text(size = 5),
        legend.key.size = unit(0.6, "lines"))

autoplot_umap_amino_kodon <- ggplot(umap_df, aes(x = Dim1, y = Dim2, color = Most_Frequent_Codon)) +
  geom_point(size = 2, alpha = 0.6) +
  labs(title = "UMAP - Most Frequent Codon (Individual Genomes)",
       x = "UMAP Dimension 1", y = "UMAP Dimension 2") +
  theme_minimal() +
  theme(legend.position = "right",
        legend.text = element_text(size = 5),
        legend.key.size = unit(0.6, "lines"))

# --- Klasteryzacja (Etap 1) - bardziej konserwatywna (zakres k 2-5) ---
set.seed(123)
sil_tsne <- c()
for(k in 2:5){
  kmeans_k <- kmeans(tsne$Y, centers = k, nstart = 50)
  ss <- silhouette(kmeans_k$cluster, dist(tsne$Y))
  sil_tsne[k-1] <- mean(ss[, 3])
}
optimal_k_tsne <- which.max(sil_tsne) + 1

set.seed(123)
sil_umap <- c()
for(k in 2:5){
  kmeans_k <- kmeans(umap_result$layout, centers = k, nstart = 50)
  ss <- silhouette(kmeans_k$cluster, dist(umap_result$layout))
  sil_umap[k-1] <- mean(ss[, 3])
}
optimal_k_umap <- which.max(sil_umap) + 1

# Klasteryzacja dla t-SNE
set.seed(123)
kmeans_tsne <- kmeans(tsne$Y, centers = optimal_k_tsne, nstart = 50)
tsne_df$Cluster <- as.factor(kmeans_tsne$cluster)

# Klasteryzacja dla UMAP
set.seed(123)
kmeans_umap <- kmeans(umap_result$layout, centers = optimal_k_umap, nstart = 50)
umap_df$Cluster <- as.factor(kmeans_umap$cluster)

# Wizualizacja z klasterami (Etap 1)
tsne_cluster_plot <- ggplot(tsne_df, aes(x = Dim1, y = Dim2, color = Cluster)) +
  geom_point(size = 2, alpha = 0.6) +
  geom_mark_hull(aes(fill = Cluster), alpha = 0.2, show.legend = FALSE) +
  labs(title = "t-SNE Clustering (Individual groups)",
       x = "t-SNE Dimension 1", y = "t-SNE Dimension 2") +
  theme_minimal()

umap_cluster_plot <- ggplot(umap_df, aes(x = Dim1, y = Dim2, color = Cluster)) +
  geom_point(size = 2, alpha = 0.6) +
  geom_mark_hull(aes(fill = Cluster), alpha = 0.2, show.legend = FALSE) +
  labs(title = "UMAP Clustering (Individual groups)",
       x = "UMAP Dimension 1", y = "UMAP Dimension 2") +
  theme_minimal()

# Potrójne wykresy dla Etap 1
gridExtra::grid.arrange(
  autoplot_tsne_total,
  autoplot_tsne_amino_kodon,
  tsne_cluster_plot,
  autoplot_umap_total,
  autoplot_umap_amino_kodon,
  umap_cluster_plot,
  ncol = 3, nrow = 2
)
}

umap_clusters_acc <- umap_df %>%
  group_by(Cluster) %>%
  summarise(Accession = unique(Acc), .groups = "drop")

```

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 20}
{# Wczytanie danych
# Wczytanie danych
Type_Per_Acc <- read_table("~/Praktyki_Bioinfa/tRNA_by_type/Type_Per_Acc_UHGG.txt", 
                           col_names = c("Count", "Acc"))

# Przetwarzanie kolumn Acc i Type (zachowujemy pełny Type jako Amino_Codon)
Type_Per_Acc$Type <- Type_Per_Acc %>%
  pull(Acc) %>% 
  str_split("_", simplify = TRUE) %>%
  .[, 3:4] %>% 
  apply(1, function(x) paste(x, collapse = "_"))  # Pełny Type (np. Met_CAT)

Type_Per_Acc$Acc <- Type_Per_Acc %>%
  pull(Acc) %>% 
  str_split("_", simplify = TRUE) %>%
  .[, 1:2] %>% 
  apply(1, function(x) paste(x, collapse = "_"))

# Filtrowanie unknown
Type_Per_Acc <- Type_Per_Acc %>% 
  filter(Type != "unknown_NNN")

# Wydobycie pełnego Type jako identyfikator (Amino_Codon) i Aminokwasu
df <- Type_Per_Acc %>% 
  mutate(Type_Full = Type,  # Zachowujemy pełny Type
         Amino = sub("_.*", "", Type))  # Wyodrębnienie Aminokwasu (np. Met z Met_CAT)

{
# Etap 2: Analiza po typach tRNA
# Tworzymy macierz Type × Acc (counts)
mat_type <- tidyr::pivot_wider(df, 
                               id_cols = Type, 
                               names_from = Acc, 
                               values_from = Count, 
                               values_fill = list(Count = 0))

# Zachowujemy informację o Amino dla kolorowania
amino_info_type <- df %>% 
  select(Type, Amino) %>% 
  distinct()

# Normalizacja CPM (counts per million)
mat_numeric_type <- mat_type %>%
  select(-Type) %>%
  mutate(across(everything(), as.numeric)) %>%
  mutate(across(everything(), ~ . / sum(.) * 1e6))  # CPM normalization

# Macierz z rownames
mat_matrix_type <- as.matrix(mat_numeric_type)
rownames(mat_matrix_type) <- mat_type$Type

# --- PCA jako preprocessing ---
pca_result_type <- prcomp(mat_matrix_type, scale. = TRUE)
pca_data_type <- pca_result_type$x[, 1:20]  # Użycie 20 głównych składowych

# --- Przygotowanie danych bez duplikatów ---
duplicate_rows_type <- duplicated(pca_data_type)
mat_nodup_type <- pca_data_type[!duplicate_rows_type, , drop = FALSE]
type_nodup <- rownames(mat_matrix_type)[!duplicate_rows_type]

# Filtrowanie metadanych dla nodup Type
amino_info_nodup_type <- amino_info_type %>% 
  filter(Type %in% type_nodup)

# --- t-SNE dla typów ---
set.seed(123)
tsne_type <- Rtsne(mat_nodup_type, perplexity = 20, theta = 0.5, max_iter = 1000, check_duplicates = FALSE)

# Data frame dla t-SNE
tsne_df_type <- data.frame(
  Dim1 = tsne_type$Y[, 1],
  Dim2 = tsne_type$Y[, 2],
  Type = type_nodup
)

# Dołączanie Amino
tsne_df_type <- tsne_df_type %>% 
  left_join(amino_info_nodup_type, by = "Type")

# --- UMAP dla typów ---
set.seed(123)
umap_type <- umap(mat_nodup_type, n_neighbors = 15, min_dist = 0.05, metric = "cosine", n_epochs = 200)

# Data frame dla UMAP
umap_df_type <- data.frame(
  Dim1 = umap_type$layout[, 1],
  Dim2 = umap_type$layout[, 2],
  Type = type_nodup
)

# Dołączanie Amino
umap_df_type <- umap_df_type %>% 
  left_join(amino_info_nodup_type, by = "Type")

# --- Wizualizacja po Total tRNA Count (Etap 2) ---
# Obliczenie Total tRNA na typ
tRNA_total_type <- df %>%
  group_by(Type) %>%
  summarise(Total_tRNA = sum(Count), .groups = "drop")

tRNA_total_type_nodup <- tRNA_total_type %>% filter(Type %in% type_nodup)

tsne_df_type$Total_tRNA <- tRNA_total_type_nodup$Total_tRNA
umap_df_type$Total_tRNA <- tRNA_total_type_nodup$Total_tRNA

autoplot_tsne_total_type <- ggplot(tsne_df_type, aes(x = Dim1, y = Dim2, color = Total_tRNA)) +
  geom_point(size = 3, alpha = 0.7) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "t-SNE - Total tRNA Count (AA-codon types)",
       x = "t-SNE Dimension 1", y = "t-SNE Dimension 2") +
  theme_minimal() +
  theme(legend.position = "right")

autoplot_umap_total_type <- ggplot(umap_df_type, aes(x = Dim1, y = Dim2, color = Total_tRNA)) +
  geom_point(size = 3, alpha = 0.7) +
  scale_color_gradient(low = "blue", high = "red") +
  labs(title = "UMAP - Total tRNA Count (AA-codon types)",
       x = "UMAP Dimension 1", y = "UMAP Dimension 2") +
  theme_minimal() +
  theme(legend.position = "right")

# --- Wizualizacja po Amino-Kodon (Etap 2) ---
autoplot_tsne_amino_type <- ggplot(tsne_df_type, aes(x = Dim1, y = Dim2, color = Amino)) +
  geom_point(size = 3, alpha = 0.7) +
  labs(title = "t-SNE - Amino-codon (AA-codon types)",
       x = "t-SNE Dimension 1", y = "t-SNE Dimension 2") +
  theme_minimal() +
  theme(legend.position = "right",  # Przeniesienie legendy na górę
        legend.text = element_text(size = 4),  # Mniejszy tekst
        legend.key.size = unit(0.7, "lines"))

autoplot_umap_amino_type <- ggplot(umap_df_type, aes(x = Dim1, y = Dim2, color = Amino)) +
  geom_point(size = 3, alpha = 0.7) +
  labs(title = "UMAP - Amino-codon (AA-codon types)",
       x = "UMAP Dimension 1", y = "UMAP Dimension 2") +
  theme_minimal() +
  theme(legend.position = "right",
        legend.text = element_text(size = 4),
        legend.key.size = unit(0.7, "lines"))

# --- Klasteryzacja (Etap 2) - bardziej konserwatywna (zakres k 2-4) ---
set.seed(123)
sil_tsne_type <- c()
for(k in 2:4){
  kmeans_k <- kmeans(tsne_type$Y, centers = k, nstart = 50)
  ss <- silhouette(kmeans_k$cluster, dist(tsne_type$Y))
  sil_tsne_type[k-1] <- mean(ss[, 3])
}
optimal_k_tsne_type <- which.max(sil_tsne_type) + 1

set.seed(123)
sil_umap_type <- c()
for(k in 2:4){
  kmeans_k <- kmeans(umap_type$layout, centers = k, nstart = 50)
  ss <- silhouette(kmeans_k$cluster, dist(umap_type$layout))
  sil_umap_type[k-1] <- mean(ss[, 3])
}
optimal_k_umap_type <- which.max(sil_umap_type) + 1

# Klasteryzacja dla t-SNE
set.seed(123)
kmeans_tsne_type <- kmeans(tsne_type$Y, centers = optimal_k_tsne_type, nstart = 50)
tsne_df_type$Cluster <- as.factor(kmeans_tsne_type$cluster)

# Klasteryzacja dla UMAP
set.seed(123)
kmeans_umap_type <- kmeans(umap_type$layout, centers = optimal_k_umap_type, nstart = 50)
umap_df_type$Cluster <- as.factor(kmeans_umap_type$cluster)

# Wizualizacja z klasterami (Etap 2)
tsne_cluster_plot_type <- ggplot(tsne_df_type, aes(x = Dim1, y = Dim2, color = Cluster)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_mark_hull(aes(fill = Cluster), alpha = 0.2, show.legend = FALSE) +
  labs(title = "t-SNE Clustering (AA-codon types grups)",
       x = "t-SNE Dimension 1", y = "t-SNE Dimension 2") +
  theme_minimal()

umap_cluster_plot_type <- ggplot(umap_df_type, aes(x = Dim1, y = Dim2, color = Cluster)) +
  geom_point(size = 3, alpha = 0.7) +
  geom_mark_hull(aes(fill = Cluster), alpha = 0.2, show.legend = FALSE) +
  labs(title = "UMAP Clustering (AA-codon type grups)",
       x = "UMAP Dimension 1", y = "UMAP Dimension 2") +
  theme_minimal()

# Potrójne wykresy dla Etap 2
gridExtra::grid.arrange(
  autoplot_tsne_total_type,
  autoplot_tsne_amino_type,
  tsne_cluster_plot_type,
  autoplot_umap_total_type,
  autoplot_umap_amino_type,
  umap_cluster_plot_type,
  ncol = 3, nrow = 2
)
}

# --- Ekstrakcja pełnych typów tRNA dla każdego klastra ---
# --- Ekstrakcja pełnych typów tRNA dla każdego klastra ---

# Dla t-SNE
tsne_clusters_full <- tsne_df_type %>%
  group_by(Cluster) %>%
  summarise(
    Amino = toString(unique(Amino)),   # Łączenie unikalnych aminokwasów w ciąg
    Type_Full = toString(unique(Type)) # Łączenie unikalnych typów w ciąg
  )

# Dla UMAP
umap_clusters_full <- umap_df_type %>%
  group_by(Cluster) %>%
  summarise(
    Amino = toString(unique(Amino)),   # Łączenie unikalnych aminokwasów w ciąg
    Type_Full = toString(unique(Type)) # Łączenie unikalnych typów w ciąg
  )
}
```

<details>

<summary>UHGG Tables with tRNA types per cluster and per acc</summary>

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 16}
datatable(tsne_clusters_full,
          filter = "top",
          options = list(max(levels(umap_clusters_acc$Cluster)),
                         autoWidth = TRUE))

datatable(
  umap_clusters_acc,
  filter = "top",
  options = list(
    pageLength = 10,
    autoWidth = TRUE,
    scrollX = TRUE
  )
)
```

</details>

</details>

</details>

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

<details>

<summary>UHGG database analysis</summary>

I extended R script and removed tRNAscan-se from pipeline and performed analysis on [UHGG database](https://www.ebi.ac.uk/metagenomics/genome-catalogues/human-gut-v2-0-2) (31.07.2025)

I used lftp to fetch whole database including specifically genomes and pan-genomes only

```{bash eval=FALSE}
lftp -c "open https://ftp.ebi.ac.uk; mirror --parallel=4 --include-glob=*.fna --exclude-glob=*.txt --exclude-glob=*gz --exclude-glob=*.tar.gz --verbose /pub/databases/metagenomics/mgnify_genomes/human-gut/v2.0.2/species_catalogue/ ./species_catalogue"
```

<details>

<summary>Pipeline</summary>

```{bash eval=FALSE}
#!/bin/bash

# Base path (script location)
BASE_DIR=$(dirname "$0")
GENOMES_DIR="$BASE_DIR/genomes"
TEMP_DIR="$BASE_DIR/temp"
RESULTS_DIR="$BASE_DIR/results"
CODETTA_DIR="$BASE_DIR/codetta"

# Create directories if they don't exist
mkdir -p "$TEMP_DIR"
mkdir -p "$RESULTS_DIR"

# Iterate over all .fna files
for genome_path in "$GENOMES_DIR"/*.fna; do
    # Filename without path and extension
    genome_file=$(basename "$genome_path")
    genome_name="${genome_file%.fna}"

    echo "Processing: $genome_name"

    # Temporary files
    aragorn_out="$TEMP_DIR/${genome_name}_aragorn.out"
    codon_summary="$TEMP_DIR/${genome_name}_codons_summary.txt"
    codetta_file="$TEMP_DIR/${genome_name}_codetta_results.txt"

    # Run ARAGORN
    aragorn -gcstd -w -t -o "$aragorn_out" "$genome_path"

    # Run Codetta
    python3.10 "$CODETTA_DIR/codetta.py" -e 1e-50 --results_summary "$codetta_file" "$genome_path"

    # Extract codons from ARAGORN
    if [[ -f "$aragorn_out" ]]; then
        grep -v "^#" "$aragorn_out" | \
        awk '{ gsub(/\(|\)/, "", $5); codon=toupper($5); if (codon ~ /^[ATGC]{3}$/) print codon }' \
        > "$TEMP_DIR/${genome_name}_codons_aragorn.txt"
    else
        echo "Missing file $aragorn_out"
    fi

    # Count unique codons
    cat "$TEMP_DIR/${genome_name}_codons_aragorn.txt" \
        | sort | uniq -c | awk '{print $2","$1}' > "$codon_summary"

    # Run R script for visualization and analysis
    histogram_file="$RESULTS_DIR/${genome_name}_codon_histogram.png"
    unused_codons="$RESULTS_DIR/${genome_name}_unused_codons.txt"

    if [[ -f "$codon_summary" && -f "$codetta_file" ]]; then
        Rscript codon_plot.R "$codon_summary" "$unused_codons" "$histogram_file" "$genome_name" "$codetta_file"
    else
        echo "Missing file $codon_summary or $codetta_file for visualization."
    fi

    # Remove all temporary files for this genome
    rm -f "$TEMP_DIR/${genome_name}"_*
    find "$GENOMES_DIR" -type f -name "${genome_name}*" ! -name "${genome_name}.fna" ! -name "${genome_name}.fasta" -delete


    echo "Saved results for $genome_name in results/ directory"
done

```

</details>

<details>

<summary>R script</summary>

```{r eval=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 16}
library(ggplot2)
library(glue)
# Define standard genetic code
basic_code <- c(
  "UUU" = "F", "UUC" = "F", "UUA" = "L", "UUG" = "L",
  "UCU" = "S", "UCC" = "S", "UCA" = "S", "UCG" = "S",
  "UAU" = "Y", "UAC" = "Y", "UAA" = "*", "UAG" = "*",
  "UGU" = "C", "UGC" = "C", "UGA" = "*", "UGG" = "W",
  "CUU" = "L", "CUC" = "L", "CUA" = "L", "CUG" = "L",
  "CCU" = "P", "CCC" = "P", "CCA" = "P", "CCG" = "P",
  "CAU" = "H", "CAC" = "H", "CAA" = "Q", "CAG" = "Q",
  "CGU" = "R", "CGC" = "R", "CGA" = "R", "CGG" = "R",
  "AUU" = "I", "AUC" = "I", "AUA" = "I", "AUG" = "M",
  "ACU" = "T", "ACC" = "T", "ACA" = "T", "ACG" = "T",
  "AAU" = "N", "AAC" = "N", "AAA" = "K", "AAG" = "K",
  "AGU" = "S", "AGC" = "S", "AGA" = "R", "AGG" = "R",
  "GUU" = "V", "GUC" = "V", "GUA" = "V", "GUG" = "V",
  "GCU" = "A", "GCC" = "A", "GCA" = "A", "GCG" = "A",
  "GAU" = "D", "GAC" = "D", "GAA" = "E", "GAG" = "E",
  "GGU" = "G", "GGC" = "G", "GGA" = "G", "GGG" = "G"
)

# Function to create genome-specific genetic code from Codetta CSV output
create_genome_specific_code <- function(codetta_file, base_code) {
  # Read Codetta output as CSV
  codetta_data <- read.csv(codetta_file, header = TRUE, stringsAsFactors = FALSE)

  # Extract inferred_gencode (assuming one row per file)
  genetic_code_str <- codetta_data$inferred_gencode[1]

  # Define codon order as per Codetta (UUU, UUC, ..., GGG)
  codons <- c("UUU", "UUC", "UUA", "UUG", "UCU", "UCC", "UCA", "UCG",
              "UAU", "UAC", "UAA", "UAG", "UGU", "UGC", "UGA", "UGG",
              "CUU", "CUC", "CUA", "CUG", "CCU", "CCC", "CCA", "CCG",
              "CAU", "CAC", "CAA", "CAG", "CGU", "CGC", "CGA", "CGG",
              "AUU", "AUC", "AUA", "AUG", "ACU", "ACC", "ACA", "ACG",
              "AAU", "AAC", "AAA", "AAG", "AGU", "AGC", "AGA", "AGG",
              "GUU", "GUC", "GUA", "GUG", "GCU", "GCC", "GCA", "GCG",
              "GAU", "GAC", "GAA", "GAG", "GGU", "GGC", "GGA", "GGG")

  # Split genetic code string into individual amino acids
  inferred_aa <- strsplit(genetic_code_str, "")[[1]]

  # Create new genetic code matrix
  genome_code <- base_code
  for (i in 1:length(codons)) {
    inferred <- inferred_aa[i]
    standard <- base_code[codons[i]]
    # Update only if inferred is not "?" and different from standard
    if (inferred != "?" && inferred != standard) {
      genome_code[codons[i]] <- inferred
    }
  }

  return(genome_code)
}

# Function to map tRNA anticodons to mRNA codons
map_tRNA_to_mRNA <- function(antykodony_DNA, species = "Organism", genetic_code) {
  tRNA_RNA <- toupper(chartr("T", "U", antykodony_DNA))
  tRNA_RNA <- tRNA_RNA[nchar(tRNA_RNA) == 3 & grepl("^[ACGU]+$", tRNA_RNA)]

  wobble_map_small <- list(
    "G" = c("U","C"),
    "C" = c("G"),
    "A" = c("U"),
    "U" = c("A","G")
  )

  wobble_map_big <- list(
    "G" = c("U","C"),
    "C" = c("G","A"),
    "A" = c("U","C","A"),
    "U" = c("A","G","U","C"),
    "I" = c("A","U","C"),
    "*" = c("A","U","C", "G")
  )

  bazy <- c("U", "C", "A", "G")
  kodony_mRNA <- as.vector(outer(outer(bazy, bazy, paste0), bazy, paste0))

  rna_comp <- function(x) chartr("ACGU", "UGCA", x)

  calculate_recognizable <- function(wobble_map) {
    rozpoznawalne <- unique(unlist(lapply(tRNA_RNA, function(anty) {
      a1 <- substr(anty, 1, 1)  # 5'
      a2 <- substr(anty, 2, 2)
      a3 <- substr(anty, 3, 3)  # 3'
      Filter(function(kodon) {
        k1 <- substr(kodon, 1, 1)  # 5'
        k2 <- substr(kodon, 2, 2)
        k3 <- substr(kodon, 3, 3)  # 3'
        k3 %in% wobble_map[[a1]] &&
          k2 == rna_comp(a2) &&
          k1 == rna_comp(a3)
      }, kodony_mRNA)
    })))

    if(genetic_code[15] == "*") {
      rozpoznawalne <- rozpoznawalne[!rozpoznawalne %in% c("UGA")]
    }
    if(genetic_code[11] == "*") {
      rozpoznawalne <- rozpoznawalne[!rozpoznawalne %in% c("UAA")]
    }
    if(genetic_code[12] == "*") {
      rozpoznawalne <- rozpoznawalne[!rozpoznawalne %in% c("UAG")]
    }

    return(rozpoznawalne)
  }

  rozpoznawalne_big <- calculate_recognizable(wobble_map_big)
  rozpoznawalne_small <- calculate_recognizable(wobble_map_small)

  df <- data.frame(
    codon = kodony_mRNA,
    b1 = substr(kodony_mRNA, 1, 1),
    b2 = substr(kodony_mRNA, 2, 2),
    b3 = substr(kodony_mRNA, 3, 3),
    aa = genetic_code[kodony_mRNA],
    status_big = ifelse(kodony_mRNA %in% rozpoznawalne_big, "Found", "Missing"),
    status_small = ifelse(kodony_mRNA %in% rozpoznawalne_small, "Found", "Missing")
  )

  df$b1 <- factor(df$b1, levels = rev(bazy))
  df$b2 <- factor(df$b2, levels = bazy)
  df$b3 <- factor(df$b3, levels = bazy)
  df$label <- paste0(df$codon, "\n", df$aa)

  p <- ggplot(df, aes(x = b3, y = b1, fill = status_big)) +
    geom_tile(color = "black", width = 0.95, height = 0.95) +
    geom_text(aes(label = label), size = 3.2, fontface = "bold", lineheight = 0.8) +
    facet_wrap(~b2, nrow = 2, labeller = label_both) +
    scale_fill_manual(values = c("Found" = "lightblue", "Missing" = "pink")) +
    theme_minimal(base_size = 14) +
    theme(
      panel.grid = element_blank(),
      strip.text = element_text(size = 14),
      axis.title = element_text(size = 14),
      plot.title = element_text(hjust = 0.5)
    ) +
    labs(
      title = glue("Map of uncodable amino acids in {species}"),
      x = "Second nucleotide", y = "First nucleotide", fill = "Status"
    )

  return(list(
    wykres = p,
    rozpoznawalne_kodony_big = sort(rozpoznawalne_big),
    brakujace_kodony_big = sort(setdiff(kodony_mRNA, rozpoznawalne_big)),
    rozpoznawalne_kodony_small = sort(rozpoznawalne_small),
    brakujace_kodony_small = sort(setdiff(kodony_mRNA, rozpoznawalne_small)),
    tRNA_RNA = tRNA_RNA,
    df_results = df
  ))
}

# Read command line arguments
args <- commandArgs(trailingOnly = TRUE)
codon_summary_file <- args[1]
unused_codons_file <- args[2]
histogram_file <- args[3]
species_name <- args[4]
codetta_file <- args[5]

# Read codon summary data
codon_data <- read.csv(codon_summary_file, header = FALSE, col.names = c("Codon", "Count"))

# Map DNA codons to RNA
codon_data$Codon <- toupper(chartr("T", "U", codon_data$Codon))

# Create genome-specific genetic code from Codetta output
genome_specific_code <- create_genome_specific_code(codetta_file, basic_code)

# Use genome-specific genetic code in analysis
result <- map_tRNA_to_mRNA(codon_data$Codon, species = species_name, genetic_code = genome_specific_code)

# Save plot
ggsave(histogram_file, plot = result$wykres, width = 10, height = 8, dpi = 300)

# Save unused codons
writeLines(
  c(
    "=== Wobble Big ===",
    result$brakujace_kodony_big,
    "",
    "=== Wobble Small ===",
    result$brakujace_kodony_small,
    "",
    "Genetic code:",
    paste(genome_specific_code, collapse = "")
  ),
  unused_codons_file
)
```

</details>

Currently pipeline is still running and results will be added later. (02.08.2025)

</details>

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

<details>

<summary>tRNA to Proteome relation</summary>

Main goal of this section is to explore possibility of protein expression range prediction based of sole genome using tRNA genes as main markers.

<details>

<summary>tRNA genes and amino acid counts</summary>

I made a semi-automated script for genome fetching, tRNA gene extraction and proteome amino acid counting. Using this data I made R statistic calculations and plotted results.

<details>

<summary>Bash</summary>

```{bash eval=FALSE}
#!/bin/bash

# Argumenty
ACCESSION=$1
SPECIES=$2

# Pobierz dane z NCBI
datasets download genome accession "$ACCESSION" --include protein,genome --filename "dataset_${SPECIES}.zip"

# Rozpakuj
unzip -o "dataset_${SPECIES}.zip" -d "ncbi_dataset_${SPECIES}"

# Plik proteomowy (faa)
FAA_FILE=$(find "ncbi_dataset_${SPECIES}" -name "*.faa" | head -n 1)

# Oblicz liczność aminokwasów i zapisz
grep -v "^>" "$FAA_FILE" \
  | fold -w1 \
  | grep -E '^[A-Z]$' \
  | sort \
  | uniq -c \
  | sort -nr > "aa_${SPECIES}.txt"

# Plik genomowy (fna)
FNA_FILE=$(find "ncbi_dataset_${SPECIES}" -name "*genomic.fna" | head -n 1)

# Uruchom aragorn
aragorn -w -t "$FNA_FILE" > "aragorn_${SPECIES}_pre.txt"
gawk 'BEGIN { skip=0 } /^>/ { skip=1; next } skip { skip=0; next } { print }' aragorn_${SPECIES}_pre.txt > aragorn_${SPECIES}.txt

rm aragorn_${SPECIES}_pre.txt
rm -r "ncbi_dataset_${SPECIES}/"
rm "dataset_${SPECIES}.zip"
```

</details>

<details>

<summary>Single genome analysis</summary>

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 20}
analyze_tRNA_aa_single <- function(species_name) {
  # Ścieżki do plików
  aa_path <- glue("~/Praktyki_Bioinfa/Transcript:tRNA/aa_{species_name}.txt")
  aragorn_path <- glue("~/Praktyki_Bioinfa/Transcript:tRNA/aragorn_{species_name}.txt")
  
  # Wczytaj dane
  aa_Proteom <- read_table(aa_path, col_names = c("Count", "AA"))
  aragorn_tRNA <- read_table(aragorn_path, col_names = FALSE)
  
  # Policz liczbę tRNA dla każdego aminokwasu
  tRNA_count <- aragorn_tRNA %>% 
    group_by(X2 = X2) %>%
    summarise(Count = n(), .groups = "drop")
  
  # Usuń SeC z obu data frame
  aa_Proteom <- aa_Proteom %>% filter(AA != "U")
  tRNA_count <- tRNA_count %>% filter(X2 != "tRNA-SeC")
  
  # Oblicz procenty
  aa_Proteom <- aa_Proteom %>%
    mutate(AAperc = Count / sum(Count) * 100)
  
  tRNA_count <- tRNA_count %>%
    mutate(perc = Count / sum(Count) * 100)
  
  # Mapowanie do jednoliterowych kodów aminokwasów
  aa_map <- c(
    Ala = "A", Arg = "R", Asn = "N", Asp = "D", Cys = "C",
    Gln = "Q", Glu = "E", Gly = "G", His = "H", Ile = "I",
    Leu = "L", Lys = "K", Met = "M", Phe = "F", Pro = "P",
    SeC = "U", Ser = "S", Thr = "T", Trp = "W", Tyr = "Y",
    Val = "V"
  )
  
  # Przypisz kody aminokwasów
  tRNA_count$AA <- aa_map[sub("tRNA-", "", tRNA_count$X2)]
  
  # Połącz dane według AA
  full_joined <- full_join(aa_Proteom, tRNA_count, by = "AA", suffix = c("_Proteom", "_tRNA"))
  
  # BOXPLOT base R
  boxplot(full_joined$perc, full_joined$AAperc,
          names = c("tRNA count %", "Proteome %"),
          main = glue("Amino Acid Usage vs tRNA Abundance: {species_name}"),
          ylab = "Percentage",
          col = c("skyblue", "tomato"))
  
  # Testy statystyczne
  diff_vec <- full_joined$perc - full_joined$AAperc
  shap_test <- shapiro.test(diff_vec)
  ttest <- t.test(full_joined$perc, full_joined$AAperc, paired = TRUE)
  
  # Oblicz różnicę i uporządkuj AA wg różnicy (tRNA - Proteome)
  full_joined <- full_joined %>%
    mutate(diff = perc - AAperc) %>%
    arrange(desc(diff))
  
  # Ustaw AA jako faktor w tej kolejności
  full_joined$AA <- factor(full_joined$AA, levels = full_joined$AA)

  long_df <- full_joined %>%
    select(AA, `Proteome` = AAperc, `tRNA` = perc) %>%
    pivot_longer(cols = c("Proteome", "tRNA"), names_to = "Source", values_to = "Percent")
  
  # Wykres ggplot
  barplot <- ggplot(long_df, aes(x = AA, y = Percent, fill = Source)) +
    geom_col(position = position_dodge(width = 0.7), width = 0.6) +
    scale_fill_manual(values = c("Proteome" = "tomato", "tRNA" = "skyblue")) +
    labs(title = glue("Amino Acid Usage vs tRNA Abundance in {species_name}"),
         x = "Amino Acid", y = "Percentage") +
    theme_minimal(base_size = 14) +
    theme(legend.title = element_blank())
  
  print(barplot)
  # Zwróć jako listę
  return(list(
    data = full_joined,
    shapiro = shap_test,
    ttest = ttest,
    barplot = barplot
  ))
}

Ecoli_analysis <- analyze_tRNA_aa_single("Ecoli")
knitr::kable(
  Ecoli_analysis$data,
  caption = "E. coli Amino Acid Usage vs tRNA Abundance",
  format = "html",
  escape = FALSE
)
```

</details>

<details>

<summary>Multiple genomes analysis</summary>

```{r echo=FALSE, message=FALSE, warning=FALSE,fig.show="hide"}
analyze_tRNA_aa <- function(species_name) {
  # Ścieżki do plików
  aa_path <- glue("~/Praktyki_Bioinfa/Transcript:tRNA/auto/aa_{species_name}.txt")
  aragorn_path <- glue("~/Praktyki_Bioinfa/Transcript:tRNA/auto/aragorn_{species_name}.txt")
  
  # Wczytaj dane
  aa_Proteom <- read_table(aa_path, col_names = c("Count", "AA"))
  aragorn_tRNA <- read_table(aragorn_path, col_names = FALSE)
  
  # Policz liczbę tRNA dla każdego aminokwasu
  tRNA_count <- aragorn_tRNA %>% 
    group_by(X2 = X2) %>%
    summarise(Count = n(), .groups = "drop")
  
  # Usuń SeC z obu data frame
  aa_Proteom <- aa_Proteom %>% filter(AA != "U")
  tRNA_count <- tRNA_count %>% filter(X2 != "tRNA-SeC")
  
  # Oblicz procenty
  aa_Proteom <- aa_Proteom %>%
    mutate(AAperc = Count / sum(Count) * 100)
  
  tRNA_count <- tRNA_count %>%
    mutate(perc = Count / sum(Count) * 100)
  
  # Mapowanie do jednoliterowych kodów aminokwasów
  aa_map <- c(
    Ala = "A", Arg = "R", Asn = "N", Asp = "D", Cys = "C",
    Gln = "Q", Glu = "E", Gly = "G", His = "H", Ile = "I",
    Leu = "L", Lys = "K", Met = "M", Phe = "F", Pro = "P",
    SeC = "U", Ser = "S", Thr = "T", Trp = "W", Tyr = "Y",
    Val = "V"
  )
  
  # Przypisz kody aminokwasów
  tRNA_count$AA <- aa_map[sub("tRNA-", "", tRNA_count$X2)]
  
  # Połącz dane według AA
  full_joined <- full_join(aa_Proteom, tRNA_count, by = "AA", suffix = c("_Proteom", "_tRNA"))
  
  # BOXPLOT base R
  boxplot(full_joined$perc, full_joined$AAperc,
          names = c("tRNA count %", "Proteome %"),
          main = glue("Amino Acid Usage vs tRNA Abundance: {species_name}"),
          ylab = "Percentage",
          col = c("skyblue", "tomato"))
  
  # Testy statystyczne
  diff_vec <- full_joined$perc - full_joined$AAperc
  shap_test <- shapiro.test(diff_vec)
  ttest <- t.test(full_joined$perc, full_joined$AAperc, paired = TRUE)
  
  # Oblicz różnicę i uporządkuj AA wg różnicy (tRNA - Proteome)
  full_joined <- full_joined %>%
    mutate(diff = perc - AAperc) %>%
    arrange(desc(diff))
  
  # Ustaw AA jako faktor w tej kolejności
  full_joined$AA <- factor(full_joined$AA, levels = full_joined$AA)

  long_df <- full_joined %>%
    select(AA, `Proteome` = AAperc, `tRNA` = perc) %>%
    pivot_longer(cols = c("Proteome", "tRNA"), names_to = "Source", values_to = "Percent")
  
  # Wykres ggplot
  barplot <- ggplot(long_df, aes(x = AA, y = Percent, fill = Source)) +
    geom_col(position = position_dodge(width = 0.7), width = 0.6) +
    scale_fill_manual(values = c("Proteome" = "tomato", "tRNA" = "skyblue")) +
    labs(title = glue("Amino Acid Usage vs tRNA Abundance in {species_name}"),
         x = "Amino Acid", y = "Percentage") +
    theme_minimal(base_size = 14) +
    theme(legend.title = element_blank())
  
  print(barplot)
  # Zwróć jako listę
  return(list(
    data = full_joined,
    shapiro = shap_test,
    ttest = ttest,
    barplot = barplot
  ))
}

sp_list <- c("Ecoli","Salmonella","Shigella","Klebsiella","Yersinia","Pseudomonas","Bacillus","Staphylococcus","Listeria","Campylobacter")

combined_data <- lapply(sp_list, analyze_tRNA_aa)
combined_data <- bind_rows(lapply(combined_data, function(x) x$data), .id = "Species")
```

```{r echo=FALSE, message=FALSE, warning=FALSE, fig.height= 16,fig.width= 20}
combined_data <- combined_data %>%
  mutate(diff = perc - AAperc) %>%
  arrange(desc(diff))

combined_long <- combined_data %>%
  select(AA, `Proteome` = AAperc, `tRNA` = perc, Species) %>%
  pivot_longer(cols = c("Proteome", "tRNA"), names_to = "Source", values_to = "Percent")

combined_boxplot <- combined_long %>%
  group_by(AA, Source) %>%
  summarise(Percent = list(Percent), .groups = "drop") %>%
  unnest(Percent)

ggplot(combined_boxplot, aes(x = AA, y = Percent, fill = Source)) +
  geom_boxplot(position = position_dodge(width = 0.75), width = 0.6) +
  scale_fill_manual(values = c("Proteome" = "coral", "tRNA" = "skyblue")) +
  labs(title = "Boxplot of Amino Acid Usage vs tRNA Abundance Across Species",
       x = "Amino Acid", y = "Percentage") +
  theme_minimal(base_size = 14) +
  theme(legend.title = element_blank())

combined_data %>%
  mutate(AA = factor(AA, levels = unique(AA[order(diff)]))) %>%
  ggplot(aes(x = AA, y = diff, fill = Species)) +
  geom_bar(stat = "identity", position = "dodge") +
  labs(title = "Difference in tRNA vs Proteome Amino Acid Usage",
       x = "Amino Acid", y = "Difference (tRNA - Proteome)") +
  theme_minimal(base_size = 14) +
  scale_fill_brewer(palette = "Set3") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

combined_avg <- combined_data %>%
  group_by(AA) %>%
  summarise(
    avg_perc = mean(perc, na.rm = TRUE),
    avg_AAperc = mean(AAperc, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(diff = avg_perc - avg_AAperc) %>% 
  arrange(desc(diff))

combined_avg %>%
  mutate(AA = factor(AA, levels = unique(AA[order(diff)]))) %>%
  ggplot(aes(x = AA, y = diff)) +
  geom_bar(stat = "identity", fill = "lightgreen",alpha = 0.7) +
  labs(title = "Average Difference in tRNA vs Proteome Amino Acid Usage",
       x = "Amino Acid", y = "Average Difference (tRNA - Proteome)") +
  theme_minimal(base_size = 14) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
```

</details>

</details>

</details>

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

<details>

<summary>Case study of data errors</summary>

Errors in data can lead to numerous issues. For example in [Thermus aquaticus](https://www.ncbi.nlm.nih.gov/datasets/genome/?taxon=271) taxon there are 3 annotated genomes. I used them in my pipeline to produce ACPMs in search of tRNA genes. One of them (GCA_049415045.1) is lacking coverage for 4 codons two of which are Tyrosine codons. This would mean that this organism don't use this amino acid at all, but of course it was not the case. My main suspicion is bad quality of Metagenome assembly. Two other genomes from this taxon have full sets of ACPM. This should be considered as a warning that data can be reliable only in great numbers and not single genome studies. (24.07.2025)

</details>

</details>

</details>

</details>
