-
Notifications
You must be signed in to change notification settings - Fork 133
/
single_r.Rmd
160 lines (115 loc) · 3 KB
/
single_r.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
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
---
title: "R Notebook"
output: html_notebook
---
```{r}
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("SingleR")
```
```{r}
install.packages('Seurat')
```
```{r}
#very basic Seurat preprocessing
prep_data <- function(Data.path){
Raw_data <- Read10X(data.dir = Data.path)
seuset_data <- CreateSeuratObject(counts = Raw_data, min.cells = 3, min.features = 200)
seuset_data[["percent.mt"]] <- PercentageFeatureSet(seuset_data, pattern = "mt-")
lb <- quantile(seuset_data[["nFeature_RNA"]]$nFeature_RNA, probs = 0.01)
ub <- quantile(seuset_data[["nFeature_RNA"]]$nFeature_RNA, probs = 0.99)
seuset_data <- seuset_data[, seuset_data[["nFeature_RNA"]] > lb & seuset_data[["nFeature_RNA"]] < ub & seuset_data[["percent.mt"]] < 15]
seuset_data <- NormalizeData(object = seuset_data, verbose = FALSE)
seuset_data <- FindVariableFeatures(object = seuset_data, nfeatures = 3000, verbose = FALSE, selection.method = 'vst')
seuset_data <- ScaleData(seuset_data, verbose = FALSE)
seuset_data <- RunPCA(seuset_data, npcs = 20, verbose = FALSE)
seuset_data <- FindNeighbors(seuset_data, dims = 1:20)
seuset_data <- FindClusters(seuset_data, resolution = 0.3)
seuset_data <- RunUMAP(seuset_data, reduction = "pca", dims = 1:20)
return(seuset_data)
}
```
```{r}
data <- prep_data("Lung1/outs/filtered_feature_bc_matrix")
```
```{r}
a <- DimPlot(data, reduction = "umap", label=TRUE)
png("./umap_u.png", res = 250, width = 1500, height = 1500)
print(a)
dev.off()
#a + b + c
a
```
```{r}
library(SingleR)
```
```{r}
ref <- celldex::MouseRNAseqData()
```
```{r}
results <- SingleR(test = as.SingleCellExperiment(data), ref = ref, labels = ref$label.main)
```
```{r}
data$singlr_labels <- results$labels
```
```{r}
DimPlot(data, reduction = 'umap', group.by = 'singlr_labels', label = TRUE)
```
```{r}
FeaturePlot(data, features = c("Ptprc", "Cd3e"))
```
```{r}
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("scRNAseq")
```
```{r}
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("scuttle")
```
```{r}
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("TabulaMurisData")
```
```{r}
library(ExperimentHub)
```
```{r}
eh <- ExperimentHub()
```
```{r}
query(eh, "TabulaMurisData")
```
```{r}
eh[['EH1617']]
```
```{r}
lung_ref <- eh[['EH1617']]
lung_ref <- lung_ref[,lung_ref$tissue == 'Lung']
lung_ref <- lung_ref[,!is.na(lung_ref$cell_ontology_class)]
```
```{r}
lung_ref
```
```{r}
library(scuttle)
```
```{r}
lung_ref <- logNormCounts(lung_ref)
```
```{r}
results <- SingleR(test = as.SingleCellExperiment(data), ref = lung_ref, labels = lung_ref$cell_ontology_class)
```
```{r}
data$singlr_label <- results$labels
```
```{r}
a <- DimPlot(data, reduction = "umap", group.by = 'singlr_label', label = FALSE)
png("./umap_l.png", res = 250, width = 2500, height = 1500)
print(a)
dev.off()
#a + b + c
a
```