Язык R и его применение в биоинформатике

class: center, middle, inverse, title-slide

# Язык R и его применение в биоинформатике
### Анна Валяева
### 17.12.2021

---

class: inverse, center, middle

# scRNA-seq

# РНК-секвенирование единичных клеток

---
# Bulk *vs* single cell

.small[*[Benchmarking scRNA-seq imputation tools with respect to network inference highlights deficits in performance at high levels of sparsity](https://doi.org/10.1101/2021.04.02.438193)*]

---
# История scRNA-seq

.small[*[Exponential scaling of single-cell RNA-seq in the past decade](https://www.nature.com/articles/nprot.2017.149)*]

---
# Разделение клеток

.small[*[Single-cell RNA sequencing technologies and bioinformatics pipelines](https://www.nature.com/articles/s12276-018-0071-8)*]

---
# Разделение клеток

- droplet-based методы (10x Chromium, Drop-seq)

- plate-based методы (Smart-seq2)

.small[*[A Single-Cell Transcriptional Roadmap of the Mouse and Human Lymph Node Lymphatic Vasculature](http://dx.doi.org/10.3389/fcvm.2020.00052)*]

---
# Каждой клетке - по баркоду!

## Кажому транскрипту - по UMI [Unique Molecular Identifier]*

*Не во всех протоколах используется UMI

.small[*[Single-cell RNA sequencing technologies and bioinformatics pipelines](https://www.nature.com/articles/s12276-018-0071-8)*]

---
# Препроцессинг данных

## От прочтений до экспрессии

- картирование прочитанных фрагментов на референс
- приписывание прочтений генам
- приписывание прочтений клеткам по баркодам (демультиплексирование) 
- подсчет уникальных молекул РНК (дедупликация UMI)

Для препроцессинга данных 10x используется Cell Ranger.

---
# Cell Ranger

---
# Препроцессинг данных

## От экспрессии до кластеров клеток

.left-column[
- [Seurat](https://satijalab.org/seurat/index.html)
- [Scanpy](https://scanpy.readthedocs.io/en/stable/)
- [Scater](https://bioconductor.org/packages/release/bioc/vignettes/scater/inst/doc/overview.html)
- ...
]

.right-column[
<img src="data:image/png;base64,#img/2021-12-17/pipeline.png" width="130%" style="display: block; margin: auto;" />
]

---
# Импорт данных

```r
library(Seurat)

# PBMC датасет
# В папке 3 файла: barcodes.tsv, genes.tsv, matrix.mtx
pbmc.data <- Read10X(data.dir = "data/pbmc3k/filtered_gene_bc_matrices/hg19/")
```

---
# Seurat объект

```r
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc
```

```
An object of class Seurat 
13714 features across 2700 samples within 1 assay 
Active assay: RNA (13714 features, 0 variable features)
```

---
# Контроль качества

- число детектированных генов (nFeature):
  - мало - плохие клетки или пустые капли
  - много - 2 или более клеток в капле (дублеты) 
- число транскриптов (nCount)
- процент прочтений от митохондриальных генов:
  - умирающие клетки

```r
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
```

---
# QC графики

```r
VlnPlot(pbmc, features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), ncol = 3)
```

```r
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
```

---
# Нормализация

`\(ncount_{ij} = log1p(\frac{count_{ij}}{\sum_{i=1}^{n}count_{ij}}10^4)\)`

```r
pbmc <- NormalizeData(pbmc, normalization.method = "LogNormalize", scale.factor = 10000)
```

Нормализованная экспрессия хранится в отдельном "слоте" - `pbmc[["RNA"]]@data`.

Альтернатива - [sctransform](https://satijalab.org/seurat/articles/sctransform_vignette.html).

---
# Feature seelction

## Отбираем гены с наиболее вариабельной экспрессией (highly variable features)

```r
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)

# топ 10
top10 <- head(VariableFeatures(pbmc), 10)

LabelPoints(VariableFeaturePlot(pbmc), points = top10, repel = TRUE)
```

---
# Scaling

## перед пониженим размерности

```r
pbmc <- ScaleData(pbmc)
```

---
# Понижение размерности

## PCA

```r
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))

DimPlot(pbmc, reduction = "pca")
```

---
# Размерность датасета

```r
ElbowPlot(pbmc)
```

---
# Кластеризация

```r
# строим граф k-ближайших соседей
pbmc <- FindNeighbors(pbmc, dims = 1:10)

# ищем кластеры
pbmc <- FindClusters(pbmc, resolution = 0.5)
```

```
Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck

Number of nodes: 2638
Number of edges: 95965

Running Louvain algorithm...
Maximum modularity in 10 random starts: 0.8723
Number of communities: 9
Elapsed time: 0 seconds
```

---
# t-SNE / UMAP

```r
pbmc <- RunUMAP(pbmc, dims = 1:10)

DimPlot(pbmc, reduction = "umap")
```

---
# Аннотация кластеров

ДЭ ген в кластере определяется по разнице в среднем уровне экспрессии этого гена внутри кластера относительно остальных кластеров.

```r
# ищем маркеры для всех кластеров
pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
pbmc.markers %>%
  group_by(cluster) %>%
  slice_max(n = 2, order_by = avg_log2FC)
```

```
# A tibble: 18 x 7
# Groups:   cluster [9]
       p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene    
       <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>   
 1 1.17e- 83       1.33 0.435 0.108 1.60e- 79 0       CCR7    
 2 1.74e-109       1.07 0.897 0.593 2.39e-105 0       LDHB    
 3 0               5.57 0.996 0.215 0         1       S100A9  
 4 0               5.48 0.975 0.121 0         1       S100A8  
 5 7.99e- 87       1.28 0.981 0.644 1.10e- 82 2       LTB     
 6 2.61e- 59       1.24 0.424 0.111 3.58e- 55 2       AQP3    
 7 0               4.31 0.936 0.041 0         3       CD79A   
 8 9.48e-271       3.59 0.622 0.022 1.30e-266 3       TCL1A   
 9 4.93e-169       3.01 0.595 0.056 6.76e-165 4       GZMK    
10 1.17e-178       2.97 0.957 0.241 1.60e-174 4       CCL5    
11 3.51e-184       3.31 0.975 0.134 4.82e-180 5       FCGR3A  
12 2.03e-125       3.09 1     0.315 2.78e-121 5       LST1    
13 6.82e-175       4.92 0.958 0.135 9.36e-171 6       GNLY    
14 1.05e-265       4.89 0.986 0.071 1.44e-261 6       GZMB    
15 1.48e-220       3.87 0.812 0.011 2.03e-216 7       FCER1A  
16 1.67e- 21       2.87 1     0.513 2.28e- 17 7       HLA-DPB1
17 3.68e-110       8.58 1     0.024 5.05e-106 8       PPBP    
18 7.73e-200       7.24 1     0.01  1.06e-195 8       PF4     
```

---
# Аннотация кластеров

- найденные ДЭ гены вручную сравниваются с известными маркерами клеточных типов

```r
FeaturePlot(pbmc, features = c("MS4A1", "LYZ"))
```

---
# Аннотация кластеров

- автоматическая аннотация кластеров и использованием референсного датасета, например [Azimuth](https://azimuth.hubmapconsortium.org/)

---
# Аннотация кластеров

```r
new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono", "NK", "DC", "Platelet")
names(new.cluster.ids) <- levels(pbmc)
pbmc <- RenameIdents(pbmc, new.cluster.ids)
DimPlot(pbmc, reduction = "umap", label = TRUE, pt.size = 0.5) + NoLegend()
```

---
# Гармонизация данных

## Batch effect correction

.left-column[
- [CCA Seurat](https://satijalab.org/seurat/articles/integration_introduction.html)
- [RPCA Seurat](https://satijalab.org/seurat/articles/integration_rpca.html)
- [Harmony](https://portals.broadinstitute.org/harmony/articles/quickstart.html)
- [LIGER](https://github.com/welch-lab/liger)
- [conos](https://github.com/kharchenkolab/conos)
- ...
]

.right-column[
<img src="data:image/png;base64,#img/2021-12-17/batch.PNG" width="90%" style="display: block; margin: auto;" />
]

---
# Что дальше?

.small[*[Current best practices in single-cell RNA-seq analysis: a tutorial](https://www.embopress.org/doi/full/10.15252/msb.20188746)*]

---
# Коммуникация между клетками

.left-column[
- [CellChat](http://www.cellchat.org/)
- [CellPhoneDB](https://github.com/ventolab/CellphoneDB)
- [ICELLNET](https://github.com/soumelis-lab/ICELLNET)
- ...
]

.right-column[
<img src="data:image/png;base64,#img/2021-12-17/cellchat.PNG" width="100%" style="display: block; margin: auto;" />

.small[*[Deciphering cell–cell interactions and communication from gene expression](https://www.nature.com/articles/s41576-020-00292-x)*]

]

---
# Импутация данных

**Задача:** избавиться от многочисленных нулей в матрице экспрессии и реконструировать "реальные" профили экспрессии генов.

.pull-left[
## Откуда нули?

- **Биологические причины** (настоящие “нули”): случайные флуктуации в экспрессии гена
- **Технические причины:** низкий процент молекул детектируется (5-15% транскриптома) - ***dropout effect***

Существует много методов импутации (восстановления нулей).
]

.pull-right[
<img src="data:image/png;base64,#img/2021-12-17/dropout.PNG" width="100%" style="display: block; margin: auto;" />

.small[*[Bayesian approach to single-cell differential expression analysis](https://www.nature.com/articles/nmeth.2967)*]

]

---
# Методы импутации

.small[*[A systematic evaluation of single-cell RNA-sequencing imputation methods](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-020-02132-x)*]

---
# MAGIC

.pull-left[
**Идея:** по мере “дифференцировки” клеток экспрессия генов согласованно изменяется, поэтому нужно хитро искать “соседей” и, ориентируясь на них, восстанавливать значения экспрессии.

]

.pull-right[
<img src="data:image/png;base64,#img/2021-12-17/magic.jpg" width="70%" style="display: block; margin: auto;" />
]

.small[*[Recovering Gene Interactions from Single-Cell Data Using Data Diffusion](https://doi.org/10.1016/j.cell.2018.05.061)*]

---
# Деконволюция bulk RNA-seq данных

---
# Пространственная транскриптомика

.small[*[Spatially Resolved Transcriptomes—Next Generation Tools for Tissue Exploration](http://dx.doi.org/10.1002/bies.201900221)*]

---
# Мультиомиксные данные

.small[*[Single-cell multiomics: technologies and data analysis methods](https://doi.org/10.1038/s12276-020-0420-2)*]

---
# Что почитать и посмотреть

- [Сайт Seurat и множество туториалов](https://satijalab.org/seurat/index.html)
- [Analysis of single cell RNA-seq data](https://www.singlecellcourse.org/index.html)
- [New Advances in Single-Cell and Spatial Genomics (2021)](https://youtu.be/uFRe-UMUlMQ)
- [Current best practices in single-cell RNA-seq analysis: a tutorial](https://www.embopress.org/doi/full/10.15252/msb.20188746)