Single cell RNA seq notes

0.sigle_cell_rnaseq_analysis.Rmd

---
output:
  pdf_document: default
  html_document: default
---
scRNA-seq is useful for detecting changes in cell sub-populations in different conditions (or detecting sub-populations de novo), but not much else. It can also be used to detect some molecular mechanisms, but not really. Please don't do random experiments hoping to see "something".

There must be multiple biological samples from the population of interest, each with the many number of measured cells.

There are many types of single cells analyses, (e.g. "tissue-wide", spatial, which resolves the cells based on their position in the tissue, etc...), so you must really be sure of what the target is (the entity that needs to be measured), again, the bio question needs to be extremely clear.

RNA quality for sc is really important, and samples (ligated to beads) cannot be frozen before sc, so you cannot re-analyze the same batch. Which samples are included in the single run is very important.

---
You get multiple output files from the seq. You seq a barcode for each cell and just a tiny portion of the 3' of the RNA transcript. You can then map these reads to the transcripts, and deconvolute the cells. The program used this is `cellranger`.
Warning: cellranger devs delete the old versions - so docker is highly recommended.
Cellranger basically uses STAR and does some quality checks.
Use version 3 or later since they increased sensitivity a lot.

A scRNA seq chip can sequence ~10k cells per well, and you get 4 wells. Each well really has 100k beads, but most have to be empty as we need to get beads with just a single cell each. Each chip costs 10k ish euro. Due to the very high cost you can multiplex the samples with RNA oligo tags bound to lipids. You mix the cells from each sample with one type of tag, and that gets sequenced. So you can demultiplex a single sample from the bulk of the cells.

This barcoding can also increase the % of occupied beads in each chip: now we can detect if multiple cells got connected to the same bead, and we can discard these beads.

The minimum num of reads per cell is 25k. A good amount of reads to see all genes is 50k. So the RNA-seq after the bead ligation is still really expensive (often ~2k per sample).

Cellranger discards all non-exon reads, so we lose a ton of them. (? maybe chromium is the quality control tool, similar to fastQC?). If a lot of reads map to introns, we might have gotten DNA contamination, so the experiment might have technically failed.

---
Output count matrixes are very highly 0-inflated, so we often defalte them to a sparse matrix (usually .mtx). You have three output files:
- matrix.mtx (with the counts)
- barcodes (with the cell barcodes)
- genes (with ensembl genes)

The MTX file looks like this:
```
#_of_rows #_of_columsn #_of_cells
row_# col_# count
row_# col_# count
row_# col_# count
row_# col_# count
row_# col_# count
...
```

In later versions the "genes" table is named "features", with gene names included to other metadata, like ATAC-Seq, called peaks, etc etc

If you have an old "genes" table, you just need to add a new col with the ID "Gene Expression" to tell generanger that those are gene expression data

## Converting from sparse to dense matrix
You can go back from sparse to dense matrices with `h5tocsv` in the library `rCASC``. This is the library that wraps docker calls for greater reproducibility.
Reproducibility is essential since the downstream scRNA-seq analyses may take multiple years to be performed, and we might need to re-run the initial analyses.

```{r}
library(rCASC) # This must be a library call... :(
# I'm gonna copy-paste this lib call in every chunk just so we can run them
# independently.

here <- "~/data/scRNAseq_v6.1/exercise1/GSM4679532_P01"
# This has very basic input data

rCASC::h5tocsv(
  group = "docker",
  file = file.path(here, "matrix.mtx.gz"),
  type = "10xgenomics",
  version = 5
)

# This fetches the docker container automatically, with the version that you specify.
```

rCASC makes a `container.txt` file with the list of used containers.
WE should have gotten a `matrix.csv`. This file has ~500 cells, and it's already 41 Mb!
For ref, the compressed mtx was 3.2 Mb.

## Cell quality

To evaluate cell quality, we count the number of reads of ribosomal genes, and mitochondrial genes. In short, ribo genes should be stable, and mito genes should be low. Apoptosis causes upregulation of mito genes. Less than ~10% of all detected genes should be mito (in human, 5% for mouse).

To say that some gene x is expressed and active in this cell, Calogero suggests to use a thresholds of at least 3 reads on that gene.
We can therefore map each cell by 3 dimensions: # of called genes (with 3 reads or more), % of mito and % of ribo.
Cells with low # of called genes, with too high ribo or too high mito are discarded, as the low information content just makes them cluster together - in a meaningless cluster.
If a lot of cells are high in mito, the tissue might not be normal. This might be expected - the cells might be tumors, or stressed, or treated to cause apoptosis.

So, choose some good thresholds and filter the calls (from all samples) with them, so we have uniform populations of cells from each sample.

```{r}
library(rCASC)
here <- "~/data/scRNAseq_v6.1/exercise2"

# This data is old, so we must use an old, compatible annotation.

rCASC::mitoRiboUmi(
  group = "docker",
  # The scratch folder is useful if you need to use a fast disk while working,
  # and using a slow disk afterwards.
  scratch.folder = "~/scratch",
  file = file.path(here, "matrix.csv"), # the decompressed matrix
  separator = ",",
  gtf.name = file.path(here, "Homo_sapiens.GRCh38.99.gtf"),
  # This is the ENSEMBL biotype of interest, e.g. protein_coding
  bio.type = "protein_coding",
  # This is the reads per gene parameter from before. 3 works well.
  umiXgene = 3
)
```

Why ensembl? They are much more reproducible as they keep all their assemblies online, forever (probably). Plus, they distribute the gene symbols in their GTF files, so they are stable in time. In the new rCASC, this is no longer needed as they use info from the previous cellranger runs.

If you have to filter, filter for low counts of genes mainly. High mito genes might be phisiologically relevant.

We should have gotten a `Ribo_mito.pdf` file with the plot (I've downloaded it).

In one of their multiplexed experiments, they noticed some samples that were very bad, and some samples that were very good. Why did all the bad samples cluster in the same barcode on the same chip? They didn't know. Since the experimenters only wanted healthy cells, they just discarded these bad runs, but we can't do this a priori.

When demultiplexing, consider keeping the real origin of the sample in its barcode, so you never lose the original batch (which might carry an effect).

Sometimes you even see clustering at the level of % of ribo and % of mito!

## Annotation and Filtering
Should we filter out mitochondrial and ribosomal genes? They are the ones that we see the most, so they probably mask the real effects that we want to see.
However, sometimes the mito/ribo genes *are* interesting, and we want to keep it.
Think about the experimental question.

You can even cluster based on mito/ribo genes first, separate the subpopulations, and *then* filter them out, running downstream analyses.

```{r}
library(rCASC)
here <- "~/data/scRNAseq_v6.1/exercise3"

rCASC::scannobyGtf(
  group = "docker",
  file = file.path(here, "matrix.csv"),
  gtf.name = "Homo_sapiens.GRCh38.99.gtf",
  biotype = "protein_coding",

  # Skip removing all mito genes?
  mt = FALSE,
  # Skip removing all ribo genes?
  ribo.proteins = FALSE,
  # Same as before
  umiXgene = 3,

  # Filter out the % of ribo and % of mito
  # Keep only cells inside this range.
  riboStart.percentage = 0,
  riboEnd.percentage = 100,

  mitoStart.percentage = 0,
  mitoEnd.percentage = 100,

  # Keep only cells with at least this # of called genes
  thresholdGenes = 100
)

```

This makes three files:
- matrix_annotated_genes.pdf with plots of the results of filtering (which I have downloaded)
- filtered_annotated_matrix.csv with the filtered and annotated matrix.csv (:))
- `filteredStatistics.txt`

## Seurat cell cycle

Seurat can try to infer the cell cycle of cells. Why is this important? If cells are duplicating, the signal from the cell cycle is strong, and clustering can be biased based on that. So, colouring cells based on their cycle is important when detecting cell clusters.

Plus, it cn be interested - e.g. undifferentiated, staminal cells might be more actively dividing rather than the rest.

Why seurat? It's one of the few that does this, and it seems to work. How does it do it? The code is black magic, so we don't actually know.

```{r}
library(rCASC)
here <- "~/data/scRNAseq_v6.1/exercise4"

# The args are self-explainatory.
rCASC::seurat_ccycle(
  group = "docker",
  scratch.folder = "~/scratch",
  file = file.path(here, "filtered_annotated_matrix.csv"),
  separator = ",",
  seed = 0
)
```
As output, we get a matrix (`filtered_annotated_matrix_cellCycle`) annotated with the predicted cell cycle and the corresponding column name in the original matrix.

Seurat never fails - it just guesses if something goes horribly wrong. Nice (?).
We usually make a new count matrix with modified colnames that have the barcode PLUS the cell cycle.

```{r}
here <- "~/data/scRNAseq_v6.1/exercise4"

original <- read_csv(
  file.path(here, "filtered_annotated_matrix.csv"), lazy = TRUE,
  show_col_types = FALSE
) %>%
  column_to_rownames("...1")
cycle <- read_csv(
  file.path(here, "filtered_annotated_matrix_cellCycle.csv"),
  col_names = c("label", "cycle"),
  show_col_types = FALSE
) %>%
  column_to_rownames("label")

# Paste the new colnames
new_colnames <- paste0(
  colnames(original), "_",
  cycle[colnames(original), "cycle"]
)
# Check visually that everything went ok
cat("\n") # Get off my back
colnames(original) %>% head() %>% paste(collapse = " ") %>% cat()
cat("\n")
new_colnames %>% head() %>% paste(collapse = " ") %>% cat()

colnames(original) <- new_colnames

# use write.csv to keep row names in a col with no name.
write.csv(original, file = file.path(here, "preprocessed_matrix.csv"))

rm(list = c("cycle", "original", "new_colnames"))
```

Note the '.1' in the colnames. This is from cellranger. It is useful to discern different cells in different experiments that are bundled together, and thus might have colliding barcodes.

## Clustering
A simple PCA on raw counts usually doesn't lead to clusters. Normally, raw counts are first log2, then projected and clusters are detected.

However - how do we do dimensionality reduction with the zero-inflated data? We need an *inputation method*, of which there are many.
One of them is SAVER, which seems to be slightly better than the competition.
There is some evidence that imputation doesn't do much, and only SAVER gives some improvement on the actual clustering.

Like DESeq2, it uses gene-level information to infer the real counts (instead of the zeroes).
Downside is that the matrix is now not sparse anymore, so we can't compress it.

Differences between populations may be masked by larger variances, like differences in the actual patient. This can be overstepped by some sort of post-doc test after clustering.
To do this, generally the n clusters are bootstrapped and DEAs are run between the bootstrapped samples. Then, the number of gens that come out after ever DEA is checked, and only those genes that are found over and over again as DEGs are considered as "true" degs.

Many clustering algos are based on K-means (or similar) on the transformed space (or similar).
The usual issue is to set a meaningful K. The idea they had was to run k-meres with some k. Then, delete 10% of cells and run it again. This bootstrapping is repeated many times (~40).
If the clusters in which each cells are more or less the same, then the k is good. If not, the k is probably bad. A "good" cluster is one that has ~50% of cells which always fall in the same cluster (i.e. "stable cells").

--> Skip to exercise 7

Clustering and especially this bootstrapping method is really computationally expensive, both in CPUs and RAM.

```{r}
library(rCASC)
here <- "~/data/scRNAseq_v6.1/exercise7"

# How many PCA dimensions should we consider for clustering?
# This function shows us the SCREE plot
rCASC::seuratPCAEval(
  group = "docker",
  scratch.folder = "~/scratch",
  file = file.path(here, "annotated_BE1500.csv"),
  separator = ",",
  # Is the matrix already log 10? 0 -> no, 1 -> yes
  logTen = 0,
  # Seed of random number generator
  seed = 0,
  # is the matrix in sparse format?
  sparse = FALSE,
  # The matrix is dense, so this is NULL
  format = "NULL"
)

```
This takes a very, very long time due to the action of SAVER (probably).
This makes a scree plot, basically, so we can check where the variance of the PCs stops dropping too much. The choice is still arbitrary, tho.
To me, 8 sounds good.

We can now run clustering
```{r}
library(rCASC)
here <- "~/data/scRNAseq_v6.1/exercise8"

rCASC::seuratBootstrap(
  # From here ---->
  group = "docker",
  scratch.folder = "~/scratch",
  file = file.path(here, "annotated_BE1500.csv"),
  separator = ",",
  logTen = 0,
  seed = 0,
  sparse = FALSE,
  format = "NULL",
  # <---- To here, same as before
  # how many bootstraps?
  nPerm = 4,
  # How many permutations?
  permAtTime = 1,
  # How many cells should be deleted in each bootstrap? (%)
  percent = 10,
  # How many PCA dimensions to use? We determined this in the previous step.
  pcaDimensions = 8,
  # Seurat resolution -> the smaller, the smaller the clusters, and so we have a
  # more granular check on the k
  resolution = 0.8
)


```

We can recreate the plot we got in the previous step easily as it saved a file with the xy coordinates of each point. This way, we can layer on more information.

```{r}
library(rCASC)
here <- "~/data/scRNAseq_v6.1/exercise9"

# We need to re-run the cell cycle analysis, so we can label the plot with them
rCASC::seurat_ccycle(
  group = "docker",
  scratch.folder = "~/scratch",
  file = file.path(here, "annotated_BE1500.csv"),
  separator = ",",
  seed = 0
)

```
```{r}
# We now add this information on the "plotting coordinates" file:
coords <- read_csv(
  file.path(here, "Results/annotated_BE1500/6/annotated_BE1500_clustering.output.csv"), lazy = TRUE,
  show_col_types = FALSE
) %>%
  select(all_of(c("cellName", "Belonging_Cluster", "xChoord", "yChoord")))
cycle <- read_csv(
  file.path(here, "annotated_BE1500_cellCycle.csv"),
  col_names = c("label", "cycle"),
  show_col_types = FALSE
)

coords <- merge(coords, cycle, by.x = "cellName", by.y = "label")
coords$cell_type <- sapply(coords$cellName, function(x){str_split_i(x, "_", 2)})

# We can now plot again, with our added information
ggplot(coords, aes(x = xChoord, y = yChoord)) +
  geom_point(aes(color = cell_type), size = 0.8) +
  theme_minimal()


rm(list = c("cycle", "original", "new_colnames"))

```
See how the giant central cluster is actually two cell lines? It's basically impossible to separate. The green points near the pink cluster on the left, as well as the weird color dots in otherwise uniform clusters are probably technical artifacts: the demultiplexing process did not go as intended.

```{r}
ggplot(coords, aes(x = xChoord, y = yChoord)) +
  geom_point(aes(color = cycle), size = 0.8) +
  theme_minimal()
```

See how some clusters are "stripey" based on the cell cycle?

## Making pseudo-bulks

Sometimes we want to make a dataset that is pseudo-bulked from a cluster of cells, to run some differential expression analysis, for instance. There are functions for that.

```{r}
library(rCASC)
here <- "~/data/scRNAseq_v6.1/exercise9"

rCASC::bulkClusters(
  group = "docker",
  scratch.folder = "~/scratch",
  file = file.path(here, "annotated_BE1500.csv"),
  separator = ",",
  cl = file.path(here, "Results/annotated_BE1500/6/annotated_BE1500_clustering.output.csv")
)
```

This makes a series of files with the pattern `annotated_BE1500_bulk*`: `Column` has one column per cluster, with Z scores computed in each column, `Row` is similar to columns, but with Z-scores calculated on rows, and `log2`, with simply summed counts in that bulk cluster, scaled as in `log2(x+1)`.

---
Why would we need to run DEAs on the clusters? Well, you often want to:
- Check if the clustering worked well;
- Find genes that are "signature" for that cluster;
- Find genes that can be then used in other tech (like FACS) to sort out the cells in that cluster;

---

## Splitting dataset based on the cluster
We might want to keep the single cell identity but just split the data based on the clusters, for downstream analysis (e.g. subclustering).

```{r}
library(rCASC)
here <- "~/data/scRNAseq_v6.1/exercise9"

rCASC::splitClusters(
  group = "docker",
  scratch.folder = "~/scratch",
  file = file.path(here, "annotated_BE1500.csv"),
  separator = ",",
  nCluster = 6
)
# This expects the "Results/**/6/" folder, and it automatically finds the output data
```

This makes the Results/**/6/ClusterSplitted folder with each splitted matrix.

## Finding marker genes in the clusters
A tool that can find these marker genes in each cluster is COMET. It is based on the minimal hypergeometric test.
I don't reallyb understand the how, but whatever.

COMET might die for some clusters due to ~Inf memory requirements if it cannot find suitable candidate gene markers. This is informative: this cluster might be very heterogeneous, and might be worthwhile sub-clutering this cluster.

COMET also checks for "negated" genes: genes that are expressed in all other clusters except this one. This is also a sign of non-informative cluster, and you might need to subcluster.

```{r}
library(rCASC)
here <- "~/data/scRNAseq_v6.1/exercise9"

rCASC::cometsc(
  group = "docker",
  scratch.folder = "~/scratch",
  file = file.path(here, "annotated_BE1500.csv"),
  separator = ",",
  # Percentage of genes marked as "expressed" that must be in the top of the list
  # it helps filter out artifacts. 0.15 is good.
  X = 0.15,
  # Number of combinations of genes to be considered. Increasing this is immensly expensive
  # 2 is good
  K = 2,
  # Are these raw counts?
  counts="True", # Needs a string for some reason...
  # Skip making plots? Makes this faster.
  skipvis = "False",
  # This finds the iput dir automatically, like the other call.
  nCluster = "6",
  # How many threads to spawn? It's good to set this equal to nCluster
  threads = 6
)
```

COMET takes AGES - often you need to let it run overnight.
The pre-calculated output is in exercise 10: Results/**/6/outputdata (for gene lists) and outputviz for plots. The singletons are one-gene markers. The pairs are the 2 (k=2) genes that, when combined, successfully mark a single cluster.

The "True positive" is the rate of this gene in this cluster. The "True negative" is the rate of this gene not inside this cluster. The optimal is to have 1 in both TP and TN (the gene is in all cells in this cluster AND it's never found in any other cell).


With the output from COMET, we often want to make a heatmap to show us if the detected marker genes are actually marking the clusters well.

```{r}
library(rCASC)
here <- "~/data/scRNAseq_v6.1/exercise10"
res_dir <- file.path(here, "Results/annotated_BE1500/6/outputdata")

cluster_data <- list()
for (i in 1:6) {
  cluster_data[[paste0("cluster_", i)]] <- read_csv(
    file.path(res_dir, paste0("cluster_", i, "_singleton_positive_markers_ranked.csv")),
    show_col_types = FALSE
  )
}
# Grab the bulk expression data to make the heatmap

expression <- read_csv(
  "/data/scRNAseq_v6.1/exercise9/annotated_BE1500_bulklog2.csv",
  show_col_types = FALSE
)
```

We can now plot the heatmap:

```{r}
# Find the top n genes that should cluster each cluster well.
n_genes_per_cluster <- 20
top_genes <- lapply(cluster_data, function(x) {
  # these are already sorted
  x$gene_1[1:n_genes_per_cluster]
})

all_genes <- unlist(top_genes)

# Keep just our top genes

top_markers <- expression %>% filter(`...1` %in% all_genes) %>% column_to_rownames("...1")

heatmap(
  as.matrix(top_markers),
  col = heat.colors(256)
)
```


With these lists of genes you can do a lot of interesting things:
- USe some enrichment techniques
- Query graphs of gene relations, like Kinases -> TFs -> downstream genes. We see the downstreams, but we might be able to "go back" to the kinases and TFs that caused these genes to be changed
- We can use these markers and check other data where the cell line is known to try and determine which cell types are actually part of this cluster
- Try to keep in mind the biological question that was posed at the start when handling all kinds of data.
- Detect sub-populations that are not resolved at the level of the initial PCA
- Find only surface markers (for FACS), with the GO

## Annotating genes
There are tools that eat each cluster and try to label it with various features.
One of these tools is scATOMIC, which uses a series of nested random forests to annotate the types of cells in the samples.

Usually, it's best to run annotations on the clusters separately.

```{r}
library(rCASC)
here <- "~/data/scRNAseq_v6.1/exercise10"

rCASC::scAtomic(
  group="docker",
  scratch.folder = "~/scratch",
  file=file.path(here, "annotated_BE1500.csv"),
  separator = ","
)

```


This nets you an excel file with the annotations from all the layers that scAtomic finds, with confidence annotation for each layer.

## Counting cell types
After we got our nice annotations, we can do a bunch of stuff. For instance, we can count the number of cells per type per cluster. This way, we can get an overview of the nature of each cluster.

```{r}
library(rCASC)
here <- "~/data/scRNAseq_v6.1/exercise13"

# Load the cluster-specific data
cluster_data <- list()
for (i in 1:7) {
  d <- read_csv(
    file.path(here, "scAtomic", paste0("wtc1cls", i), paste0("annotation_s1abr_wt_c1_Cluster", i, ".csv")),
    show_col_types = FALSE
  )
  d$cluster <- paste0("cluster_", i)
  d$pan_cancer_cluster <- NULL # Why are you here, just in sample 5?
  cluster_data[[paste0("cluster_", i)]] <- d
}

# Merge them in a big table
big_cell_types <- Reduce(function(x, y) {rbind(x,y)}, cluster_data)

tab <- big_cell_types %>% select(all_of(c("cluster", "layer_5"))) %>% table()

round((tab / rowSums(tab)) * 100, 2) %>% t()
```
We can also plot the types in our PCA plot!

```{r fig.width=5, fig.height=3.5}
coords <- read_csv(
  file.path(here, "s1abr_wt_c1_clustering.output.csv"), lazy = TRUE,
  show_col_types = FALSE
) %>%
  select(all_of(c("cellName", "Belonging_Cluster", "xChoord", "yChoord")))
coords <- merge(coords, big_cell_types, by.x = "cellName", by.y = "...1")

# We can now plot again, with our added information
# This does not work in a loop and I doin't really want to find out why.
ggplot(coords, aes(x = xChoord, y = yChoord)) +
geom_point(aes(color = layer_1), size = 1.2) +
theme_minimal() +
  scale_color_brewer(palette="Set2")

ggplot(coords, aes(x = xChoord, y = yChoord)) +
geom_point(aes(color = layer_2), size = 1.2) +
theme_minimal() +
  scale_color_brewer(palette="Set2")

ggplot(coords, aes(x = xChoord, y = yChoord)) +
geom_point(aes(color = layer_3), size = 1.2) +
theme_minimal() +
  scale_color_brewer(palette="Set3")

ggplot(coords, aes(x = xChoord, y = yChoord)) +
geom_point(aes(color = layer_4), size = 1.2) +
theme_minimal() +
  scale_color_brewer(palette="Set3")

ggplot(coords, aes(x = xChoord, y = yChoord)) +
geom_point(aes(color = layer_5), size = 1.2) +
theme_minimal()

ggplot(coords, aes(x = xChoord, y = yChoord)) +
geom_point(aes(color = layer_6), size = 1.2) +
theme_minimal()

```

We can also see the same heatmap as before:

```{r}
res_dir <- file.path(here, "cometsc/outputdata")

cluster_data <- list()
for (i in 1:6) {
  cluster_data[[paste0("cluster_", i)]] <- read_csv(
    file.path(res_dir, paste0("cluster_", i, "_singleton_positive_markers_ranked.csv")),
    show_col_types = FALSE
  )
}
# Grab the bulk expression data to make the heatmap

expression <- read_csv(
  "/data/scRNAseq_v6.1/exercise13/pseudoBulks/s1abr_wt_c1_bulklog2.csv",
  show_col_types = FALSE
)

# Find the top n genes that should cluster each cluster well.
n_genes_per_cluster <- 10
top_genes <- lapply(cluster_data, function(x) {
  # these are already sorted
  x$gene_1[1:n_genes_per_cluster]
})

all_genes <- unlist(top_genes)

# Keep just our top genes

top_markers <- expression %>% filter(`...1` %in% all_genes) %>% column_to_rownames("...1")

heatmap(
  as.matrix(top_markers),
  col = RColorBrewer::brewer.pal(256, "Spectral")
)
```

## Combining multiple batches
If you have multiple experiments, you might need to merge the batches into a single file (comBAT style). Seurat uses a method where it selects ~2k genes as "anchors" in each computed cluster, and moves the other genes to remove the batch effects.

To use this function:
```{r}
library(rCASC)
here <- "~/data/scRNAseq_v6.1/exercise14"

rCASC::seuratIntegration(
  group = "docker",
  scratch.folder = "~/scratch",
  file1 = file.path(here, "normal", "annotated_GSM3516673_normal.csv"),
  separator1 = ",",
  file2 = file.path(here, "tumor", "annotated_GSM3516672_tumor.csv"),
  separator2 = ",",
  seed = 0,
  # This is the seurat resolution. To start, use the same that you did when you ran
  # the individual clusterings
  k = 0.1
)

```

hum... I think this crashes before it can copy the results.
They are in exercise15, however.

cell_cycle_strat.png

cell_type.png

Ribo_mito.pdf