forked from satijalab/seurat
-
Notifications
You must be signed in to change notification settings - Fork 0
/
interaction.R
949 lines (932 loc) · 31.1 KB
/
interaction.R
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
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
#' @include seurat.R
NULL
globalVariables(names = 'cell.name', package = 'Seurat', add = TRUE)
#' Merge Seurat Objects
#'
#' Merge two Seurat objects
#'
#' @param object1 First Seurat object to merge
#' @param object2 Second Seurat object to merge
#' @param project Project name (string)
#' @param min.cells Include genes with detected expression in at least this
#' many cells
#' @param min.genes Include cells where at least this many genes are detected
#' @param is.expr Expression threshold for 'detected' gene
#' @param do.normalize Normalize the data after merging. Default is TRUE.
#' If set, will perform the same normalization strategy as stored for the first object
#' @param scale.factor If normalizing on the cell level, this sets the scale factor.
#' @param do.scale In object@@scale.data, perform row-scaling (gene-based
#' z-score). FALSE by default, so run ScaleData after merging.
#' @param do.center In object@@scale.data, perform row-centering (gene-based
#' centering). FALSE by default
#' @param names.field For the initial identity class for each cell, choose this
#' field from the cell's column name
#' @param names.delim For the initial identity class for each cell, choose this
#' delimiter from the cell's column name
#' @param add.cell.id1 String to be appended to the names of all cells in object1
#' @param add.cell.id2 String to be appended to the names of all cells in object2
#'
#' @return Merged Seurat object
#'
#' @import Matrix
#' @importFrom dplyr full_join filter
#'
#' @export
#'
#' @examples
#' # Split pbmc_small for this example
#' pbmc1 <- SubsetData(object = pbmc_small, cells.use = [email protected][1:40])
#' pbmc1
#' pbmc2 <- SubsetData(object = pbmc_small, cells.use = [email protected][41:80])
#' pbmc2
#' # Merge pbmc1 and pbmc2 into one Seurat object
#' pbmc_merged <- MergeSeurat(object1 = pbmc1, object2 = pbmc2)
#' pbmc_merged
#'
MergeSeurat <- function(
object1,
object2,
project = NULL,
min.cells = 0,
min.genes = 0,
is.expr = 0,
do.normalize = TRUE,
scale.factor = 1e4,
do.scale = FALSE,
do.center = FALSE,
names.field = 1,
names.delim = "_",
add.cell.id1 = NULL,
add.cell.id2 = NULL
) {
if (length(x = [email protected]) < 2) {
stop("First object provided has an empty raw.data slot. Adding/Merging performed on raw count data.")
}
if (length(x = [email protected]) < 2) {
stop("Second object provided has an empty raw.data slot. Adding/Merging performed on raw count data.")
}
if (!missing(x = add.cell.id1)) {
[email protected] <- paste(add.cell.id1,[email protected], sep = "_")
colnames(x = [email protected]) <- paste(
add.cell.id1,
colnames(x = [email protected]),
sep = "_"
)
rownames(x = [email protected]) <- paste(
add.cell.id1,
rownames(x = [email protected]),
sep = "_"
)
}
if (!missing(x = add.cell.id2)) {
[email protected] <- paste(add.cell.id2,[email protected], sep = "_")
colnames(x = [email protected]) <- paste(
add.cell.id2,
colnames(x = [email protected]),
sep = "_"
)
rownames(x = [email protected]) <- paste(
add.cell.id2,
rownames(x = [email protected]),
sep = "_"
)
}
if (any([email protected] %in% [email protected])) {
stop("Duplicate cell names, please provide 'add.cell.id1' and/or 'add.cell.id2' for unique names")
}
merged.raw.data <- RowMergeSparseMatrices(
mat1 = [email protected][,[email protected]],
mat2 = [email protected][,[email protected]]
)
project <- SetIfNull(x = project, default = [email protected])
[email protected]$cell.name <- rownames(x = [email protected])
[email protected]$cell.name <- rownames(x = [email protected])
merged.meta.data <- suppressMessages(
suppressWarnings(
full_join(x = [email protected], y = [email protected])
)
)
merged.object <- CreateSeuratObject(
raw.data = merged.raw.data,
project = project,
min.cells = min.cells,
min.genes = min.genes,
is.expr = is.expr,
normalization.method = NULL,
scale.factor = scale.factor,
do.scale = FALSE,
do.center = FALSE,
names.field = names.field,
names.delim = names.delim
)
if (do.normalize) {
normalization.method.use = GetCalcParam(
object = object1,
calculation = "NormalizeData",
parameter = "normalization.method"
)
scale.factor.use = GetCalcParam(
object = object1,
calculation = "NormalizeData",
parameter = "scale.factor"
)
if (is.null(normalization.method.use)) {
normalization.method.use="LogNormalize"
scale.factor.use=10000
}
merged.object <- NormalizeData(
object = merged.object,
assay.type = "RNA",
normalization.method=normalization.method.use,
scale.factor=scale.factor.use
)
}
if (do.scale | do.center) {
merged.object <- ScaleData(
object = merged.object,
do.scale = do.scale,
do.center = do.center
)
}
merged.meta.data %>% filter(
cell.name %in% [email protected]
) -> merged.meta.data
rownames(x= merged.meta.data) <- [email protected]
merged.meta.data$cell.name <- NULL
[email protected] <- merged.meta.data
return(merged.object)
}
#' Add samples into existing Seurat object.
#'
#' @param object Seurat object
#' @param project Project name (string)
#' @param new.data Data matrix for samples to be added
#' @param min.cells Include genes with detected expression in at least this
#' many cells
#' @param min.genes Include cells where at least this many genes are detected
#' @param is.expr Expression threshold for 'detected' gene
#' @param do.normalize Normalize the data after merging. Default is TRUE.
#' If set, will perform the same normalization strategy as stored in the object
#' @param scale.factor scale factor in the log normalization
#' @param do.scale In object@@scale.data, perform row-scaling (gene-based z-score)
#' @param do.center In object@@scale.data, perform row-centering (gene-based
#' centering)
#' @param names.field For the initial identity class for each cell, choose this
#' field from the cell's column name
#' @param names.delim For the initial identity class for each cell, choose this
#' delimiter from the cell's column name
#' @param meta.data Additional metadata to add to the Seurat object. Should be
#' a data frame where the rows are cell names, and the columns are additional
#' metadata fields
#' @param add.cell.id String to be appended to the names of all cells in new.data. E.g. if add.cell.id = "rep1",
#' "cell1" becomes "cell1.rep1"
#'
#' @import Matrix
#' @importFrom dplyr full_join
#'
#' @export
#'
#' @examples
#' pbmc1 <- SubsetData(object = pbmc_small, cells.use = [email protected][1:40])
#' pbmc1
#' pbmc2 <- SubsetData(object = pbmc_small, cells.use = [email protected][41:80])
#' pbmc2_data <- pbmc2@data
#' dim(pbmc2_data)
#' pbmc_added <- AddSamples(object = pbmc1, new.data = pbmc2_data)
#' pbmc_added
#'
AddSamples <- function(
object,
new.data,
project = NULL,
min.cells = 0,
min.genes = 0,
is.expr = 0,
do.normalize = TRUE,
scale.factor = 1e4,
do.scale = FALSE,
do.center = FALSE,
names.field = 1,
names.delim = "_",
meta.data = NULL,
add.cell.id = NULL
) {
if (length(x = [email protected]) < 2) {
stop("Object provided has an empty raw.data slot. Adding/Merging performed on raw count data.")
}
if (!missing(x = add.cell.id)) {
colnames(x = new.data) <- paste(
add.cell.id,
colnames(x = new.data),
sep = "_"
)
}
if (any(colnames(x = new.data) %in% [email protected])) {
stop("Duplicate cell names, please provide 'add.cell.id' for unique names")
}
combined.data <- RowMergeSparseMatrices(
mat1 = [email protected][, [email protected]],
mat2 = new.data
)
if (is.null(x = meta.data)) {
filler <- matrix(NA, nrow = ncol(new.data), ncol = ncol([email protected]))
rownames(filler) <- colnames(new.data)
colnames(filler) <- colnames([email protected])
filler <- as.data.frame(filler)
combined.meta.data <- rbind([email protected], filler)
} else {
combined.meta.data <- suppressMessages(
suppressWarnings(
full_join(x = [email protected], y = meta.data)
)
)
}
combined.meta.data$nGene <- NULL
combined.meta.data$nUMI <- NULL
if (!is.null(x = add.cell.id)) {
combined.meta.data$orig.ident <- factor(
x = combined.meta.data$orig.ident,
levels = c(levels(x = combined.meta.data$orig.ident), add.cell.id)
)
combined.meta.data[colnames(new.data), ] <- add.cell.id
}
project <- SetIfNull(x = project, default = [email protected])
new.object <- CreateSeuratObject(
raw.data = combined.data,
project = project,
min.cells = min.cells,
min.genes = min.genes,
is.expr = is.expr,
scale.factor = scale.factor,
do.scale = F,
do.center = F,
names.field = names.field,
names.delim = names.delim
)
if (do.normalize) {
normalization.method.use = GetCalcParam(
object = object,
calculation = "NormalizeData",
parameter = "normalization.method"
)
scale.factor.use = GetCalcParam(
object = object,
calculation = "NormalizeData",
parameter = "scale.factor"
)
if (is.null(x = normalization.method.use)) {
normalization.method.use <- "LogNormalize"
scale.factor.use <- 10000
}
new.object <- NormalizeData(
object = new.object,
assay.type = "RNA",
normalization.method = normalization.method.use,
scale.factor = scale.factor.use
)
}
if (do.scale | do.center) {
new.object <- ScaleData(
object = new.object,
do.scale = do.scale,
do.center = do.center
)
}
[email protected]$orig.ident <- NULL
[email protected] <- cbind([email protected], combined.meta.data)
return(new.object)
}
#' Return a subset of the Seurat object
#'
#' Creates a Seurat object containing only a subset of the cells in the
#' original object. Takes either a list of cells to use as a subset, or a
#' parameter (for example, a gene), to subset on.
#'
#' @param object Seurat object
#' @param cells.use A vector of cell names to use as a subset. If NULL
#' (default), then this list will be computed based on the next three
#' arguments. Otherwise, will return an object consissting only of these cells
#' @param subset.name Parameter to subset on. Eg, the name of a gene, PC1, a
#' column name in object@@meta.data, etc. Any argument that can be retreived
#' using FetchData
#' @param ident.use Create a cell subset based on the provided identity classes
#' @param ident.remove Subtract out cells from these identity classes (used for filtration)
#' @param accept.low Low cutoff for the parameter (default is -Inf)
#' @param accept.high High cutoff for the parameter (default is Inf)
#' @param do.center Recenter the new object@@scale.data
#' @param do.scale Rescale the new object@@scale.data. FALSE by default
#' @param max.cells.per.ident Can be used to downsample the data to a certain max per cell ident. Default is inf.
#' @param random.seed Random seed for downsampling
#' @param do.clean Only keep object@@raw.data and object@@data. Cleans out most
#' other slots. Can be useful if you want to start a fresh analysis on just a
#' subset of the data. Also clears out stored clustering results in
#' object@@meta.data (any columns containing "res"). Will by default subset the
#' raw.data slot.
#' @param subset.raw Also subset object@@raw.data
#' @param \dots Additional arguments to be passed to FetchData (for example,
#' use.imputed=TRUE)
#'
#' @return Returns a Seurat object containing only the relevant subset of cells
#'
#' @importFrom stats complete.cases
#'
#' @export
#'
#' @examples
#' pbmc1 <- SubsetData(object = pbmc_small, cells.use = [email protected][1:40])
#' pbmc1
#'
SubsetData <- function(
object,
cells.use = NULL,
subset.name = NULL,
ident.use = NULL,
ident.remove = NULL,
accept.low = -Inf,
accept.high = Inf,
do.center = FALSE,
do.scale = FALSE,
max.cells.per.ident = Inf,
random.seed = 1,
do.clean = FALSE,
subset.raw,
...
) {
data.use <- NULL
cells.use <- SetIfNull(x = cells.use, default = [email protected])
if (!is.null(x = ident.use)) {
ident.use <- setdiff(ident.use, ident.remove)
cells.use <- WhichCells(object, ident.use)
}
if ((is.null(x = ident.use)) && ! is.null(x = ident.remove)) {
ident.use <- setdiff(unique(object@ident), ident.remove)
cells.use <- WhichCells(object, ident.use)
}
if (! is.null(x = subset.name)) {
data.use <- FetchData(object, subset.name, ...)
if (length(x = data.use) == 0) {
return(object)
}
subset.data <- data.use[, subset.name]
pass.inds <- which(x = (subset.data > accept.low) & (subset.data < accept.high))
cells.use <- rownames(data.use)[pass.inds]
}
cells.use <- intersect(x = cells.use, y = [email protected])
cells.use <- WhichCells(
object = object,
cells.use = cells.use,
max.cells.per.ident = max.cells.per.ident,
random.seed = random.seed
)
[email protected] <- cells.use
object@data <- object@data[, cells.use]
if(! is.null(x = [email protected])) {
if (length(x = colnames(x = [email protected]) > 0)) {
[email protected][, cells.use]
complete.cases([email protected]), # Row
cells.use # Columns
]
}
}
if (do.scale) {
object <- ScaleData(
object = object,
do.scale = do.scale,
do.center = do.center
)
complete.cases([email protected]), # Row
cells.use # Column
]
}
object@ident <- drop.levels(x = object@ident[cells.use])
if (length(x = object@dr) > 0) {
for (i in 1:length(object@dr)) {
if(length(object@dr[[i]]@cell.embeddings) > 0){
object@dr[[i]]@cell.embeddings <- object@dr[[i]]@cell.embeddings[cells.use, ,drop = FALSE]
}
}
}
# handle multimodal casess
if (! .hasSlot(object = object, name = "assay")) {
object@assay <- list()
}
if (length(object@assay) > 0) {
for(i in 1:length(object@assay)) {
if ((! is.null(x = object@assay[[i]]@raw.data)) && (ncol(x = object@assay[[i]]@raw.data) > 1)) {
object@assay[[i]]@raw.data <- object@assay[[i]]@raw.data[, cells.use]
}
if ((! is.null(x = object@assay[[i]]@data)) && (ncol(x = object@assay[[i]]@data) > 1)) {
object@assay[[i]]@data <- object@assay[[i]]@data[, cells.use]
}
if ((! is.null(x = object@assay[[i]]@scale.data)) && (ncol(x = object@assay[[i]]@scale.data) > 1)) {
object@assay[[i]]@scale.data <- object@assay[[i]]@scale.data[, cells.use]
}
}
}
#[email protected][email protected][cells.use, ]
# [email protected] <- data.frame([email protected][cells.use,])
# colnames(x = [email protected])[1] <- "nGene"
# rownames(x = [email protected]) <- colnames(x = object@data)
[email protected] <- data.frame([email protected][cells.use,])
#[email protected]=data.frame([email protected][cells.use,]); colnames([email protected])[1]="nGene"; rownames([email protected])=colnames(object@data)
if (do.clean){
calcs.to.keep <- c("CreateSeuratObject", "NormalizeData", "ScaleData")
[email protected] <- [email protected][calcs.to.keep]
[email protected] <- vector()
[email protected] <- data.frame()
[email protected] <- list()
object@snn <- as(matrix(), 'dgCMatrix')
[email protected] <- matrix()
object@misc <- NULL
object@kmeans <- NULL
object@dr <- list()
if(missing(subset.raw)) {
subset.raw <- TRUE
}
[email protected][, sapply(colnames([email protected]), function(x){grepl("res", x)})] <- NULL
}
if(!missing(subset.raw)){
if(subset.raw){
[email protected] <- [email protected][, cells.use]
}
}
return(object)
}
#' Reorder identity classes
#'
#' Re-assigns the identity classes according to the average expression of a particular feature (i.e, gene expression, or PC score)
#' Very useful after clustering, to re-order cells, for example, based on PC scores
#'
#' @param object Seurat object
#' @param feature Feature to reorder on. Default is PC1
#' @param rev Reverse ordering (default is FALSE)
#' @param aggregate.fxn Function to evaluate each identity class based on (default is mean)
#' @param reorder.numeric Rename all identity classes to be increasing numbers starting from 1 (default is FALSE)
#' @param \dots additional arguemnts (i.e. use.imputed=TRUE)
#'
#' @return A seurat object where the identity have been re-oredered based on the average.
#'
#' @export
#'
#' @examples
#' head(x = pbmc_small@ident)
#' pbmc_small <- ReorderIdent(object = pbmc_small)
#' head(x = pbmc_small@ident)
#'
ReorderIdent <- function(
object,
feature = "PC1",
rev = FALSE,
aggregate.fxn = mean,
reorder.numeric = FALSE,
...
) {
ident.use <- object@ident
data.use <- FetchData(object = object, vars.all = feature, ...)[, 1]
revFxn <- Same
if (rev) {
revFxn <- function(x) {
return(max(x) + 1 - x)
}
}
names.sort <- names(
x = revFxn(
sort(
x = tapply(
X = data.use,
INDEX = (ident.use),
FUN = aggregate.fxn
)
)
)
)
ident.new <- factor(x = ident.use, levels = names.sort, ordered = TRUE)
if (reorder.numeric) {
ident.new <- factor(
x = revFxn(
rank(
tapply(
X = data.use,
INDEX = as.numeric(x = ident.new),
FUN = mean
)
)
)[as.numeric(ident.new)],
levels = 1:length(x = levels(x = ident.new)),
ordered = TRUE
)
}
names(x = ident.new) <- names(x = ident.use)
object@ident <- ident.new
return(object)
}
#' Access cellular data
#'
#' Retreives data (gene expression, PCA scores, etc, metrics, etc.) for a set
#' of cells in a Seurat object
#'
#' @param object Seurat object
#' @param vars.all List of all variables to fetch
#' @param cells.use Cells to collect data for (default is all cells)
#' @param use.imputed For gene expression, use imputed values. Default is FALSE
#' @param use.scaled For gene expression, use scaled values. Default is FALSE
#' @param use.raw For gene expression, use raw values. Default is FALSE
#'
#' @return A data frame with cells as rows and cellular data as columns
#'
#' @export
#'
#' @examples
#' pc1 <- FetchData(object = pbmc_small, vars.all = 'PC1')
#' head(x = pc1)
#'
FetchData <- function(
object,
vars.all = NULL,
cells.use = NULL,
use.imputed = FALSE,
use.scaled = FALSE,
use.raw = FALSE
) {
cells.use <- SetIfNull(x = cells.use, default = [email protected])
data.return <- data.frame(row.names = cells.use)
data.expression <- as.matrix(x = data.frame(row.names = cells.use))
if (length(which(c(use.imputed, use.scaled, use.raw))) > 1) {
stop("Can only set one of the following to TRUE: use.imputed, use.scaled, use.raw")
}
slot.use <- "data"
# if any vars passed are genes, subset expression data
gene.check <- vars.all %in% rownames(object@data)
if (use.scaled) {
slot.use <- "scale.data"
gene.check <- vars.all %in% rownames([email protected])
}
if (use.raw) {
slot.use <- "raw.data"
}
if (any(gene.check)) {
if (use.imputed) {
gene.check <- vars.all %in% rownames(object@imputed)
if (length(object@imputed) == 0) {
stop ("Imputed expression values not calculated yet.")
}
data.expression <- t(object@imputed[vars.all[gene.check], cells.use, drop = FALSE])
} else {
data.expression <- GetAssayData(object, assay.type = "RNA", slot = slot.use)
data.expression <- t(data.expression[vars.all[gene.check], cells.use, drop = FALSE])
}
if (all(gene.check)) {
return(as.matrix(x = data.expression))
}
}
# now check for multimodal data
if (length(x = object@assay) > 0) {
data.types <- names(x = object@assay)
for (data.type in data.types) {
all_data <- (GetAssayData(
object = object,
assay.type = data.type,
slot = slot.use
))
genes.include <- intersect(x = vars.all, y = rownames(x = all_data))
data.expression <- cbind(
data.expression,
t(x = all_data[genes.include, , drop = FALSE])
)
}
}
var.options <- c("meta.data", "mix.probs", "gene.scores")
if (length(x = names(x = object@dr)) > 0) {
dr.options <- names(x = object@dr)
dr.names <- paste0("dr$", names(x = object@dr), "@key")
dr.names <- sapply(
X = dr.names,
FUN = function(x) {
return(eval(expr = parse(text = paste0("object@", x))))
}
)
names(x = dr.names) <- dr.options
var.options <- c(var.options, dr.names)
}
[email protected][,"ident"] <- object@ident[rownames(x = [email protected])]
for (my.var in vars.all) {
data.use=data.frame()
if (my.var %in% colnames(data.expression)) {
data.use <- data.expression
} else {
for(i in var.options) {
if (all(unlist(x = strsplit(x = my.var, split = "[0-9]+")) == i)) {
eval(
expr = parse(
text = paste0(
"data.use <- object@dr$",
names(x = var.options[which(i == var.options)]),
"@cell.embeddings"
)
)
)
colnames(x = data.use) <- paste0(i, 1:ncol(x = data.use))
break
}
}
}
if (my.var %in% colnames([email protected])) {
data.use <- [email protected][, my.var, drop = FALSE]
}
if (ncol(x = data.use) == 0) {
stop(paste("Error:", my.var, "not found"))
}
cells.use <- intersect(x = cells.use, y = rownames(x = data.use))
if (! my.var %in% colnames(x = data.use)) {
stop(paste("Error:", my.var, "not found"))
}
data.add <- data.use[cells.use, my.var]
if (is.null(x = data.add)) {
stop(paste("Error:", my.var, "not found"))
}
data.return <- cbind(data.return, data.add)
}
colnames(x = data.return) <- vars.all
rownames(x = data.return) <- cells.use
return(data.return)
}
#' FastWhichCells
#' Identify cells matching certain criteria (limited to character values)
#' @param object Seurat object
#' @param group.by Group cells in different ways (for example, orig.ident). Should be a column name in [email protected]
#' @param subset.value Return cells matching this value
#' @param invert invert cells to return.FALSE by default
#'
#' @export
#'
#' @examples
#' FastWhichCells(object = pbmc_small, group.by = 'res.1', subset.value = 1)
#'
FastWhichCells <- function(object, group.by, subset.value, invert = FALSE) {
object <- SetAllIdent(object = object, id = group.by)
cells.return <- WhichCells(object = object, ident = subset.value)
if (invert) {
cells.return <- setdiff(x = [email protected], y = cells.return)
}
return(cells.return)
}
#' Identify cells matching certain criteria
#'
#' Returns a list of cells that match a particular set of criteria such as
#' identity class, high/low values for particular PCs, ect..
#'
#' @param object Seurat object
#' @param ident Identity classes to subset. Default is all identities.
#' @param ident.remove Indentity classes to remove. Default is NULL.
#' @param cells.use Subset of cell names
#' @param subset.name Parameter to subset on. Eg, the name of a gene, PC1, a
#' column name in object@@meta.data, etc. Any argument that can be retreived
#' using FetchData
#' @param accept.low Low cutoff for the parameter (default is -Inf)
#' @param accept.high High cutoff for the parameter (default is Inf)
#' @param accept.value Returns all cells with the subset name equal to this value
#' @param max.cells.per.ident Can be used to downsample the data to a certain max per cell ident. Default is inf.
#' @param random.seed Random seed for downsampling
#'
#' @return A vector of cell names
#'
#' @export
#'
#' @examples
#' WhichCells(object = pbmc_small, ident = 2)
#'
WhichCells <- function(
object,
ident = NULL,
ident.remove = NULL,
cells.use = NULL,
subset.name = NULL,
accept.low = -Inf,
accept.high = Inf,
accept.value = NULL,
max.cells.per.ident = Inf,
random.seed = 1
) {
set.seed(seed = random.seed)
cells.use <- SetIfNull(x = cells.use, default = [email protected])
ident <- SetIfNull(x = ident, default = unique(x = object@ident))
ident <- setdiff(x = ident, y = ident.remove)
if (! all(ident %in% unique(x = object@ident))) {
bad.idents <- ident[! (ident %in% unique(x = object@ident))]
stop(paste("Identity :", bad.idents, "not found. "))
}
cells.to.use <- character()
for (id in ident) {
cells.in.ident <- object@ident[cells.use]
cells.in.ident <- names(x = cells.in.ident[cells.in.ident == id])
cells.in.ident <- cells.in.ident[! is.na(x = cells.in.ident)]
if (length(x = cells.in.ident) > max.cells.per.ident) {
cells.in.ident <- sample(x = cells.in.ident, size = max.cells.per.ident)
}
cells.to.use <- c(cells.to.use, cells.in.ident)
}
cells.use <- cells.to.use
if (! is.null(x = subset.name)){
subset.name <- as.character(subset.name)
data.use <- FetchData(
object = object,
vars.all = subset.name,
cells.use = cells.use
)
if (length(x = data.use) == 0) {
stop(paste("Error : ", id, " not found"))
}
subset.data <- data.use[, subset.name, drop = F]
if(! is.null(x = accept.value)) {
pass.inds <- which(x = subset.data == accept.value)
} else {
pass.inds <- which(x = (subset.data > accept.low) & (subset.data < accept.high))
}
cells.use <- rownames(x = data.use)[pass.inds]
}
return(cells.use)
}
#' Switch identity class definition to another variable
#'
#' @param object Seurat object
#' @param id Variable to switch identity class to (for example, 'DBclust.ident', the output
#' of density clustering) Default is orig.ident - the original annotation pulled from the cell name.
#'
#' @return A Seurat object where object@@ident has been appropriately modified
#'
#' @export
#'
#' @examples
#' head(x = pbmc_small@ident)
#' pbmc_small <- SetAllIdent(object = pbmc_small, id = 'orig.ident')
#' head(x = pbmc_small@ident)
#'
SetAllIdent <- function(object, id = NULL) {
id <- SetIfNull(x = id, default = "orig.ident")
if (id %in% colnames(x = [email protected])) {
cells.use <- rownames(x = [email protected])
ident.use <- [email protected][, id]
object <- SetIdent(
object = object,
cells.use = cells.use,
ident.use = ident.use
)
}
return(object)
}
#' Rename one identity class to another
#'
#' Can also be used to join identity classes together (for example, to merge clusters).
#'
#' @param object Seurat object
#' @param old.ident.name The old identity class (to be renamed)
#' @param new.ident.name The new name to apply
#'
#' @return A Seurat object where object@@ident has been appropriately modified
#'
#' @export
#'
#' @examples
#' head(x = pbmc_small@ident)
#' pbmc_small <- RenameIdent(
#' object = pbmc_small,
#' old.ident.name = 0,
#' new.ident.name = 'cluster_0'
#' )
#' head(x = pbmc_small@ident)
#'
RenameIdent <- function(object, old.ident.name = NULL, new.ident.name = NULL) {
if (! old.ident.name %in% object@ident) {
stop(paste("Error : ", old.ident.name, " is not a current identity class"))
}
new.levels <- old.levels <- levels(x = object@ident)
# new.levels <- old.levels
if (new.ident.name %in% old.levels) {
new.levels <- new.levels[new.levels != old.ident.name]
}
if(! (new.ident.name %in% old.levels)) {
new.levels[new.levels == old.ident.name] <- new.ident.name
}
ident.vector <- as.character(x = object@ident)
names(x = ident.vector) <- names(object@ident)
ident.vector[WhichCells(object = object, ident = old.ident.name)] <- new.ident.name
object@ident <- factor(x = ident.vector, levels = new.levels)
return(object)
}
#' Set identity class information
#'
#' Stashes the identity in data.info to be retrieved later. Useful if, for example, testing multiple clustering parameters
#'
#' @param object Seurat object
#' @param save.name Store current object@@ident under this column name in object@@meta.data. Can be easily retrived with SetAllIdent
#'
#' @return A Seurat object where object@@ident has been appropriately modified
#'
#' @export
#'
#' @examples
#' head(x = [email protected])
#' pbmc_small <- StashIdent(object = pbmc_small, save.name = 'cluster.ident')
#' head(x = [email protected])
#'
StashIdent <- function(object, save.name = "oldIdent") {
[email protected][, save.name] <- as.character(x = object@ident)
return(object)
}
#' Set identity class information
#'
#' Sets the identity class value for a subset (or all) cells
#'
#' @param object Seurat object
#' @param cells.use Vector of cells to set identity class info for (default is
#' all cells)
#' @param ident.use Vector of identity class values to assign (character
#' vector)
#'
#' @return A Seurat object where object@@ident has been appropriately modified
#'
#' @importFrom gdata drop.levels
#'
#' @export
#'
#' @examples
#' cluster2 <- WhichCells(object = pbmc_small, ident = 2)
#' pbmc_small@ident[cluster2]
#' pbmc_small <- SetIdent(
#' object = pbmc_small,
#' cells.use = cluster2,
#' ident.use = 'cluster_2'
#' )
#' pbmc_small@ident[cluster2]
#'
SetIdent <- function(object, cells.use = NULL, ident.use = NULL) {
cells.use <- SetIfNull(x = cells.use, default = [email protected])
if (length(x = setdiff(x = cells.use, y = [email protected]) > 0)) {
stop(paste(
"ERROR : Cannot find cells ",
setdiff(x = cells.use, y = [email protected])
))
}
ident.new <- setdiff(x = ident.use, y = levels(x = object@ident))
object@ident <- factor(
x = object@ident,
levels = unique(
x = c(
as.character(x = object@ident),
as.character(x = ident.new)
)
)
)
object@ident[cells.use] <- ident.use
object@ident <- drop.levels(x = object@ident)
return(object)
}
#' Add Metadata
#'
#' Adds additional data for single cells to the Seurat object. Can be any piece
#' of information associated with a cell (examples include read depth,
#' alignment rate, experimental batch, or subpopulation identity). The
#' advantage of adding it to the Seurat object is so that it can be
#' analyzed/visualized using FetchData, VlnPlot, GenePlot, SubsetData, etc.
#'
#' @param object Seurat object
#' @param metadata Data frame where the row names are cell names (note : these
#' must correspond exactly to the items in object@@cell.names), and the columns
#' are additional metadata items.
#' @param col.name Name for metadata if passing in single vector of information
#'
#' @return Seurat object where the additional metadata has been added as
#' columns in object@@meta.data
#'
#' @export
#'
#' @examples
#' cluster_letters <- LETTERS[pbmc_small@ident]
#' pbmc_small <- AddMetaData(
#' object = pbmc_small,
#' metadata = cluster_letters,
#' col.name = 'letter.idents'
#' )
#' head(x = [email protected])
#'
AddMetaData <- function(object, metadata, col.name = NULL) {
if (typeof(x = metadata) != "list") {
metadata <- as.data.frame(x = metadata)
if (is.null(x = col.name)) {
stop("Please provide a name for provided metadata")
}
colnames(x = metadata) <- col.name
}
cols.add <- colnames(x = metadata)
meta.add <- metadata[rownames(x = [email protected]), cols.add]
if (all(is.null(x = meta.add))) {
stop("Metadata provided doesn't match the cells in this object")
}
[email protected][, cols.add] <- meta.add
return(object)
}