<figure><img src="image_1.png" /></figure>

HAL
open science

# scFeatureFilter: correlation-based feature filtering for single-cell RNA-seq

Guillaume Devailly

► To cite this version:

Guillaume Devailly. scFeatureFilter: correlation-based feature filtering for single-cell RNA-seq. Workshop: Bioinfo/Biostat, Jan 2018, Toulouse, France. hal-01801978

HAL Id: hal-01801978

https://hal.science/hal-01801978

Submitted on 5 Jun 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
scFeatureFilter:
correlation-based feature filtering for single-cell RNA-seq

@G_Devailly

2018/01/17
Introduction

* Single cell RNA-sequencing is increasingly popular.

* scRNA-seq is noisier than bulk RNA-seq.

* Filtering of noisy, lowly-expressed features* is common.

<footnote>*Feature: gene or transcript</footnote>
* Using spike-in RNA information

* Arbitrary filtering:

[...] on a filtered data set, where we retain only genes with an estimated TPM above 1 in more than 25% of the considered cells. (1)

Genes with less than 5 reads and expressed in less than 10 cells were removed. (2)

Here, low-abundance genes are defined as those with an average count below a filter threshold of 1 [count]. (3)

Genes were filtered, keeping 15,633 out of 26,178 genes that were expressed in at least 5 out of 1,919 sequenced cells (RPKM ≥ 10) and for which cells with expression came from at least two different embryos. (4)

<footnote n="1">1 Soneson & Robinson, bioRxiv, 2017</footnote>

<footnote n="2">2 Stevant et al., bioRxiv, 2017</footnote>

<footnote n="3">3 Lun et al., F1000Research, 2016</footnote>

<footnote n="4">4 Petropulos et al., Cell, 2016</footnote>
Introduction

* No standard threshold for filtering.

* Same threshold might not be of the *same stringency* in different datasets, notably across species*.

* Can we do better?

<footnote>* See Mansoki et al., Comput Biol Chem., 2016</footnote>
scFeatureFilter

* R package

* Available on GitHub: github.com/gdevailly/scFeatureFilte

* Accepted in Bioconductor

* Might help to set a relevant expression threshold for feature filtering.

library(scFeatureFilter)

Need R ≥ 3.5 (or edit DESCRIPTION to depends: R ≥ 3.4)
## Example datasets:

### 32 scRNA-seq of human embryonic stem cells (Yan et al., Nat Struct Mol Biol, 2013.)

```
dim(scData_hESC)
## [1] 60468 33
```

```
scData_hESC
```

<table><thead><tr><th>gene</th><th>cell_1</th><th>cell_2</th><th>cell_3</th><th>cell_4</th></tr></thead><tbody><tr><td>ENSG00000000003.13</td><td>55.33</td><td>35.98</td><td>53.68</td><td>31.95</td></tr><tr><td>ENSG00000000005.5</td><td>0.00</td><td>0.00</td><td>0.13</td><td>0.00</td></tr><tr><td>ENSG00000000419.11</td><td>53.97</td><td>55.47</td><td>41.87</td><td>110.75</td></tr><tr><td>ENSG00000000457.12</td><td>0.92</td><td>0.22</td><td>0.65</td><td>0.87</td></tr></tbody></table>

Expression matrices can be either `data.frame`, `tibble`, `matrix` or `SingleCellExperiment`.
## Mean-variance exploration:

calculate_cvs(scData_hESC, max zeros = 0.75)[1:4, 1:5]
## # A tibble: 4 x 5
## geneName mean sd cv
## <chr> <dbl> <dbl> <dbl>
## 1 ENSG00000000003.13 70.022328 73.332127 1.0472678
## 2 ENSG00000000005.5 1.000291 2.527535 2.5267987
## 3 ENSG00000000419.11 80.991847 71.534922 0.8832361
## 4 ENSG00000000457.12 1.590995 1.804686 1.1343132
## # ... with 1 more variable: GSM922224_hESCpassage0_Cell4_0 <dbl>

**max zeros:** maximum proportion of 0 value for a feature to be kept
## Mean-variance exploration:

calculate_cvs(scData_hESC, max zeros = 0.75) %>
plot_mean-variance(colourByBin = FALSE) +
annotation_logticks(sides = "l")

<figure><img src="image_3.png" /></figure>
Binning of the genes:

scData_hESC %>
calculate_cvs %>
define_top_genes(window_size = 100) %>
bin_scdata(window_size = 1000) %>
plot_mean-variance() +
annotation_logticks(sides = "l")

<figure><img src="image_4.png" /></figure>
Assumptions

<figure><img src="image_5.png" /></figure>

* high expression = less technical variation

* biological variation = transcription module

* correlation of genes belonging to the same transcription module
## Correlation of the data:

* A reference set of genes: the top bin

* Three control sets of genes: shuffling of the expression values of the reference set

Each gene in each bin is correlated against each gene in the **reference set**, and each gene in the **3 control sets**:

```
corDistrib <- scData_hESC %>%
calculate_cvs %>%
define_top_genes(window_size = 100) %>%
bin_scdata(window_size = 1000) %>%
correlate_windows(n_random = 3)
```

## Mean expression of last top gene: 2114.40221875

## Number of windows: 18
Correlation of the data

corDens <- correlations_to_densities corDistrib)
plot_correlations_distributions corDens, facet_ncol = 6) +
scale_x continuous breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1"))

<figure><img src="image_6.png" /></figure>
Correlation of the data

```
metrics <- get_meanaro median corDistrib)
```

```
plot_correlations_distributions corDens, metrics = metrics, facet_ncol = 6) +
scale_x continuous breaks = c(0, 0.5, 1), labels = c("0", "0.5", "1"))
```

<figure><img src="image_7.png" /></figure>
Threshold decision

plotMetric metrics, show\_ctrl = FALSE, show\_threshold = FALSE)

<figure><img src="image_8.png" /></figure>
Threshold decision

plotMetric metrics, show\_ctrl = TRUE, show\_threshold = FALSE)

<figure><img src="image_9.png" /></figure>
Threshold decision

plotMetric metrics, show\_ctrl = TRUE, show\_threshold = TRUE, threshold = 2

<figure><img src="image_10.png" /></figure>
Getting back the filtered expression matrix:

```
binned_data <- scData_hESC %>%
calculate_cvs %>%
define_top_genes(window_size = 100) %>%
bin_scdata(window_size = 1000)
```

```
determine_bin_cutoff.metrics, threshold = 2)
## [1] 9
```

```
filtered_data <- filter_expression_table(
binned_data,
bin_cutoff = determine_bin_cutoff.metrics)
```

```
nrow(scData_hESC)
## [1] 60468
```

```
nrow(binned_data)
## [1] 17929
```

```
nrow(filter_data)
## [1] 7442
```
A shortcut:

filtered_data <- sc feature_filter(scData_hESC)
## Mean expression of last top gene: 2114.40221875
## Number of windows: 18

dim(scData_hESC)
## [1] 60468 33

dim(filtered_data)
## [1] 7442 33
Testing scFeatureFilter
Testing scFeatureFilter on more datasets:

16 datasets from ConquerDB (human and mouse)

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>
scFeatureFilter on bulk RNA-seq:

48 replicated bulk RNA-seq (yeast) (Gierlinski et al., Bioinformatics, 2015)

<figure><img src="image_15.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_17.png" /></figure>

bin: 13

<figure><img src="image_18.png" /></figure>

bin: 14
Limits of scFeatureFilter
## Lots of parameters?

<table><thead><tr><td>parameter</td><td>description</td><td>default</td></tr></thead><tbody><tr><td>max zeros</td><td>maximum proportion of 0 for a feature to be kept</td><td>0.75</td></tr><tr><td>top window size</td><td>size of the reference set</td><td>100</td></tr><tr><td>other window size</td><td>size of the other bins</td><td>1000</td></tr><tr><td>threshold</td><td>stringency of the selection</td><td>2</td></tr></tbody></table>
# Lots of parameters?

The method is robust to other_windows_size:

<figure><img src="image_19.png" /></figure>
## Robust to max zeros?

high proportion of 0s =~ low expression

* mostly not in the reference set

* more abundant in the low expression bins

* less abundant in the high expression bins

## threshold is a feature

threshold: More or less stringency depending of the use cases and user preference.
top_windows_size

top_windows_size can have massive impact:

* if **too big**: Risk of selecting everything.

* if **too small**: Might not capture enough biological variation.

~100 seems to be a sweet spot on mouse and human data.
top_windows_size

Average auto-correlation of the top window depending on its size:

plot_top_windows_autocor(calculatedcvs(scData_hESC))

<figure><img src="image_20.png" /></figure>
Other limits?

* Tested on a dataset with 1378 cells. Scalable until when?

* Not designed nor tested for 10x genomics scRNA-seq.
Conclusion:

scFeatureFilter: an R package for less arbitrary threshold selection

We are looking for feedback:

• Usefull?

• Overkilled?

• Broken assumptions?

• Better existing methods?
Thanks

Anagha Joshi

Angeles Arzalluz-Luque

Anna Mantsoki

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

ROSLIN

<figure><img src="image_23.png" /></figure>

THE UNIVERSITY
of EDINBURGH

<figure><img src="image_24.png" /></figure>

BBSRC

<figure><img src="image_25.png" /></figure>

CASCADE-FELLOWS