Task_3_JDK.Rmd

---
title: "Praktinė užduotis Nr. 2"
author: "Jaroslav Rutkovskij, Danielė Stasiūnaitė, Karolis Augustauskas (JDK)"
date: "`r format(Sys.time(), '%d/%m/%y')`"
output:
  pdf_document: default
  html_document:
    df_print: paged
subtitle: Klasterizavimas
---

```{r Figure_size, include=FALSE}
knitr::opts_chunk$set(fig.height = 13, fig.width = 16)
```

```{r Preparations, warning=FALSE, include=FALSE}
library(pacman)
p_load(Signac, dplyr, heatmaply, reshape2, ggplot2, data.table, matrixStats,
       WGCNA, ggpubr, tidyverse, palmerpenguins, Rtsne, methylclockData,
       methylclock, ComplexHeatmap, hrbrthemes, dendextend, RColorBrewer,
       Heatplus, gplots, heatmap3, patchwork, circlize)

enableWGCNAThreads()
set.seed(1973)
```

```{r Data_access, warning=FALSE, include=FALSE}
beta <- readRDS("dat8.Rds")               # Modification value matrix
pData_col <- readRDS("pData.Rds")         # Sample (column) matrix
annot_row <- readRDS("annotation.Rds")    # Position (row) matrix
info <- read.csv("data.csv")
```

## Pagrindinės užduotys

Prieš atliekant duomenų hierarchinį klasterizavimą buvo sukurti nauji
mėginių pavadinimai, pakeičiant originalius pavadinimus, kurie buvo ilgi
ir neinformatyvūs. Pavadinimai buvo pakeisti, vadovaujantis žemiau
aprašytu principu:

 1. Nurodomas donoro numeris. Pvz.: ***n37***.

 2. Po apatinio brūkšnio nurodomas laikas, kada buvo imta biopsija (gali
    būti ***t1*** arba ***t2***, kur ***t1*** nurodo, kad biopsija buvo
    imta pirmaisiais metais, o ***t2*** nurodo, kad biopsija buvo paimta
    praėjus 10 metų po pirmosios biopsijos ėmimo).

 3. Po apatinio brūkšnio nurodoma, iš kurios žarnos (proksimalinės ar
    distalinės) buvo imama biopsija. Žarnų tipai užrašomi sutrumpinus
    pilnus pavadinimus iki keturių raidžių:\
    **proximal** $\Rightarrow$ ***prox***,\
    **distal** $\Rightarrow$ ***dist***.

 4. Po apatinio brūkšnio trijų raidžių kombinacija nurodoma, ar
    tiriamasis vartojo aspiriną:\
    ***nnu*** $\Rightarrow$ **aspirinas nevartojamas (non-user)**,\
    ***lng*** $\Rightarrow$ **aspirinas vartojamas (longterm user)**.

 5. Po apatinio brūkšnio nurodomas moters amžius.

Pilnas naujai sukurto pavadinimo pavyzdys: ***n10_t1_prox_nnu_55***,
kur:

 - **n10** - donoro numeris;
 - **t1** - biopsija imta pirmaisiais tyrimo metais;
 - **prox** - mėginys gautas iš proksimalinės žarnos;
 - **nnu** - tiriamasis nevartojo asprino;
 - **55** - tiriamojo amžius.

Mėginių vardų pavyzdys realiems duomenims:

```{r Filename_generation, echo=FALSE, warning=FALSE}
samples <- read.table("samples.txt", header = TRUE) %>% setDT()
samples <- samples[samples$Basename %in% colnames(beta),]
# Sort samples by beta
samples <- samples[order(match(samples$Basename, colnames(beta))),]
samples$graph_name <- paste0("n", samples$donor, "_", samples$timepoint, "_",
                               substring(samples$clonic_location, 1, 4), "_",
                               paste0(substring(samples$aspirin, 1, 1),
                                      substring(samples$aspirin, 3, 4)), "_",
                               samples$age)

colnames(beta) <- samples$graph_name[match(colnames(beta), samples$Basename)]

head(beta[, 1:2], n = 3)
write.csv(samples, "samples2.csv")
```

\newpage

# 1. Mėginių hierarchinis klasterizavimas

```{r Task_1_1, echo=FALSE, message=FALSE, warning=FALSE, fig.show='hold', out.width='90%'}
d <- dist(1 - cor(beta))
clust <- hclust(d, method = "complete")

# Creating color groups:
clonic_location <- labels2colors((samples[samples$graph_name == colnames(beta),]$clonic_location),
                                 colorSeq = c("royalblue", "royalblue4"))
timepoint <- labels2colors((samples[samples$graph_name == colnames(beta),]$timepoint),
                           colorSeq = c("darkolivegreen1", "darkolivegreen4"))
aspirin <- labels2colors((samples[samples$graph_name == colnames(beta),]$aspirin),
                         colorSeq = c("salmon", "seagreen1"))
has_polyp <- labels2colors((samples[samples$graph_name == colnames(beta),]$polyp),
                           colorSeq = c("honeydew2", "lightcoral"))
groups <- cutree(clust, 3)

# Plotting dendrogram:
plotDendroAndColors(clust, cbind(clusters = labels2colors(groups),
                                 clonic_location, has_polyp, aspirin,
                                 timepoint), main = "Cluster dendrogram",
                    ylab = "Aukštis", hang = -1)

plot(NULL, xaxt = 'n', yaxt = 'n', bty = 'n', ylab = '', xlab = '', xlim = 0:1,
     ylim=0:1)

legend("left", pch = 16, bty = 'n', xpd = TRUE,
       legend = c("Proximal","Distal", "t1", "t2", "Aspirin_user",
                  "Aspirin_nonuser", "Polyp_yes", "Polyp_no" ),
       col = c("royalblue4", "royalblue", "darkolivegreen1", "darkolivegreen4",
               "salmon", "seagreen1", "lightcoral", "honeydew2"))
```

Atliekant šią užduotį išskyrėme 3 klasterius. Iš dendrogramos galime
matyti 2 ryškius klasterius su distalinės žarnos ir proksimalinės žarnos
mėginiais. Išskirtame viduriniame klasteryje matosi mėginiai, kurie
turėjo polipų (polyp_yes), daugiau buvo ilgalaikių aspirino naudotojų
(aspirin_user), taip pat daug mėginių paimtų vėlesniu laikotarpiu (t2 -
paėmus mėginį po 10 metų). Daugiau smulkesnių susiskirstymų neišskyrėme.

\newpage

# 2. Mėginių atvaizdavimas "heatmap" pavidalu

Kadangi **beta** matrica turi ***`r dim(beta)[1]`*** eilutes bei
***`r dim(beta)[2]`*** stulpelius, pavaizduoti visas matricos reikšmes
bei, atsižvelgus į spalvų intensyvumą ir atspalvį, daryti išvadas apie
duomenų grupes - klasterius, naudojantis ***heatmap()*** funkcija yra
negalima, iš eilutes aprašančios matricos buvo atrinktos tik tam tikros
eilutės.

**Eilučių atrinkimas buvo atliktas, remiantis žemiau pažingsniui
aprašyta** **logika:**

1.  Iš pradžių iš eilutes aprašančios matricos buvo surasti genai (jų
    pavadinimai, remiantis eilutes aprašančios matricos
    'UCSC_RefGene_Name' stulpeliu), kurie pasikartodavo tarp dviejų
    chromosomų. Gauti bendrų genų dažniai pavaizduoti lentele:

```{r Task_2_1, echo=FALSE, message=FALSE, warning=FALSE}
# Determining unique 'chr' column values (chromosome names) from the matrix
# that describes positions (rows):
chr_names <- unique(annot_row$chr)     
genes_only <- c()

# Loop which checks whether there are common genes between chromosome pairs:
for (i in 1:length(chr_names)) {
  for (y in 1:length(chr_names)) {
    if (i == y) { next }
    else {
      names1 <- annot_row[annot_row$chr == chr_names[i],
                          'UCSC_RefGene_Name'] %>% unique()
      names2 <- annot_row[annot_row$chr == chr_names[y],
                          'UCSC_RefGene_Name'] %>% unique()
      genes <- intersect(names1, names2)
      if (genes[1] == '' & length(genes) == 1) { next }
      else {
        # print(paste0(i, " ir ", y, ": ", genes))
        genes_only <- c(genes_only, genes)
      }
    }
  }
}
```

```{r Task_2_2, include=FALSE}
# Removing empty elements ('') from the gene list:
genes_only <- genes_only[! genes_only %in% c('')]

# Creating gene expression frequency table:
print(table(genes_only))
```

2.  Turint dažnius buvo suformuotas objektas, suskaičiuojantis konkrečių
    genų dažnius chromosomose.

```{r Task_2_3, echo=FALSE}
gene_list <- c(unique(genes_only))
gene_cnts <- data.frame(chr = chr_names)  # gene_conts --> gene count dataframe

for (g in 1:length(gene_list)) {
  gene <- c()
  for (j in 1:length(chr_names)) {
    annot_chr <- annot_row[annot_row$chr == chr_names[j], ]
    annot_gene <- annot_chr[annot_chr$UCSC_RefGene_Name == gene_list[g], ] %>% 
                              data.matrix()
    #if (length(annot_gene) == 0) {next}
    # print(paste0(gene_list[g], ", ", chr_names[j], " = ", length(annot_gene)))
    gene <- c(gene, length(annot_gene))
  }
  gene_cnts[ , ncol(gene_cnts) + 1] <- gene
  colnames(gene_cnts)[ncol(gene_cnts)] <- gene_list[g]
}
```

3.  Iš matricos (ji buvo gauta iš data.frame objekto) buvo pašalintos
    eilutės (chromosomos), kuriose nebuvo nustatytas nei vienas genas
    (visos reikšmės eilutėje buvo lygios 0), ir buvo sukurtos
    stulpelinės diagramos, vizualizuojančios, kokiose chromosomose tam
    tikro geno buvo nustatyta daugiausiai.

Modifikacijos įverčiams atvaizduoti ***heatmap*** pavidalu buvo
pasirinkta chromosoma, turinti didžiausią tam tikro geno raišką.

```{r Task_2_4, echo=FALSE, warning=FALSE}
gene_cmt <- data.matrix(gene_cnts)     # gene_cmt --> gene count matrix
gene_cmt <- gene_cmt[ , 2:length(colnames(gene_cmt))]

# fltrd_gene_cnts --> filtered gene count matrix (rows that contained only
# only zero values were removed. It means that gene expression in
# particular chromosome was not identified):
fltrd_gene_cnts <- gene_cmt
rownames(fltrd_gene_cnts) <- c(chr_names)
fltrd_gene_cnts <- fltrd_gene_cnts[rowSums(fltrd_gene_cnts[]) > 0, ]
gene_cnts_melt <- reshape2::melt(fltrd_gene_cnts)

plot1 <- ggplot(gene_cnts_melt, aes(x = Var2, y = value))+
            geom_bar(stat = 'identity', fill = "#123e68")+
            ylim(0, 1700) +
            facet_wrap(~Var1) +
            labs(title = "Gene occurrence", x = "Gene",
                 y = "Occurrence") +
            geom_text(aes(label = value), vjust = -0.5, color = "black",
                      size = 3) +
            theme(axis.text.x = element_text(angle = 45, vjust = 0.5),
                  plot.title = element_text(hjust = 0.5, face = "bold"))
print(plot1)
```

Remiantis aukščiau pavaizduotų stulpelinių diagramų duomenimis,
***heatmap*** buvo pasirinkta atvaizduoti **9 chromosomos** ir
***MIR1268A geno*** modifikacijos duomenis.

```{r Task_2_5, warning=FALSE, include=FALSE}
chr <- "chr9"
gene <- "MIR1268A"

# Extracting rows from 'annotation' matrix that describe 8th chromosome:
slctd_chr <- annot_row[annot_row$chr == chr, ] 

# Extracting rows from the 'slctd_chr' matrix which has 'MIR1268A' values
# described in 'UCSC_RefGene_Name' column:
slctd_gene <- slctd_chr[slctd_chr$UCSC_RefGene_Name == gene, ] %>%
                  data.matrix()

# Selecting values from the main modification matrix - beta matrix based on
# the 'slctd_gene' matrix rownames:
beta_slctd_gene <- (beta[rownames(slctd_gene), ])
normalized_beta <- (beta_slctd_gene - rowMeans(beta_slctd_gene)) /
                        rowSds(beta_slctd_gene)

timepoint <- unlist(lapply(colnames(beta_slctd_gene), function(x) {
  if(grepl("t1", x)) '#0D2855'         # dark blue
  else if(grepl("t2", x)) '#3471DA'    # light blue
}))

clonic_lct <- unlist(lapply(colnames(beta_slctd_gene), function(x) {
  if(grepl("dist", x)) '#7D0F1A'       # dark red
  else if(grepl("prox", x)) '#FD2E43'  # light red
}))

aspirin <- unlist(lapply(colnames(beta_slctd_gene), function(x) {
  if(grepl("lng", x)) '#2E5F15'        # dark green
  else if(grepl("nnu", x)) '#5EDA21'   # light green
}))

ages <- samples$age
color_options <- labels2colors(ages)

age_colors <- unlist(lapply(colnames(beta_slctd_gene), function(x) { 
for (age in 1:length(colnames(beta_slctd_gene))) {
  if(grepl(paste0("_", ages[age]), x)) return(color_options[age])
  #print(paste(ages[age], "-->", color_options[age]))
}
}))

annotation_col <- cbind(timepoint, clonic_lct, aspirin, age_colors)
```

\pagebreak

Žemiau atvaizduojami modifikacijos įverčiai ***heatmap*** grafiko
pavidalu prieš atliekant normalizavimą:

```{r Task_2_6, echo=FALSE, message=FALSE, warning=FALSE}
heatmap3(beta_slctd_gene, col = magma(length(beta_slctd_gene)),
         cexRow = 0.8, cexCol = 0.5, margins = c(10, 5), scale = "row",
         ColSideColors = annotation_col, showRowDendro = F,
         main = "Representation of modification values before normalization",
         distfun = function(x) as.dist(1 - cor(t(x), use = "pa")),
         method = "complete", hclustfun = hclust,
         legendfun = function(x) plot.new())

legend("left", legend = c("t1 - first year", "t2 - after 10 years",
                              "distal colon", "proximal colon",
                              "longterm user", "non-user"),
       col = c("#0D2855", "#3471DA", "#7D0F1A", "#FD2E43", "#2E5F15",
               "#5EDA21"), lty= 1, lwd = 10, bty = "n", cex = 0.6,
       title = "Color meanings")
```

\pagebreak

Žemiau atvaizduojami modifikacijos įverčiai ***heatmap*** grafiko
pavidalu atlikus normalizavimą:

```{r Task_2_7, echo=FALSE, message=FALSE, warning=FALSE}
heatmap3(normalized_beta, col = magma(length(normalized_beta)),
         cexRow = 0.8, cexCol = 0.5, margins = c(10, 5), scale = "row",
         ColSideColors = annotation_col, showRowDendro = F,
         main = "Representation of modification values after normalization",
         distfun = function(x) as.dist(1 - cor(t(x), use = "pa")),
         method = "complete", hclustfun = hclust,
         legendfun = function(x) plot.new())

legend("left", legend = c("t1 - first year", "t2 - after 10 years",
                              "distal colon", "proximal colon",
                              "longterm user", "non-user"),
       col = c("#0D2855", "#3471DA", "#7D0F1A", "#FD2E43", "#2E5F15",
               "#5EDA21"), lty= 1, lwd = 10, bty = "n", cex = 0.6,
       title = "Color meanings")
```

Apibendrinus gautus rezultatus galima padaryti išvadą, kad net ir
sumažinus modifikacijos pozicijų skaičių, nagrinėjant tik devintos
chromosomos ***MIR1268A*** geną, negalima įžvelgti aiškių grupių, kurios
gali būti analizuojamos toliau, ieškant bendrumų tarp grupės objektų.

Siekiant gauti ***heatmap*** grafiką, kuriame galima matyti aiškiai
susidariusias grupes (klasterius), buvo nuspręsta iš modifikacijos
įverčių matricos atrinkti 30 variabiliausių pozicijų ir sukurti
**heatmap** grafiką būtent šioms pozicijoms, pavaizduojant ***heatmap***
grafikus prieš atliekant normalizavimą bei atlikus jį.

```{r Task_2_8, echo=FALSE, message=FALSE, warning=FALSE, fig.show='hold', out.width='70%', fig.align = 'center'}
beta_nw <- as.data.frame(beta)
beta_nw$row_var = rowVars(as.matrix(beta_nw[, c(-1)]))

beta_var <- tail(beta_nw[order(beta_nw$row_var), ], n = 30)
beta_var$row_var <- NULL

beta_norm_var <- (data.matrix(beta_var) - rowMeans(data.matrix(beta_var))) /
                              rowSds(data.matrix(beta_var))

col_rnorm = colorRamp2(c(-3, 0, 3), c("#FFFFFF", "#62A5F2", "#073366"))

h1 <- Heatmap(data.matrix(beta_var), cluster_rows = TRUE, col = col_rnorm,
              cluster_columns = FALSE, show_heatmap_legend = FALSE,
              show_column_names = FALSE,
              column_title = "Before standardization")

h2 <- Heatmap(beta_norm_var, cluster_rows = TRUE, col = col_rnorm,
              cluster_columns = FALSE, column_title = "After standardization",
              show_column_names = FALSE,
              heatmap_legend_param = list(title = "Color key"))

plot(h1 + h2)
```

Aukščiau pavaizduotuose ***heatmap*** grafikuose aiškiai matomas spalvų
pasikeitimas, atlikus modifikacijos pozicijų normalizavimą taip, jog
eilučių vidurkis būtų lygus 0, o standartinis nuokrypis būtų lygus 1.

Remiantis "Color key" skale galima pastebėti, jog modifikacijos
pozicijos, kurios prieš atliekant standartizavimą buvo vizualizuojamos
šviesesne spalva, atlikus standartizavimą buvo pavaizduotos dar
šviesesne spalva (reikšmės buvo sumažintos). Tuo tarpu modifikacijos
pozicijos, kurios prieš normalizavimą buvo vaizduojamos tamsesne spalva,
po normalizavimo įgijo dar didesnę reikšmę ir buvo atvaizduotos dar
tamsesne spalva.

\pagebreak

# 3. Tikro ir nuspėto amžiaus bei senėjimo "pagreitėjimo" nustatymas

Žemiau vaizduojami skirtingi **t1** laiko žymės laikrodžiai su
skirtingais žarnų tipais:

```{r Task_3_1, echo=FALSE, message=FALSE, warning=FALSE}
age.all <- DNAmAge(beta, clocks = c("Horvath", "Hannum", "Levine"))
# Creating objects that save needed 'timepoint' and 'clonic location' data:
vec <- colnames(beta)
info.clean <-samples[samples$graph_name %in% vec, ]
info.t1 <- info.clean[info.clean$timepoint == "t1", ]
info.t1.proximal <- info.t1[info.t1$clonic_location == "proximal", ]
info.t1.distal <- info.t1[info.t1$clonic_location == "distal", ]
info.t2 <- info.clean[info.clean$timepoint == "t2", ]
info.t2.proximal <- info.t1[info.t1$clonic_location == "proximal", ]
info.t2.distal <- info.t1[info.t1$clonic_location == "distal", ]

# Using age.all object to extract 't1'/'t2' and 'distal'/'proximal':
vec.t1 <- info.t1.proximal$graph_name
age.t1.proximal <- age.all[age.all$id %in% vec.t1, ]
vec.t1 <- info.t1.distal$graph_name
age.t1.distal <- age.all[age.all$id %in% vec.t1, ]
vec.t2 <- info.t2.proximal$graph_name
age.t2.proximal <- age.all[age.all$id %in% vec.t2, ]
vec.t2 <- info.t2.distal$graph_name
age.t2.distal <- age.all[age.all$id %in% vec.t2, ]

# Representing the clocks with different timepoints and clonic locations.
# Timepoint - t1:
p1.Horvath.proximal <- plotDNAmAge(age.t1.proximal$Horvath, info.t1.proximal$age,
                                   tit = "Horvath method (proximal)")
p1.Hannum.proximal <- plotDNAmAge(age.t1.proximal$Hannum, info.t1.proximal$age,
                                  tit = "Hannum method (proximal)")
p1.Levine.proximal <- plotDNAmAge(age.t1.proximal$Levine, info.t1.proximal$age,
                                  tit = "Levine method (proximal)")

p1.Horvath.distal <- plotDNAmAge(age.t1.distal$Horvath, info.t1.distal$age,
                                 tit = "Horvath method (distal)")
p1.Hannum.distal <- plotDNAmAge(age.t1.distal$Hannum, info.t1.distal$age,
                                tit = "Hannum method (distal)")
p1.Levine.distal <- plotDNAmAge(age.t1.distal$Levine, info.t1.distal$age,
                                tit = "Levine method (distal)")

p1.Horvath.proximal + p1.Horvath.distal + p1.Hannum.proximal +
  p1.Hannum.distal + p1.Levine.proximal + p1.Levine.distal +

  plot_layout(ncol = 2)
```

\pagebreak

Žemiau vaizduojami skirtingi **t2** laiko žymės laikrodžiai su
skirtingais žarnų tipais:

```{r Task_3_2, echo=FALSE, message=FALSE, warning=FALSE}
# Timepoint - t2:
p2.Horvath.proximal <- plotDNAmAge(age.t2.proximal$Horvath, info.t2.proximal$age,
                                   tit = "Horvath method (proximal)")
p2.Hannum.proximal <- plotDNAmAge(age.t2.proximal$Hannum, info.t2.proximal$age,
                                  tit = "Hannum method (proximal)")
p2.Levine.proximal <- plotDNAmAge(age.t2.proximal$Levine, info.t2.proximal$age,
                                  tit = "Levine method (proximal)")

p2.Horvath.distal <- plotDNAmAge(age.t2.distal$Horvath, info.t2.distal$age,
                                 tit = "Horvath method (distal)")
p2.Hannum.distal <- plotDNAmAge(age.t2.distal$Hannum, info.t2.distal$age,
                                tit = "Hannum method (distal)")
p2.Levine.distal <- plotDNAmAge(age.t2.distal$Levine, info.t2.distal$age,
                                tit = "Levine method (distal)")

p2.Horvath.proximal + p2.Horvath.distal + p2.Hannum.proximal +
  p2.Hannum.distal + p2.Levine.proximal + p2.Levine.distal +

  plot_layout(ncol = 2)
```

\pagebreak

```{r Task_3_3, echo=FALSE, message=FALSE, warning=FALSE}
knitr::opts_chunk$set(fig.height = 13, fig.width = 16)
age.GSE <- DNAmAge(beta, clocks = c("Horvath", "Hannum","Levine"),
                   age = info.clean$age)

Acc <- data.frame(Horvarth = abs(mean(age.GSE$ageAcc.Horvath)),
                   Hannum = abs(mean(age.GSE$ageAcc.Hannum)),
                   Levine = abs(mean(age.GSE$ageAcc.Levine)))

# The most precise clock:
MostAcc <- colnames(Acc)[max.col(-Acc, ties.method = "first")]

# The least precise clock:
LeastAcc <- colnames(Acc)[max.col(Acc, ties.method = "first")]
```

Parodoma, kuris iš trijų naudojamų laikrodžių yra tiksliausias ir kuris
yra mažiausiai tikslus. Naudododami *DNAmAge* funkciją pamatome, kad
tiksliausias yra ***Horvarth laikrodis***, kurio nuklydimo vidurkis yra
***`r MostAcc`***, o mažiausiai tikslus yra ***Levine laikrodis***,
kurio vidurkis yra lygus ***`r LeastAcc`***.

Žemiau yra vaizduojama koreliacija tarp visų trijų laikrodžių.

```{r Task_3_4, echo=FALSE, message=FALSE, warning=FALSE}
plotCorClocks(age.GSE)
```

\pagebreak

Žemiau esančiuose grafikuose vaizduojama, kaip nuspėti amžiai skiriasi
tarp aspiriną vartojančių ir jo nevartojančių.

```{r Task_3_5, echo=FALSE, message=FALSE, warning=FALSE}
tit = "Horvath method"
xx = "DNA Methylation Age"
yy <- "Chronological Age"
my.formula <- y ~ x

df <- data.frame(x = age.t1.proximal$Horvath, y = info.t1.proximal$age,
                 aspirin = info.t1.proximal$aspirin)

p1 <- ggplot(data = df, aes(x = x, y = y)) + 
          geom_smooth(method = "lm", se = FALSE, color = "black",
                      formula = my.formula) +
          labs(title = tit, x = xx, y = yy) +
          aes(x, y, colour = aspirin) +
          geom_point()

df <- data.frame(x = age.t2.proximal$Horvath, y = info.t2.proximal$age,
                 aspirin = info.t2.proximal$aspirin)

p2 <- ggplot(data = df, aes(x = x, y = y)) + 
          geom_smooth(method = "lm", se = FALSE, color = "black",
                      formula = my.formula) +
          labs(title = tit, x = xx, y = yy) +
          aes(x, y, colour = aspirin) +
          geom_point()

plot(p1 + p2)
```

\pagebreak

## Papildoma užduotis

### t-SNE klasterizavimo metodo pritaikymas

Kaip ir pirmosios užduoties atveju diagramoje matome du ryškius
klasterius: vieną su distalinės žarnos mėginiais (apskritimai), kitą su
proksimalinės žarnos mėginiais (trikampiai).

```{r Extra_1, echo=FALSE, message=FALSE, warning=FALSE, fig.show='hold'}
beta_meta <- samples %>%
    select(clonic_location, timepoint, aspirin) %>% mutate(ID=row_number())

tSNE_fit <- Rtsne(d)

# The t-SNE result object contains two t-SNE components that we are
# interested in. We can extract the components and save it in a dataframe:
tSNE_df <- tSNE_fit$Y %>% 
  as.data.frame() %>%
  rename(tSNE1 = "V1", tSNE2 = "V2") %>%
  mutate(ID = row_number())

tSNE_df <- tSNE_df %>%
  inner_join(beta_meta, by = "ID")

tSNE_df %>%
  ggplot(aes(x = tSNE1, y = tSNE2, color = aspirin, shape = clonic_location)) +
      geom_point(size = 3) +
      ggtitle("t-SNE") +
      theme(legend.position = "right")
```