forked from hemberg-lab/scRNA.seq.course
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path22-imputation.Rmd
125 lines (103 loc) · 3.67 KB
/
22-imputation.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
---
output: html_document
---
## Imputation
```{r, echo=FALSE}
library(knitr)
opts_chunk$set(fig.align="center", echo=FALSE)
```
```{r, echo=TRUE, message=FALSE, warning=FALSE}
library(scImpute)
library(SC3)
library(scater)
library(SingleCellExperiment)
set.seed(1234567)
```
As discussed previously, one of the main challenges when analyzing scRNA-seq data is the presence of zeros, or dropouts. The dropouts are assumed to have arisen for three possible reasons:
* The gene was not expressed in the cell and hence there are no transcripts to sequence
* The gene was expressed, but for some reason the transcripts were lost somewhere prior to sequencing
* The gene was expressed and transcripts were captured and turned into cDNA, but the sequencing depth was not sufficient to produce any reads.
Thus, dropouts could be result of experimental shortcomings, and if this is the case then we would like to provide computational corrections. One possible solution is to impute the dropouts in the expression matrix. To be able to impute gene expression values, one must have an underlying model. However, since we do not know which dropout events are technical artefacts and which correspond to the transcript being truly absent, imputation is a difficult challenges.
To the best of our knowledge, there are currently two different imputation methods available: MAGIC [@Van_Dijk2017-bh] and scImpute [@Li2017-tz]. Since [MAGIC](https://github.com/pkathail/magic) is only available for Python or Matlab, we will only cover [scImpute](https://github.com/Vivianstats/scImpute) in this chapter.
### scImpute
To test `scImpute`, we use the default parameters and we apply it to the Deng dataset that we have worked with before. scImpute takes a .csv or .txt file as an input:
```{r}
deng <- readRDS("deng/deng-reads.rds")
write.csv(counts(deng), "deng.csv")
scimpute(
count_path = "deng.csv",
infile = "csv",
outfile = "txt",
out_dir = "./",
Kcluster = 10,
ncores = 2
)
```
Now we can compare the results with original data by considering a PCA plot
```{r}
res <- read.table("scimpute_count.txt")
colnames(res) <- NULL
res <- SingleCellExperiment(
assays = list(logcounts = log2(as.matrix(res) + 1)),
colData = colData(deng)
)
plotPCA(
res,
colour_by = "cell_type2"
)
```
Compare this result to the original data in Fig. X, Chapter Y. What are the most significant differences?
To evaluate the impact of the imputation, we use `SC3` to cluster the imputed matrix
```{r}
res <- sc3_estimate_k(res)
metadata(res)$sc3$k_estimation
res <- sc3(res, ks = 10, gene_filter = FALSE)
sc3_plot_consensus(
res,
k = 10,
show_pdata = "cell_type2"
)
plotPCA(
res,
colour_by = "sc3_10_clusters"
)
```
Based on the PCA and the clustering results, do you think that imputation is a good idea for the Pollen dataset?
### MAGIC
```{r}
system("MAGIC.py -d deng.csv -o MAGIC_count.csv --cell-axis columns -l 1 --pca-non-random csv")
```
```{r}
res <- read.csv("MAGIC_count.csv", header = TRUE)
rownames(res) <- res[,1]
res <- res[,-1]
res <- t(res)
res <- SingleCellExperiment(
assays = list(logcounts = res),
colData = colData(deng)
)
plotPCA(
res,
colour_by = "cell_type2"
)
```
Compare this result to the original data in Fig. X, Chapter Y. What are the most significant differences?
To evaluate the impact of the imputation, we use `SC3` to cluster the imputed matrix
```{r, eval = FALSE}
res <- sc3_estimate_k(res)
metadata(res)$sc3$k_estimation
res <- sc3(res, ks = 10, gene_filter = FALSE)
sc3_plot_consensus(
res,
k = 10,
show_pdata = "cell_type2"
)
plotPCA(
res,
colour_by = "sc3_10_clusters"
)
```
### sessionInfo()
```{r echo=FALSE}
sessionInfo()
```