-
Notifications
You must be signed in to change notification settings - Fork 0
/
SUV39H1_CAR_T_integrated_premerged_CD4s.R
3293 lines (2932 loc) · 235 KB
/
SUV39H1_CAR_T_integrated_premerged_CD4s.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
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
library(Seurat)
library(ggplot2)
library(sctransform)
library(patchwork)
library(pheatmap)
library(Scillus)
library(pals)
library(RColorBrewer)
library(tidyverse)
library(paletteer)
library(EnhancedVolcano)
library(readxl)
library(ggsignif)
library(harmony)
library(Nebulosa)
library(data.table)
library(gridExtra)
library(CytoTRACE)
library(future)
plan("multiprocess", workers = 16, strategy = "multicore")
options(future.globals.maxSize = 1500 * 2024^2)
# Signatures
sign_term_trem2 <- toupper(readLines("E:/Work/Signatures/Resident Memory/Synder_et_al_2019/sign_Term_Trem_Synder_bis.txt"))
sign_memory <- toupper(read.csv("E:/Work/Signatures/Stemness/signature_memory_stemness.csv", h = F)[,1])
sign_senescence <- toupper(read.csv("E:/Work/Signatures/Senescence/sign_sag.txt", h = F)[,1])
sign_apoptosis <- toupper(read.csv("E:/Work/Signatures/Signatures Papier Leticia/Apoptosis.txt", h = F)[,1])
sign_DNA_repair <- toupper(read.csv("E:/Work/Signatures/DNA repair/sign_DNA_repair.txt", h = F)[,1])
sign_core_DNA_repair <- toupper(read.csv("E:/Work/Signatures/DNA repair/Core_DNA_damage_repair.txt", h = F)[,1])
sign_DNA_damage_repair <- toupper(read.csv("E:/Work/Signatures/DNA repair/DNA_damage_repair_Knijnenburg.txt", h = F)[,1])
# Colors setup
mycolors <- pals::kelly(18)
mycolors.CD4 <- mycolors[5:18]
# Saving
save(suv.car.t.integrated.merged.timepoints.premerged.cd4s, file = "E:/Work/A_R_files/SUV39H1_CAR_T_integrated_premerged_CD4s.Rdata")
save(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, file = "E:/Work/A_R_files/SUV39H1_CAR_T_integrated_premerged_CD4s_cleaned.Rdata")
# Merging CD4s only
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- subset(suv.car.t.integrated.merged.timepoints.premerged, subset = cd4A > 0 & cd4B > 0, invert = T)
TCR.removed <- rownames(suv.car.t.integrated.merged.timepoints.premerged.cd4s)[which(rownames(suv.car.t.integrated.merged.timepoints.premerged.cd4s) %!in% TCR.features)]
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- RunPCA(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, npcs = 50, verbose = FALSE, features = TCR.removed)
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- RunHarmony(suv.car.t.integrated.merged.timepoints.premerged.cd4s, group.by.vars = c("orig.ident", "timepoint"), project.dim = T)
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- RunUMAP(suv.car.t.integrated.merged.timepoints.premerged.cd4s, reduction = "harmony", dims = 1:50)
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- FindNeighbors(suv.car.t.integrated.merged.timepoints.premerged.cd4s, reduction = "harmony", dims = 1:50)
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- FindClusters(suv.car.t.integrated.merged.timepoints.premerged.cd4s, graph.name = "RNA_snn", resolution = 0.8)
DimPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, reduction = "umap", label = F, pt.size = 1, label.size = 4, cols = mycolors.CD4) + NoLegend() +
DimPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, reduction = "umap", group.by = "genotype", label = F, pt.size = 1, cols = c("black", "red"), label.size = 6) +
DimPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, reduction = "umap", group.by = "timepoint", label = F, pt.size = 1, label.size = 6) + NoLegend()
# Densities with cells as grey points in the background
umap <- as.data.frame(Embeddings(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, reduction = "umap"))
ggplot() +
geom_point(data = umap, aes(x = UMAP_1, y = UMAP_2), color = "grey")
# Densities by genotype and timepoint
df <- as.data.frame(Embeddings(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, reduction = "umap"))
df <- cbind(df, suv.car.t.integrated.merged.timepoints.premerged.cd4s$timepoint, suv.car.t.integrated.merged.timepoints.premerged.cd4s$genotype)
colnames(df) <- c("UMAP_1", "UMAP_2", "timepoint", "genotype")
df$timepoint <- factor(x = df$timepoint, levels = c("D8", "D28"))
df$genotype <- factor(x = df$genotype, levels = c("Mock", "KO"))
ggplot(data = df) +
geom_point(data = df, aes(x = UMAP_1, y = UMAP_2), color = "lightgrey", size = 0.5) +
stat_density_2d(aes(x = UMAP_1, y = UMAP_2, fill = after_stat(level), alpha = (..level..)), bins = 40, geom = "polygon") +
theme_classic() +
xlim(c(-8, 10)) +
ylim(c(-7, 8)) +
scale_fill_viridis_c(option = "plasma") +
facet_grid(. ~ timepoint + genotype) +
theme(
strip.text.x = element_text(
size = 20, color = "black",
face = "bold"
),
strip.text.y = element_text(
size = 20, color = "black",
face = "bold"
)
)
# Densities by genotype and timepoint
df <- as.data.frame(Embeddings(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, reduction = "umap"))
df <- cbind(df, suv.car.t.integrated.merged.timepoints.premerged.cd4s$timepoint, suv.car.t.integrated.merged.timepoints.premerged.cd4s$genotype)
colnames(df) <- c("UMAP_1", "UMAP_2", "timepoint", "genotype")
df$timepoint <- factor(x = df$timepoint, levels = c("D8", "D28"))
df$genotype <- factor(x = df$genotype, levels = c("Mock", "gSUV"))
ggplot(data = df) +
geom_point(data = df, aes(x = UMAP_1, y = UMAP_2), color = "lightgrey", size = 0.5) +
stat_density_2d(aes(x = UMAP_1, y = UMAP_2, fill = after_stat(level), alpha = (..level..)), bins = 40, geom = "polygon") +
theme_classic() +
xlim(c(-8, 10)) +
ylim(c(-7, 8)) +
scale_fill_viridis_c(option = "plasma") +
facet_grid(. ~ timepoint + genotype) +
theme(
strip.text.x = element_text(
size = 20, color = "black",
face = "bold"
),
strip.text.y = element_text(
size = 20, color = "black",
face = "bold"
)
)
# Test densities with unique scale for each
df.a <- df %>% filter(timepoint == "D8" & genotype == "Mock")
a <- plot(ggplot(data = df.a) +
geom_point(data = df.a, aes(x = UMAP_1, y = UMAP_2), color = "lightgrey", size = 0.5) +
stat_density_2d(aes(x = UMAP_1, y = UMAP_2, fill = after_stat(level), alpha = (..level..)), bins = 30, geom = "polygon") +
theme_classic() +
xlim(c(-8, 10)) +
ylim(c(-7, 8)) +
scale_fill_viridis_c(option = "plasma") +
theme(
strip.text.x = element_text(
size = 20, color = "black",
face = "bold"
),
strip.text.y = element_text(
size = 20, color = "black",
face = "bold"
)
))
df.b <- df %>% filter(timepoint == "D8" & genotype == "gSUV")
b <- plot(ggplot(data = df.b) +
geom_point(data = df.b, aes(x = UMAP_1, y = UMAP_2), color = "lightgrey", size = 0.5) +
stat_density_2d(aes(x = UMAP_1, y = UMAP_2, fill = after_stat(level), alpha = (..level..)), bins = 30, geom = "polygon") +
theme_classic() +
xlim(c(-8, 10)) +
ylim(c(-7, 8)) +
scale_fill_viridis_c(option = "plasma") +
theme(
strip.text.x = element_text(
size = 20, color = "black",
face = "bold"
),
strip.text.y = element_text(
size = 20, color = "black",
face = "bold"
)
))
df.c <- df %>% filter(timepoint == "D28" & genotype == "Mock")
c <- plot(ggplot(data = df.c) +
geom_point(data = df.c, aes(x = UMAP_1, y = UMAP_2), color = "lightgrey", size = 0.5) +
stat_density_2d(aes(x = UMAP_1, y = UMAP_2, fill = after_stat(level), alpha = (..level..)), bins = 30, geom = "polygon") +
theme_classic() +
xlim(c(-8, 10)) +
ylim(c(-7, 8)) +
scale_fill_viridis_c(option = "plasma") +
theme(
strip.text.x = element_text(
size = 20, color = "black",
face = "bold"
),
strip.text.y = element_text(
size = 20, color = "black",
face = "bold"
)
))
df.d <- df %>% filter(timepoint == "D28" & genotype == "gSUV")
d <- plot(ggplot(data = df.d) +
geom_point(data = df.d, aes(x = UMAP_1, y = UMAP_2), color = "lightgrey", size = 0.5) +
stat_density_2d(aes(x = UMAP_1, y = UMAP_2, fill = after_stat(level), alpha = (..level..)), bins = 30, geom = "polygon") +
theme_classic() +
xlim(c(-8, 10)) +
ylim(c(-7, 8)) +
scale_fill_viridis_c(option = "plasma") +
theme(
strip.text.x = element_text(
size = 20, color = "black",
face = "bold"
),
strip.text.y = element_text(
size = 20, color = "black",
face = "bold"
)
))
a + b + c + d + plot_layout(ncol = 4) & NoLegend()
# TCR features
"%!in%" <- function(x, y) !("%in%"(x, y))
TCRA <- read.csv("E:/Work/Signatures/TCR_genes/TCRA.txt", sep = "\t")
TCRA.features <- TCRA$Approved.symbol
TCRB <- read.csv("E:/Work/Signatures/TCR_genes/TCRB.txt", sep = "\t")
TCRB.features <- TCRB$Approved.symbol
TCRD <- read.csv("E:/Work/Signatures/TCR_genes/TCRD.txt", sep = "\t")
TCRD.features <- TCRD$Approved.symbol
TCRG <- read.csv("E:/Work/Signatures/TCR_genes/TCRG.txt", sep = "\t")
TCRG.features <- TCRG$Approved.symbol
TCR.features <- c(TCRA.features, TCRB.features, TCRG.features, TCRD.features)
TCR.removed <- rownames(suv.car.t.integrated.merged.timepoints.premerged.cd4s)[which(rownames(suv.car.t.integrated.merged.timepoints.premerged.cd4s) %!in% TCR.features)]
# Heatmap
Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s) <- factor(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s), levels = order)
Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s) <- factor(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s), levels = c(0:11))
markers <- FindAllMarkers(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, assay = "RNA", features = TCR.removed, only.pos = T)
markers %>%
group_by(cluster) %>%
filter(avg_log2FC > 0) %>%
top_n(n = 8, wt = avg_log2FC) -> top8
Average <- AverageExpression(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = unique(top8$gene), return.seurat = T, assays = "RNA")
Average <- ScaleData(object = Average, features = unique(top8$gene))
pdf("Heatmap_average_d8_d28_harmony_sample_CD4s_subcluster.pdf", width = 9, height = 20)
plot(DoHeatmap(Average, features = unique(top8$gene), draw.lines = F, size = 3, group.colors = mycolors.CD4) + scale_fill_viridis_c(option = "inferno"))
dev.off()
Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s) <- paste(suv.car.t.integrated.merged.timepoints.premerged.cd4s$RNA_snn_res.0.2, suv.car.t.integrated.merged.timepoints.premerged.cd4s$timepoint)
Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s) <- paste(suv.car.t.integrated.merged.timepoints.premerged.cd4s$RNA_snn_res.0.2, suv.car.t.integrated.merged.timepoints.premerged.cd4s$genotype)
Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s) <- suv.car.t.integrated.merged.timepoints.premerged.cd4s$RNA_snn_res.0.2
# Saving DEGs
write.csv(markers, "markers_CD4s.txt", quote = F, row.names = F)
# Transfer with NSCLC data
features <- SelectIntegrationFeatures(object.list = c(eleven.tils.cd3.integrated, suv.car.t.integrated.merged.timepoints.premerged.cd4s), nfeatures = 6000)
tumor.anchors <- FindTransferAnchors(
reference = eleven.tils.cd3.integrated, query = suv.car.t.integrated.merged.timepoints.premerged.cd4s,
dims = 1:30, query.assay = "RNA", reference.assay = "integrated"
)
predictions <- TransferData(
anchorset = tumor.anchors, refdata = Idents(eleven.tils.cd3.integrated),
dims = 1:30
)
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- AddMetaData(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, metadata = predictions)
DimPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, group.by = "predicted.id", label = T, repel = T, reduction = "umap", cols = paletteer_d("pals::polychrome", n = 21))
FeaturePlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = "prediction.score.max")
# Cell cycle scoring
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- CellCycleScoring(suv.car.t.integrated.merged.timepoints.premerged.cd4s,
g2m.features = cc.genes$g2m.genes,
s.features = cc.genes$s.genes, set.ident = F
)
DimPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, group.by = "Phase", label = T, label.size = 6, pt.size = 1.1)
# Nebulosa plotting
p <- plot_density(suv.car.t.integrated.merged.timepoints.premerged.cd4s, c("TCF7", "IL7R", "SELL", "KLF2"), joint = T, )
p + plot_layout(ncol = 3)
# CytoTRACE
library(CytoTRACE)
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s, cells = sample(Cells(suv.car.t.integrated.merged.timepoints.premerged.cd4s), 20000))
suv.car.t.integrated.merged.timepoints.premerged.cd4s.sub <- subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s, subset = Phase == "G1")
suv.car.t.integrated.merged.timepoints.premerged.cd4s[["ident_genotype_timepoint"]] <- paste(
suv.car.t.integrated.merged.timepoints.premerged.cd4s$RNA_snn_res.0.2, suv.car.t.integrated.merged.timepoints.premerged.cd4s$genotype,
suv.car.t.integrated.merged.timepoints.premerged.cd4s$timepoint
)
mat <- GetAssayData(suv.car.t.integrated.merged.timepoints.premerged.cd4s.sub, assay = "RNA")
mat <- as.data.frame(mat)
results.cyto <- CytoTRACE(mat, ncores = 1)
idents <- as.character(suv.car.t.integrated.merged.timepoints.premerged.cd4s.sub$ident_genotype_timepoint)
names(idents) <- names(suv.car.t.integrated.merged.timepoints.premerged.cd4s.sub$ident_genotype_timepoint)
emb <- as.data.frame(Embeddings(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s.sub, reduction = "umap"))
plotCytoTRACE(results.cyto, emb = emb, phenotype = idents)
# Monocle3
library(monocle3)
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s, idents = c("Mono-S100A8", "Mono-FCGR3A", "Mono-TREM1", "Early MAC-ISG15", "Early MAC-CXCR4", "MAC-CXCL2", "MAC-FBP1", "Mo-derived LAM-STAB1", "LAM-APOC1"))
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- RunUMAP(suv.car.t.integrated.merged.timepoints.premerged.cd4s,
reduction = "harmony",
dims = 1:50, reduction.name = "UMAP"
)
suv.car.t.integrated.merged.timepoints.premerged.cd4s.WT.D8 <- subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s, subset = timepoint == "D8" & genotype == "WT")
suv.car.t.integrated.merged.timepoints.premerged.cd4s.WT.D28 <- subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s, subset = timepoint == "D28" & genotype == "WT")
suv.car.t.integrated.merged.timepoints.premerged.cd4s.KO.D8 <- subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s, subset = timepoint == "D8" & genotype == "KO")
suv.car.t.integrated.merged.timepoints.premerged.cd4s.KO.D28 <- subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s, subset = timepoint == "D28" & genotype == "KO")
cds <- as.cell_data_set(suv.car.t.integrated.merged.timepoints.premerged.cd4s.KO.D8)
cds <- cluster_cells(cds)
p1 <- plot_cells(cds, show_trajectory_graph = FALSE)
p2 <- plot_cells(cds, color_cells_by = "partition", show_trajectory_graph = FALSE)
wrap_plots(p1, p2)
cds@clusters$UMAP$clusters <- Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.KO.D8)
cds <- learn_graph(cds, close_loop = F, use_partition = T, learn_graph_control = list(minimal_branch_len = 13))
p <- plot_cells(cds, cell_size = 1, label_cell_groups = F, label_groups_by_cluster = FALSE, label_leaves = FALSE, label_branch_points = FALSE, show_trajectory_graph = T, trajectory_graph_color = "black", group_label_size = 5, trajectory_graph_segment_size = 1.5, ) + ggtitle("WT D8")
p2 <- plot_cells(cds, cell_size = 1, label_cell_groups = F, label_groups_by_cluster = FALSE, label_leaves = FALSE, label_branch_points = FALSE, show_trajectory_graph = T, trajectory_graph_color = "black", group_label_size = 5, trajectory_graph_segment_size = 1.5, ) + ggtitle("KO D8")
p3 <- plot_cells(cds, cell_size = 1, label_cell_groups = F, label_groups_by_cluster = FALSE, label_leaves = FALSE, label_branch_points = FALSE, show_trajectory_graph = T, trajectory_graph_color = "black", group_label_size = 5, trajectory_graph_segment_size = 1.5, ) + ggtitle("WT D28")
p4 <- plot_cells(cds, cell_size = 1, label_cell_groups = F, label_groups_by_cluster = FALSE, label_leaves = FALSE, label_branch_points = FALSE, show_trajectory_graph = T, trajectory_graph_color = "black", group_label_size = 5, trajectory_graph_segment_size = 1.5, ) + ggtitle("KO D28")
p + p2 + p3 + p4
# Single genes nebulosa
library("Nebulosa")
plot_density(suv.car.t.integrated.merged.timepoints.premerged.cd4s, c("TOX", "PDCD1", "HAVCR2", "LAG3"))
plot_density(suv.car.t.integrated.merged.timepoints.premerged.cd4s, c("TCF7", "SELL", "IL7R", "KLF2"))
plot_density(suv.car.t.integrated.merged.timepoints.premerged.cd4s, c("NFE2L1", "NFE2L3", "JUN"))
plot_density(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, c("PRF1", "IFNG", "TNF", "ITGA1", "ITGA4", "ITGB1", "S1PR1", "CD69", "CTLA4", "NRP1"))
FeaturePlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, c("ITGA5", "ITGA6", "ITGAX", "ITGB1", "ITGB2", "ITGAM"), pt.size = 0.8, min.cutoff = "q3", max.cutoff = "q97", order = T) & scale_color_paletteer_c("pals::coolwarm")
plot_density(suv.car.t.integrated.merged.timepoints.premerged.cd4s, c("S1PR1", "KLF2", "CCR4", "CCR7", "ITGAE", "ITGA4"))
# ReactomeGSA
library(ReactomeGSA)
gsva_result <- analyse_sc_clusters(suv.car.t.integrated.merged.timepoints.premerged.cd4s, verbose = TRUE)
plot_gsva_heatmap(gsva_result, max_pathways = 20, margins = c(6, 25))
# Plot volcano KO vs WT in batch for each cluster
library(EnhancedVolcano)
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- NormalizeData(suv.car.t.integrated.merged.timepoints.premerged.cd4s) %>%
FindVariableFeatures() %>%
ScaleData()
Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s) <- paste(suv.car.t.integrated.merged.timepoints.premerged.cd4s$RNA_snn_res.0.2, suv.car.t.integrated.merged.timepoints.premerged.cd4s$timepoint)
pdf("volcano_CD4_KO_vs_WT_bycluster_timepoint.pdf", width = 9, height = 9)
for (i in levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s))) {
degs <- try(FindMarkers(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, ident.1 = "KO", ident.2 = "WT", subset.ident = i, group.by = "genotype", assay = "RNA"))
try(plot(EnhancedVolcano(degs,
lab = rownames(degs),
x = "avg_log2FC",
y = "p_val_adj",
drawConnectors = T,
pCutoff = 0.001,
FCcutoff = 0.3
) + ggtitle(paste("cluster", i, sep = " "))))
}
dev.off()
# Number of DEGs
suv.car.t.integrated.merged.timepoints.premerged.cd4s[["label"]] <- suv.car.t.integrated.merged.timepoints.premerged.cd4s$genotype
suv.car.t.integrated.merged.timepoints.premerged.cd4s[["cell_type"]] <- Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s)
results <- c()
for (i in levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s))) {
markers <- try(FindMarkers(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, ident.1 = "WT", ident.2 = "KO", subset.ident = i, group.by = "genotype", assay = "RNA"))
results[[i]] <- markers
}
lengths.deg <- c()
for (i in 1:length(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s))) {
lengths.deg[i] <- print(length(rownames(results[[i]])))
}
names(lengths.deg) <- levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s))
which(names(lengths.deg) %in% augur$AUC$cell_type)
lengths.deg <- as.data.frame(lengths.deg)
lengths.deg[["cell_type"]] <- rownames(lengths.deg)
ggplot(data = lengths.deg, aes(x = cell_type, y = lengths.deg)) +
geom_bar(stat = "identity", fill = "coral") +
theme_light() +
RotatedAxis()
# Number of DEGs split for D8
suv.car.t.integrated.merged.timepoints.premerged.cd4s.D8 <- subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s, subset = timepoint == "D8")
suv.car.t.integrated.merged.timepoints.premerged.cd4s.D8[["label"]] <- suv.car.t.integrated.merged.timepoints.premerged.cd4s.D8$genotype
suv.car.t.integrated.merged.timepoints.premerged.cd4s.D8[["cell_type"]] <- Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.D8)
results <- c()
for (i in levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.D8))) {
markers <- try(FindMarkers(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s.D8, ident.1 = "WT", ident.2 = "KO", subset.ident = i, group.by = "genotype", assay = "RNA"))
results[[i]] <- markers
}
lengths.deg <- c()
for (i in 1:length(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.D8))) {
lengths.deg[i] <- print(length(rownames(results[[i]])))
}
names(lengths.deg) <- levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.D8))
which(names(lengths.deg) %in% augur$AUC$cell_type)
lengths.deg <- as.data.frame(lengths.deg)
lengths.deg[["cell_type"]] <- rownames(lengths.deg)
ggplot(data = lengths.deg, aes(x = cell_type, y = lengths.deg)) +
geom_bar(stat = "identity", fill = "coral") +
theme_light() +
RotatedAxis()
# Proportions
gg_color_hue <- function(n) {
hues <- seq(15, 375, length = n + 1)
hcl(h = hues, l = 65, c = 100)[1:n]
}
suv.car.t.integrated.merged.timepoints.premerged.cd4s$genotype_timepoint <- paste(suv.car.t.integrated.merged.timepoints.premerged.cd4s$genotype, suv.car.t.integrated.merged.timepoints.premerged.cd4s$timepoint)
suv.car.t.integrated.merged.timepoints.premerged.cd4s$genotype_timepoint <- factor(x = suv.car.t.integrated.merged.timepoints.premerged.cd4s$genotype_timepoint, levels = c("Mock D8", "KO D8", "Mock D28", "KO D28"))
suv.car.t.integrated.merged.timepoints.premerged.cd4s$genotype <- factor(x = suv.car.t.integrated.merged.timepoints.premerged.cd4s$genotype, levels = c("WT", "KO"))
suv.car.t.integrated.merged.timepoints.premerged.cd4s$timepoint <- factor(x = suv.car.t.integrated.merged.timepoints.premerged.cd4s$timepoint, levels = c("D8", "D28"))
plot_stat(suv.car.t.integrated.merged.timepoints.premerged.cd4s, plot_type = "prop_fill", group_by = "Phase") + scale_fill_manual(values = mycolors.CD4)
# VP by timepoint / genotype
VlnPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, c("TOX", "PDCD1", "HAVCR2", "LAG3", "TIGIT", "CTLA4"), pt.size = 0.1, group.by = "genotype_timepoint") + NoLegend() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) &
geom_signif(
comparisons = list(c("WT D8", "KO D8"), c("WT D28", "KO D28")),
map_signif_level = TRUE, textsize = 6, test = "wilcox.test", y_position = 2
) & stat_summary(fun = median, geom = "point", size = 12, colour = "white", shape = 95)
VlnPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, c("TCF7", "SELL", "LEF1", "KLF2", "CCR7", "IL7R"), pt.size = 0.1, group.by = "genotype_timepoint") + NoLegend() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) &
geom_signif(
comparisons = list(c("WT D8", "KO D8"), c("WT D28", "KO D28")),
map_signif_level = TRUE, textsize = 6, test = "wilcox.test", y_position = 2.5
) & stat_summary(fun = median, geom = "point", size = 12, colour = "white", shape = 95)
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- DietSeurat(suv.car.t.integrated.merged.timepoints.premerged.cd4s)
# Label transfer of cycling
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- FindVariableFeatures(suv.car.t.integrated.merged.timepoints.premerged.cd4s, nfeatures = 4000)
suv.car.t.integrated.merged.timepoints.premerged.cd4s.ref <- subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s, subset = Phase == "G1")
suv.car.t.integrated.merged.timepoints.premerged.cd4s.ref <- subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s.ref, ident = c("C11-Cycling-CENPF", "C12-Cycling-TOP2A", "C13-Cycling-ATP5F1E"), invert = T)
suv.car.t.integrated.merged.timepoints.premerged.cd4s.query <- subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s, subset = Phase == "G1", invert = T)
tumor.anchors <- FindTransferAnchors(
reference = suv.car.t.integrated.merged.timepoints.premerged.cd4s.ref, query = suv.car.t.integrated.merged.timepoints.premerged.cd4s.query,
dims = 1:30, reference.assay = "RNA", query.assay = "RNA"
)
predictions <- TransferData(
anchorset = tumor.anchors, refdata = Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.ref),
dims = 1:30
)
suv.car.t.integrated.merged.timepoints.premerged.cd4s.query <- AddMetaData(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s.query, metadata = predictions)
DimPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s.query, group.by = "predicted.id", label = T, repel = T)
VlnPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s.query, "prediction.score.max", pt.size = 0)
suv.car.t.integrated.merged.timepoints.premerged.cd4s.query$seurat_clusters <- suv.car.t.integrated.merged.timepoints.premerged.cd4s.query$predicted.id
suv.car.t.integrated.merged.timepoints.premerged.cd4s.query$seurat_clusters <- factor(x = suv.car.t.integrated.merged.timepoints.premerged.cd4s.query$seurat_clusters, levels = c(0:10))
suv.car.t.integrated.merged.timepoints.premerged.cd4s.query$genotype <- factor(x = suv.car.t.integrated.merged.timepoints.premerged.cd4s.query$genotype, levels = c("WT", "KO"))
suv.car.t.integrated.merged.timepoints.premerged.cd4s.query$timepoint <- factor(x = suv.car.t.integrated.merged.timepoints.premerged.cd4s.query$timepoint, levels = c("D8", "D28"))
suv.car.t.integrated.merged.timepoints.premerged.cd4s.query$genotype_timepoint <- paste(suv.car.t.integrated.merged.timepoints.premerged.cd4s.query$genotype, suv.car.t.integrated.merged.timepoints.premerged.cd4s.query$timepoint)
suv.car.t.integrated.merged.timepoints.premerged.cd4s.query$genotype_timepoint <- factor(x = suv.car.t.integrated.merged.timepoints.premerged.cd4s.query$genotype_timepoint, levels = c("WT D8", "KO D8", "WT D28", "KO D28"))
plot_stat(suv.car.t.integrated.merged.timepoints.premerged.cd4s.query, plot_type = "prop_fill", group_by = "genotype_timepoint") + RotatedAxis() + scale_fill_manual(values = gg_color_hue(7))
plot_stat(suv.car.t.integrated.merged.timepoints.premerged.cd4s.query, plot_type = "prop_fill", group_by = "Phase") + RotatedAxis() + scale_fill_manual(values = gg_color_hue(7))
# Confirming that label transfer did work well
VlnPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s.query, c("KLF2", "TCF7", "NEAT1", "TMSB10", "IL7R", "SELL"), pt.size = 0, group.by = "predicted.id")
# Boxplots
freq_table <- as.data.frame(prop.table(
x = table(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s), suv.car.t.integrated.merged.timepoints.premerged.cd4s$orig.ident),
margin = 2
))
colnames(freq_table) <- c("cluster", "replicate", "Freq")
freq_table$condition[which(freq_table$replicate == "Mock.1")] <- "WT"
freq_table$condition[which(freq_table$replicate == "Mock.2")] <- "WT"
freq_table$condition[which(freq_table$replicate == "Mock.3")] <- "WT"
freq_table$condition[which(freq_table$replicate == "Suv.ko.1")] <- "KO"
freq_table$condition[which(freq_table$replicate == "Suv.ko.2")] <- "KO"
freq_table$condition[which(freq_table$replicate == "Suv.ko.3")] <- "KO"
freq_table$condition <- factor(freq_table$condition, levels = c("WT", "KO"))
freq_table$cluster <- factor(freq_table$cluster, levels = order)
pdf("boxplots_suvcart_CD4_WT_vs_KO.pdf", width = 10, height = 5)
ggplot(data = freq_table, aes(x = condition, y = Freq, fill = condition)) +
geom_boxplot() +
facet_grid(~cluster) +
theme_bw() +
geom_point(pch = 21, position = position_jitterdodge()) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 15)) +
theme(
axis.text.x = element_blank(),
axis.ticks.x = element_blank()
) +
scale_fill_manual(values = c("deepskyblue", "orange1")) +
theme(strip.text.x = element_text(size = 8, face = "bold", angle = 45)) +
geom_signif(
comparisons = list(c("WT", "KO")),
map_signif_level = F, textsize = 4, test = "wilcox.test", y_position = 0.23
)
dev.off()
order <- paste(rep(0:11, each = 2), rep(c("WT", "KO"), each = 1), rep(c("D8", "D28"), each = 12))
order <- paste(rep(0:11, each = 2), rep(c("D8", "D28"), each = 12))
order <- paste(rep(0:11, each = 2), rep(c("D8", "D28"), each = 1))
order <- paste(rep(0:11, each = 2), rep(c("WT", "KO"), each = 1))
# Boxplots split by timepoint
suv.car.t.integrated.merged.timepoints.premerged.cd4s.D8 <- subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s, subset = timepoint == "D8")
freq_table <- as.data.frame(prop.table(
x = table(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.D8), suv.car.t.integrated.merged.timepoints.premerged.cd4s.D8$orig.ident),
margin = 2
))
colnames(freq_table) <- c("cluster", "replicate", "Freq")
freq_table$condition[which(freq_table$replicate == "Mock.1")] <- "WT"
freq_table$condition[which(freq_table$replicate == "Mock.2")] <- "WT"
freq_table$condition[which(freq_table$replicate == "Mock.3")] <- "WT"
freq_table$condition[which(freq_table$replicate == "Suv.ko.1")] <- "KO"
freq_table$condition[which(freq_table$replicate == "Suv.ko.2")] <- "KO"
freq_table$condition[which(freq_table$replicate == "Suv.ko.3")] <- "KO"
freq_table$condition <- factor(freq_table$condition, levels = c("WT", "KO"))
order <- paste(c(0:11), "D8")
freq_table$cluster <- factor(freq_table$cluster, levels = order)
pdf("boxplots_suvcart_cd4_WT_vs_KO.pdf", width = 10, height = 5)
ggplot(data = freq_table, aes(x = condition, y = Freq, fill = condition)) +
geom_boxplot() +
facet_grid(~cluster) +
theme_bw() +
geom_point(pch = 21, position = position_jitterdodge()) +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 15)) +
theme(
axis.text.x = element_blank(),
axis.ticks.x = element_blank()
) +
scale_fill_manual(values = c("deepskyblue", "orange1")) +
theme(strip.text.x = element_text(size = 8, face = "bold", angle = 45)) +
geom_signif(
comparisons = list(c("WT", "KO")),
map_signif_level = F, textsize = 4, test = "wilcox.test", test.args = c(alternative = "two.sided"), y_position = 0.26
)
dev.off()
# FP for QC
FeaturePlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, c("CD3D", "cd4A", "cd4B", "CD4"), pt.size = 0.8, order = F, max.cutoff = "q97") & scale_color_paletteer_c("pals::coolwarm") & NoAxes()
FeaturePlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, c("nFeature_RNA", "nCount_RNA", "percent.mt", "percent.rb"), pt.size = 0.8, order = F, max.cutoff = "q97") & scale_color_paletteer_c("pals::coolwarm") & NoAxes()
# Splitting cluster 0
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- FindSubCluster(suv.car.t.integrated.merged.timepoints.premerged.cd4s, cluster = "0", resolution = 0.1, graph.name = "RNA_snn")
DimPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, group.by = "sub.cluster", reduction = "umap", label = TRUE, pt.size = 1.3, label.size = 6) + ggtitle("CD4s / Day 8 + 28 / n = 49728")
# Plotting signature scores
all_sign2 <- all_sign
all_sign2 <- paste0(all_sign, "1")
for (i in 1:length(all_sign)) {
try(suv.car.t.integrated.merged.timepoints.premerged.cd4s <- AddModuleScore(suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = list(get(all_sign[i])), name = all_sign[i], assay = "RNA"))
}
# As Featureplots
pdf("Signatures_CD4s_Fraietta.pdf", width = 12, height = 12)
for (i in 1:length(all_sign2)) {
try(plot(FeaturePlot(
object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = all_sign2[i], order = T, reduction = "umap",
pt.size = 1.5, min.cutoff = "q1", max.cutoff = "q99"
) + labs(title = all_sign[i], size = 13) + scale_color_paletteer_c("pals::coolwarm")))
}
dev.off()
# Comparison with Chen et al. data
Chen2021_long_response <- readLines("E:/Work/Signatures/Chen_et_al_2021/Long_response_all.txt")
Chen2021_long_response_top50_FC <- readLines("E:/Work/Signatures/Chen_et_al_2021/Long_response_all_top50.txt")
Chen2021_short_response <- readLines("E:/Work/Signatures/Chen_et_al_2021/short_response_all.txt")
fraietta_cr <- readLines("E:/Work/Signatures/Fraietta_et_al_2018/Fraietta_CR.txt")
fraietta_nr <- readLines("E:/Work/Signatures/Fraietta_et_al_2018/Fraietta_NR.txt")
# As violinplots between genotype and time point
all_sign <- c("Chen2021_long_response", "Chen2021_long_response_top50_FC", "Chen2021_short_response")
all_sign <- c("fraietta_cr", "fraietta_nr")
all_sign <- c("sign_ox_ph", "sign_glycolysis", "sign_FA_metabolism_GO", "sign_fatty_acid_metab_KEGG")
pdf("VP_car_t_CD4_Fra_full.pdf", width = 5, height = 6)
for (i in 1:length(all_sign2)) {
try(plot(VlnPlot(subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s, subset = timepoint == "D8"), "fraietta_nr1", group.by = "genotype", pt.size = 0.1) + labs(title = "Fraietta NR", size = 13) + ylim(-0.5, 0.2) + NoLegend() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) &
geom_signif(
comparisons = list(c("WT", "KO")),
map_signif_level = TRUE, textsize = 6, test = "wilcox.test", y_position = 0.1
) & stat_summary(fun = median, geom = "point", size = 12, colour = "white", shape = 95)))
}
dev.off()
# FP
FeaturePlot(
object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = "fraietta_nr1", split.by = "genotype_timepoint", order = T, reduction = "umap",
pt.size = 1.5, min.cutoff = "q1", max.cutoff = "q99"
) & scale_color_paletteer_c("pals::coolwarm")
Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s) <- paste(suv.car.t.integrated.merged.timepoints.premerged.cd4s$genotype, suv.car.t.integrated.merged.timepoints.premerged.cd4s$timepoint)
Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s) <- factor(x = Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s), levels = c("WT D8", "KO D8", "WT D28", "KO D28"))
Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s) <- paste(
suv.car.t.integrated.merged.timepoints.premerged.cd4s$RNA_snn_res.0.2,
suv.car.t.integrated.merged.timepoints.premerged.cd4s$genotype, suv.car.t.integrated.merged.timepoints.premerged.cd4s$timepoint
)
Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s) <- factor(x = Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s), levels = order)
order <- paste(rep(0:11, each = 4), rep(c("WT", "KO"), each = 1), rep(c("D8", "D28"), each = 2))
# muscat test
library(muscat)
suv.car.t.integrated.merged.timepoints.premerged.cd4s.ref <- subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s, ident = c("C11-Cycling-CENPF", "C12-Cycling-TOP2A", "C13-Cycling-ATP5F1E"), invert = T)
suv.car.t.integrated.merged.timepoints.premerged.cd4s.ref$seurat_clusters <- Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.ref)
sce <- as.SingleCellExperiment(suv.car.t.integrated.merged.timepoints.premerged.cd4s.ref, assay = "RNA")
sce$id <- paste(sce$timepoint, sce$orig.ident, sep = "_")
muscat.sce <- prepSCE(sce,
kid = "seurat_clusters", # I have no cluter id
gid = "genotype_timepoint", # group IDs (ctrl/stim)
sid = "id", drop = T
)
pb <- aggregateData(muscat.sce,
assay = "counts", fun = "sum",
by = c("cluster_id", "sample_id")
)
(pb_mds <- pbMDS(pb))
# sc T cell naive/stemness signature
cd4_stem_sc <- readLines("E:/Work/Signatures/Stemness/cd4_Stem_like_custom_single_cell.txt")
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- AddModuleScore(suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = list(cd4_stem_sc), name = "cd4_stem_sc", assay = "RNA")
# As Featureplots
FeaturePlot(
object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = "CD8_stem_sc1", split.by = "genotype", order = T, reduction = "umap",
pt.size = 1.5, min.cutoff = "q1", max.cutoff = "q99"
) & scale_color_paletteer_c("pals::coolwarm")
plot_density(suv.car.t.integrated.merged.timepoints.premerged.cd4s, "CD8_stem_sc1")
# sc T cell naive/stemness signature
cd4_naive_goldrath <- toupper(readLines("E:/Work/Signatures/Naive/GOLDRATH_NAIVE_VS_EFF_cd4_TCELL_UP.txt"))
cd4_naive_Guo <- readLines("E:/Work/Signatures/Guo_et_al_2018/cd4_LEF1.txt")
CD4_naive_Guo <- readLines("E:/Work/Signatures/Guo_et_al_2018/CD4_CCR7.txt")
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- AddModuleScore(suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = list(cd4_naive_goldrath), name = "cd4_naive_goldrath", assay = "RNA")
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- AddModuleScore(suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = list(cd4_naive_Guo), name = "cd4_naive_Guo", assay = "RNA")
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- AddModuleScore(suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = list(CD4_naive_Guo), name = "CD4_naive_Guo", assay = "RNA")
# As Featureplots
p <- FeaturePlot(
object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = "cd4_naive_goldrath1", order = T, reduction = "umap",
pt.size = 1, min.cutoff = "q1", max.cutoff = "q99"
) + labs(title = "cd4+ naive goldrath") + scale_color_paletteer_c("pals::coolwarm")
p1 <- FeaturePlot(
object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = "cd4_naive_Guo1", order = T, reduction = "umap",
pt.size = 1, min.cutoff = "q1", max.cutoff = "q99"
) + labs(title = "cd4_naive_Guo") + scale_color_paletteer_c("pals::coolwarm")
p2 <- FeaturePlot(
object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = "CD4_naive_Guo1", order = T, reduction = "umap",
pt.size = 1, min.cutoff = "q1", max.cutoff = "q99"
) + labs(title = "CD4_naive_Guo") + scale_color_paletteer_c("pals::coolwarm")
p + p1 + p2 + plot_layout(ncol = 3)
# Fraietta et al. 2018 FGI2C signatures
Fraietta.2C <- readxl::read_xlsx(path = "E:/Work/Signatures/Fraietta_et_al_2018/Fig2C_41591_2018_10_MOESM5_ESM.xlsx")
colnames(Fraietta.2C) <- gsub("[\r\n]", "", colnames(Fraietta.2C))
colnames(Fraietta.2C) <- make.names(colnames(Fraietta.2C))
pdf("Fraietta.2C_SUV_CD4.pdf", height = 10, width = 11)
for (i in colnames(Fraietta.2C)) {
suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned <- AddModuleScore(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, features = list(Fraietta.2C[[i]]), name = i, assay = "RNA")
# plot(FeaturePlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, paste0(i, 1), pt.size = 0.8, min.cutoff = "q3", max.cutoff = "q97", order = T) + scale_color_paletteer_c("pals::coolwarm") +ggtitle(i))
}
dev.off()
# Heatmap with signatures
colnames(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned@meta.data) <- gsub(pattern = 1, replacement = "", x = colnames(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned@meta.data))
suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned[["Fraietta"]] <- CreateAssayObject(data = t(x = FetchData(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, vars = colnames(Fraietta.2C))))
Average <- AverageExpression(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, features = colnames(Fraietta.2C), assays = "Fraietta", slot = "data", return.seurat = T)
Average <- ScaleData(object = Average, features = colnames(Fraietta.2C))
pdf("Heatmap_average_Fraietta_CD4s.pdf", width = 9, height = 12)
plot(DoHeatmap(Average, features = colnames(Fraietta.2C), draw.lines = F, size = 3) + scale_fill_viridis_c(option = "inferno"))
dev.off()
# Reorder idents
suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned$reordered.idents <- paste(
Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned),
suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned$genotype, suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned$timepoint
)
order <- paste(rep(c(
"C1-Active translation-IL7R", "C2-Stem/memory-KLF2", "C3-Stem/cytotoxic-GNLY", "C4-Resident-ITGA1", "C5-Activated-CD69", "C6-Activated-JUN", "C7-Active translation-RPS27",
"C8-Effector-KLRB1", "C9-Cytosolic stress-HSPA1B", "C10-Cytosolic stress-LncRNA", "C11-Pro-apopototic-BAX", "C12-Cycling/G2M-HSPA1B", "C13-Cycling-TOP2A", "C14-Cycling/S-PCNA"
), each = 1), rep(c("Mock", "gSUV"), each = 14), rep(c("D8", "D28"), each = 28))
suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned$reordered.idents <- factor(x = suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned$reordered.idents, levels = order)
# Redo same with only cluster 2-memory-KLF2
suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned.KLF2 <- subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, idents = "C2-Stem/memory-KLF2")
Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned.KLF2) <- suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned.KLF2$reordered.idents
colnames(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned.KLF2@meta.data) <- gsub(pattern = 1, replacement = "", x = colnames(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned.KLF2@meta.data))
suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned.KLF2[["Fraietta"]] <- CreateAssayObject(data = t(x = FetchData(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned.KLF2, vars = colnames(Fraietta.2C))))
Average <- AverageExpression(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned.KLF2, features = colnames(Fraietta.2C), assays = "Fraietta", slot = "data", return.seurat = T)
Average <- ScaleData(object = Average, features = colnames(Fraietta.2C))
plot(DoHeatmap(Average, features = colnames(Fraietta.2C), draw.lines = F, size = 4) + scale_fill_viridis_c(option = "inferno") + theme(axis.text.y = element_text(size = 12)))
# SignacX classification
library(SignacX)
labels <- Signac(suv.car.t.integrated.merged.timepoints.premerged.cd4s, num.cores = 14)
celltypes <- GenerateLabels(labels, E = suv.car.t.integrated.merged.timepoints.premerged.cd4s)
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- AddMetaData(suv.car.t.integrated.merged.timepoints.premerged.cd4s, metadata = celltypes$CellStates_novel, col.name = "CellStates_novel")
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- SetIdent(suv.car.t.integrated.merged.timepoints.premerged.cd4s, value = "CellStates_novel")
DimPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, pt.size = 0.9, label = T)
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- AddMetaData(suv.car.t.integrated.merged.timepoints.premerged.cd4s, metadata = celltypes$L8, col.name = "L8")
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- SetIdent(suv.car.t.integrated.merged.timepoints.premerged.cd4s, value = "L18")
DimPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, pt.size = 0.9, group.by = "L8")
# Label transfer Guo et al.
tumor.anchors <- FindTransferAnchors(
reference = guo.seurat, query = suv.car.t.integrated.merged.timepoints.premerged.cd4s,
dims = 1:30
)
predictions <- TransferData(
anchorset = tumor.anchors, refdata = guo.seurat$majorCluster,
dims = 1:30
)
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- AddMetaData(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, metadata = predictions)
DimPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, group.by = "predicted.id", label = T)
# RoE heatmap
library("pheatmap")
suv.car.t.integrated.merged.timepoints.premerged.cd4s$genotype_timepoint <- factor(x = suv.car.t.integrated.merged.timepoints.premerged.cd4s$genotype_timepoint, levels = c("Mock D8", "KO D8", "Mock D28", "KO D28"))
clusters.table <- table(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s), suv.car.t.integrated.merged.timepoints.premerged.cd4s$genotype_timepoint)
chi.prop <- chisq.test(clusters.table)
Roe <- chi.prop$observed / chi.prop$expected
pheatmap(Roe, cluster_rows = F, cluster_cols = F, scale = "column", display_numbers = F, fontsize = 14)
# Roe heatmap split by D8 and D28
suv.car.t.integrated.merged.timepoints.premerged.cd4s$genotype <- factor(x = suv.car.t.integrated.merged.timepoints.premerged.cd4s$genotype, levels = c("Mock", "KO"))
suv.car.t.integrated.merged.timepoints.premerged.cd4s.D8 <- subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s, subset = timepoint == "D8")
clusters.table <- table(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.D8), suv.car.t.integrated.merged.timepoints.premerged.cd4s.D8$genotype)
chi.prop <- chisq.test(clusters.table)
Roe <- chi.prop$observed / chi.prop$expected
pheatmap(Roe, cluster_rows = F, cluster_cols = F, scale = "none", display_numbers = F, fontsize = 14)
# DEenrichR plot
library(enrichR)
dbs <- listEnrichrDbs()
dbs <- c("KEGG_2021_Human", "WikiPathway_2021_Human", "MSigDB_Hallmark_2020", "GO_Biological_Process_2021", "Elsevier_Pathway_Collection", "HuBMAP_ASCT_plus_B_augmented_w_RNAseq_Coexpression")
dbs <- c("MSigDB_Hallmark_2020")
pdf(file = "Pathways_CD4s_DEenrichR_KEGG_2021_Human.pdf", width = 9, height = 8)
for (i in levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s))) {
plot(DEenrichRPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, ident.1 = "C5-Memory-KLF2", logfc.threshold = 0.2, balanced = T, assay = "RNA", enrich.database = dbs, max.genes = 300))
}
dev.off()
# DEenrichR plot for WT vs KO
library(enrichR)
dbs <- listEnrichrDbs()
dbs <- c("KEGG_2021_Human", "WikiPathway_2021_Human", "MSigDB_Hallmark_2020", "GO_Biological_Process_2021", "GO_Molecular_Function_2021", "Elsevier_Pathway_Collection", "HuBMAP_ASCT_plus_B_augmented_w_RNAseq_Coexpression")
# ROGUE index by genotype and by cluster
library(ROGUE)
expr <- suv.car.t.integrated.merged.timepoints.premerged.cd4s@assays$RNA@counts
expr <- as.matrix(expr)
rogue.res <- rogue(expr,
labels = Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s),
samples = suv.car.t.integrated.merged.timepoints.premerged.cd4s$orig.ident,
platform = "UMI",
filter = T
)
rogue.plot <- rogue.boxplot(rogue.res)
rogue.plot + RotatedAxis()
# Comparison with Leticia's data
new.symbols <- toupper(leti.integrated.sub3@assays$RNA@counts@Dimnames[[1]])
leti.integrated.sub3@assays$RNA@counts@Dimnames[[1]] <- new.symbols
new.symbols <- toupper(leti.integrated.sub3@assays$RNA@data@Dimnames[[1]])
leti.integrated.sub3@assays$RNA@data@Dimnames[[1]] <- new.symbols
new.symbols <- toupper(leti.integrated.sub3@assays$integrated@data@Dimnames[[1]])
leti.integrated.sub3@assays$integrated@data@Dimnames[[1]] <- new.symbols
Hum_gene_assay <- CreateAssayObject(GetAssayData(leti.integrated.sub3, assay = "integrated"))
leti.integrated.sub3[["hum_gene_assay"]] <- Hum_gene_assay
features <- SelectIntegrationFeatures(object.list = c(leti.integrated.sub3, suv.car.t.integrated.merged.timepoints.premerged.cd4s), features = 5000)
tumor.anchors <- FindTransferAnchors(
reference = leti.integrated.sub3, query = suv.car.t.integrated.merged.timepoints.premerged.cd4s,
dims = 1:50, reference.assay = "RNA", query.assay = "RNA", reduction = "cca", features = features
)
predictions <- TransferData(
anchorset = tumor.anchors, refdata = Idents(leti.integrated.sub3),
dims = 1:50, weight.reduction = "cca"
)
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- AddMetaData(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, metadata = predictions)
DimPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, group.by = "predicted.id", label = T, pt.size = 0.8) + scale_color_brewer(palette = "Set2")
# Cytotoxic signatures
Guo_CD4_GZMA <- readLines("E:/Work/Signatures/Guo_et_al_2018/CD4_GZMA.txt")
Guo_cd4_CX3CR1 <- readLines("E:/Work/Signatures/Guo_et_al_2018/cd4_CX3CR1.txt")
Gueguen_CD4_GZMA <- readLines("E:/Work/Signatures/Gueguen_et_al_2021/CD4-GZMA.txt")
Gueguen_CD8_FCGR3A <- readLines("E:/Work/Signatures/Gueguen_et_al_2021/cd4-FCGR3A.txt")
Biocarta_cytotoxic <- readLines("E:/Work/Signatures/Cytotoxic/BIOCARTA_TCYTOTOXIC_PATHWAY.txt")
GO_BP_cytotoxic <- readLines("E:/Work/Signatures/Cytotoxic/GOBP_T_CELL_MEDIATED_CYTOTOXICITY.txt")
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- AddModuleScore(suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = list(Guo_CD4_GZMA), name = "Guo_CD4_GZMA", assay = "RNA")
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- AddModuleScore(suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = list(Guo_cd4_CX3CR1), name = "Guo_cd4_CX3CR1", assay = "RNA")
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- AddModuleScore(suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = list(Gueguen_CD4_GZMA), name = "Gueguen_CD4_GZMA", assay = "RNA")
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- AddModuleScore(suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = list(Gueguen_CD8_FCGR3A), name = "Gueguen_CD8_FCGR3A", assay = "RNA")
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- AddModuleScore(suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = list(Biocarta_cytotoxic), name = "Biocarta_cytotoxic", assay = "RNA")
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- AddModuleScore(suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = list(GO_BP_cytotoxic), name = "GO_BP_cytotoxic", assay = "RNA")
# As Featureplots
p <- FeaturePlot(
object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = "Guo_CD4_GZMA1", order = T, reduction = "umap",
pt.size = 0.7, min.cutoff = "q1", max.cutoff = "q99"
) + labs(title = "Guo_CD4_GZMA") + scale_color_paletteer_c("pals::coolwarm")
p1 <- FeaturePlot(
object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = "Guo_cd4_CX3CR11", order = T, reduction = "umap",
pt.size = 0.7, min.cutoff = "q1", max.cutoff = "q99"
) + labs(title = "Guo_cd4_CX3CR1") + scale_color_paletteer_c("pals::coolwarm")
p2 <- FeaturePlot(
object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = "Biocarta_cytotoxic1", order = T, reduction = "umap",
pt.size = 0.7, min.cutoff = "q1", max.cutoff = "q99"
) + labs(title = "Biocarta_cytotoxic") + scale_color_paletteer_c("pals::coolwarm")
p3 <- FeaturePlot(
object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = "GO_BP_cytotoxic1", order = T, reduction = "umap",
pt.size = 0.7, min.cutoff = "q1", max.cutoff = "q99"
) + labs(title = "GO_BP_cytotoxic") + scale_color_paletteer_c("pals::coolwarm")
p4 <- FeaturePlot(
object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = "Gueguen_CD4_GZMA1", order = T, reduction = "umap",
pt.size = 0.7, min.cutoff = "q1", max.cutoff = "q99"
) + labs(title = "Gueguen_CD4_GZMA") + scale_color_paletteer_c("pals::coolwarm")
p5 <- FeaturePlot(
object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = "Gueguen_CD8_FCGR3A1", order = T, reduction = "umap",
pt.size = 0.7, min.cutoff = "q1", max.cutoff = "q99"
) + labs(title = "Gueguen_CD8_FCGR3A") + scale_color_paletteer_c("pals::coolwarm")
p + p1 + p2 + p3 + p4 + p5 + plot_layout(ncol = 3)
# As Featureplots split by genotype
p <- FeaturePlot(
object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = "Guo_CD4_GZMA1", split.by = "genotype", order = T, reduction = "umap",
pt.size = 0.7, min.cutoff = "q1", max.cutoff = "q99"
) & scale_color_paletteer_c("pals::coolwarm")
p1 <- FeaturePlot(
object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = "Guo_cd4_CX3CR11", split.by = "genotype", order = T, reduction = "umap",
pt.size = 0.7, min.cutoff = "q1", max.cutoff = "q99"
) & scale_color_paletteer_c("pals::coolwarm")
p2 <- FeaturePlot(
object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = "Biocarta_cytotoxic1", split.by = "genotype", order = T, reduction = "umap",
pt.size = 0.7, min.cutoff = "q1", max.cutoff = "q99"
) & scale_color_paletteer_c("pals::coolwarm")
p3 <- FeaturePlot(
object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = "GO_BP_cytotoxic1", split.by = "genotype", order = T, reduction = "umap",
pt.size = 0.7, min.cutoff = "q1", max.cutoff = "q99"
) & scale_color_paletteer_c("pals::coolwarm")
p4 <- FeaturePlot(
object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = "Gueguen_CD4_GZMA1", split.by = "genotype", order = T, reduction = "umap",
pt.size = 0.7, min.cutoff = "q1", max.cutoff = "q99"
) & scale_color_paletteer_c("pals::coolwarm")
p5 <- FeaturePlot(
object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, features = "Gueguen_CD8_FCGR3A1", split.by = "genotype", order = T, reduction = "umap",
pt.size = 0.7, min.cutoff = "q1", max.cutoff = "q99"
) & scale_color_paletteer_c("pals::coolwarm")
p + p1 + p2 + p3 + p4 + p5 + plot_layout(ncol = 1, nrow = 7)
VlnPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, c("Guo_CD4_GZMA1", "Guo_cd4_CX3CR11", "Biocarta_cytotoxic1", "GO_BP_cytotoxic1", "Gueguen_CD4_GZMA1", "Gueguen_CD8_FCGR3A1"), pt.size = 0.1) + NoLegend() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) &
geom_signif(
comparisons = list(c("WT D8", "KO D8"), c("WT D28", "KO D28")),
map_signif_level = TRUE, textsize = 6, test = "wilcox.test", y_position = 5
) & stat_summary(fun = median, geom = "point", size = 12, colour = "white", shape = 95)
# Comparison with Luigia's data
new.symbols <- toupper(suv.nWTp@assays$RNA@counts@Dimnames[[1]])
suv.nWTp@assays$RNA@counts@Dimnames[[1]] <- new.symbols
new.symbols <- toupper(suv.nWTp@assays$RNA@data@Dimnames[[1]])
suv.nWTp@assays$RNA@data@Dimnames[[1]] <- new.symbols
features <- SelectIntegrationFeatures(object.list = c(suv.nWTp, suv.car.t.integrated.merged.timepoints.premerged.cd4s), features = 5000)
tumor.anchors <- FindTransferAnchors(
reference = suv.nWTp, query = suv.car.t.integrated.merged.timepoints.premerged.cd4s,
dims = 1:50, reference.assay = "RNA", query.assay = "RNA", reduction = "cca", features = features
)
predictions <- TransferData(
anchorset = tumor.anchors, refdata = Idents(suv.nWTp),
dims = 1:50, weight.reduction = "cca"
)
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- AddMetaData(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s, metadata = predictions)
DimPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, group.by = "predicted.id", label = T, pt.size = 0.8) + scale_color_brewer(palette = "Set2")
# Heatmap of number of DEGs (split by timepoint)
suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned[["label"]] <- suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned$genotype
suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned[["cell_type"]] <- Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned)
suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned.D8 <- subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, subset = timepoint == "D8")
results.D8 <- c()
for (i in levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned.D8))[1:11]) {
markers <- try(FindMarkers(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned.D8, logfc.threshold = 0.2, ident.1 = "Mock", ident.2 = "gSUV", subset.ident = i, group.by = "genotype", assay = "RNA"))
results.D8[[i]] <- markers
results.D8[[i]] <- results.D8[[i]] %>% filter(p_val < 0.01)
}
lengths.deg.D8 <- c()
for (i in levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned.D8))[1:11]) {
lengths.deg.D8[i] <- print(length(rownames(results.D8[[i]])))
}
names(lengths.deg.D8) <- levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned.D8))[1:11]
lengths.deg.D8 <- as.data.frame(lengths.deg.D8)
lengths.deg.D8[["cell_type"]] <- rownames(lengths.deg.D8)
# Same for D28
suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned.D28 <- subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, subset = timepoint == "D28")
results.D28 <- c()
for (i in levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned.D8))[1:11]) {
markers <- try(FindMarkers(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned.D28, logfc.threshold = 0.2, ident.1 = "Mock", ident.2 = "gSUV", subset.ident = i, group.by = "genotype", assay = "RNA"))
try(results.D28[[i]] <- markers)
try(results.D28[[i]] <- results.D28[[i]] %>% filter(p_val < 0.01))
}
lengths.deg.D28 <- c()
for (i in levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned.D8))[1:11]) {
lengths.deg.D28[i] <- print(length(rownames(results.D28[[i]])))
}
names(lengths.deg.D28) <- levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned.D28))[1:11]
lengths.deg.D28 <- as.data.frame(lengths.deg.D28)
lengths.deg.D28[["cell_type"]] <- rownames(lengths.deg.D28)
lengths.deg.D8 <- cbind(lengths.deg.D8, lengths.deg.D28)
lengths.deg.D8[, c(2, 4)] <- NULL
df <- as.data.frame(lengths.deg.D8)
df <- data.matrix(df)
pheatmap(df, cluster_rows = T, cluster_cols = T, scale = "none", fontsize = 13)
# Total
for (i in levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned))[1:11]) {
markers <- try(FindMarkers(object = suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, logfc.threshold = 0.2, ident.1 = "Mock", ident.2 = "gSUV", subset.ident = i, group.by = "genotype", assay = "RNA"))
try(results.D28[[i]] <- markers)
try(results.D28[[i]] <- results.D28[[i]] %>% filter(p_val < 0.01))
}
lengths.deg.D28 <- c()
for (i in levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned))[1:11]) {
lengths.deg.D28[i] <- print(length(rownames(results.D28[[i]])))
}
names(lengths.deg.D28) <- levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned))[1:11]
lengths.deg.D28 <- as.data.frame(lengths.deg.D28)
lengths.deg.D28[["cell_type"]] <- rownames(lengths.deg.D28)
df <- as.data.frame(lengths.deg.D28)
df <- data.matrix(df)
df <- df[-1,-1]
pheatmap(df, cluster_rows = T, cluster_cols = T, scale = "none", fontsize = 13)
# Dimplot with cluster names shown
Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s) <- suv.car.t.integrated.merged.timepoints.premerged.cd4s$sub.cluster
suv.car.t.integrated.merged.timepoints.premerged.cd4s <- RenameIdents(suv.car.t.integrated.merged.timepoints.premerged.cd4s,
"0_0" = "C1-Naivelike-IL7R", "0_1" = "C5-Resident-ITGA1", "1" = "C4-Resident-CD69",
"2" = "C12-Cycling-TOP2A", "3" = "C3-Stem/memory-KLF2", "4" = "C13-Cycling-ATP5F1E", "5" = "C7-Cytotoxic-GNLY", "6" = "C6-Effector/ Memory-SNHG25",
"7" = "C11-Cycling-CENPF", "8" = "C8-Activated-HSPA1B", "9" = "C2-Naivelike-RPS27",
"10" = "C9-Activated-PLCG2", "11" = "C10-Pro-apopototic-BAX"
)
Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s) <- factor(x = Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s), levels = c(
"C1-Naivelike-IL7R", "C2-Naivelike-RPS27", "C3-Stem/memory-KLF2", "C4-Resident-CD69", "C5-Resident-ITGA1",
"C6-Effector/ Memory-SNHG25", "C7-Cytotoxic-GNLY", "C8-Activated-HSPA1B",
"C9-Activated-PLCG2", "C10-Pro-apopototic-BAX", "C11-Cycling-CENPF", "C12-Cycling-TOP2A", "C13-Cycling-ATP5F1E"
))
# Colors setup
mycolors <- pals::kelly(16)
mycolors.CD4 <- mycolors[4:16]
DimPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, pt.size = 0.9, label = T, repel = T, label.size = 4) + NoLegend() + scale_color_manual(values = mycolors.CD4)
# Selected FP
FeaturePlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s, c("CD8_stem_sc1", "Guo_cd4_CX3CR11", "sign_effector1", "sign_cycling1", "sign_exhausted1", "sign_ICP_complete21"), pt.size = 0.8, min.cutoff = "q1", max.cutoff = "q99", order = T) & scale_color_paletteer_c("pals::coolwarm") & NoAxes()
# VP of cytotoxic and stem in cluster 3 (CD4 and/or cd4s)
suv.car.t.integrated.merged.timepoints.premerged.cd4s.sub <- subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s, idents = c("C5-Stem/memory-KLF2"))
suv.car.t.integrated.merged.timepoints.premerged.cd4s.sub$genotype_timepoint <- factor(x = suv.car.t.integrated.merged.timepoints.premerged.cd4s.sub$genotype_timepoint, levels = c("Mock D8", "KO D8", "Mock D28", "KO D28"))
VlnPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s.sub, "CD8_stem_sc1", group.by = "genotype_timepoint", pt.size = 0) + NoLegend() + stat_summary(fun = median, geom = "point", size = 12, colour = "white", shape = 95) +
geom_signif(
comparisons = list(c("Mock D8", "KO D8"), c("Mock D28", "KO D28")),
map_signif_level = TRUE, textsize = 6, test = "wilcox.test", y_position = 0.8
)
VlnPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s.sub, "Gueguen_CD8_FCGR3A1", group.by = "genotype_timepoint", pt.size = 0) + NoLegend() + stat_summary(fun = median, geom = "point", size = 12, colour = "white", shape = 95) +
geom_signif(
comparisons = list(c("Mock D8", "KO D8"), c("Mock D28", "KO D28")),
map_signif_level = TRUE, textsize = 6, test = "wilcox.test", y_position = 0.6
)
VlnPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s.sub, "sign_ox_ph1", group.by = "genotype_timepoint", pt.size = 0) + NoLegend() + stat_summary(fun = median, geom = "point", size = 12, colour = "white", shape = 95) +
geom_signif(
comparisons = list(c("Mock D8", "KO D8"), c("Mock D28", "KO D28")),
map_signif_level = TRUE, textsize = 6, test = "wilcox.test", y_position = 0.7
)
VlnPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s.sub, "sign_glycolysis1", group.by = "genotype_timepoint", pt.size = 0) + NoLegend() + stat_summary(fun = median, geom = "point", size = 12, colour = "white", shape = 95) +
geom_signif(
comparisons = list(c("Mock D8", "KO D8"), c("Mock D28", "KO D28")),
map_signif_level = TRUE, textsize = 6, test = "wilcox.test", y_position = 0.21
)
VlnPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s.sub, "sign_FA_metabolism_GO1", group.by = "genotype_timepoint", pt.size = 0) + NoLegend() + stat_summary(fun = median, geom = "point", size = 12, colour = "white", shape = 95) +
geom_signif(
comparisons = list(c("Mock D8", "KO D8"), c("Mock D28", "KO D28")),
map_signif_level = TRUE, textsize = 6, test = "wilcox.test", y_position = 0.04
)
# Metabolism single genes
VlnPlot(suv.car.t.integrated.merged.timepoints.premerged.cd4s.sub, "PPA1", group.by = "genotype_timepoint", pt.size = 0) + NoLegend() + stat_summary(fun = median, geom = "point", size = 12, colour = "white", shape = 95) +
geom_signif(
comparisons = list(c("Mock D8", "KO D8"), c("Mock D28", "KO D28")),
map_signif_level = TRUE, textsize = 6, test = "wilcox.test", y_position = 2
)
# Scatterplot cytotoxic / stem-like
c <- plot(FeatureScatter(subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, idents = levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned))[2], subset = genotype == "Mock"), feature1 = "Gueguen_CD8_FCGR3A", feature2 = "CD8_stem_sc", group.by = "genotype")) + xlim(c(-0, 1)) + ylim(c(-1, 1.5)) + geom_hline(yintercept = 0) + geom_vline(xintercept = 0)
d <- plot(FeatureScatter(subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, idents = levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned))[2], subset = genotype == "KO"), feature1 = "Gueguen_CD8_FCGR3A", feature2 = "CD8_stem_sc", group.by = "genotype")) + xlim(c(-0, 1)) + ylim(c(-1, 1.5)) + geom_hline(yintercept = 0) + geom_vline(xintercept = 0)
e <- plot(FeatureScatter(subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, idents = levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned))[3], subset = genotype == "Mock"), feature1 = "Gueguen_CD8_FCGR3A", feature2 = "CD8_stem_sc", group.by = "genotype")) + xlim(c(-0, 1)) + ylim(c(-1, 1.5)) + geom_hline(yintercept = 0) + geom_vline(xintercept = 0)
f <- plot(FeatureScatter(subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, idents = levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned))[3], subset = genotype == "KO"), feature1 = "Gueguen_CD8_FCGR3A", feature2 = "CD8_stem_sc", group.by = "genotype")) + xlim(c(-0, 1)) + ylim(c(-1, 1.5)) + geom_hline(yintercept = 0) + geom_vline(xintercept = 0)
g <- plot(FeatureScatter(subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, idents = levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned))[13], subset = genotype == "Mock"), feature1 = "Gueguen_CD8_FCGR3A", feature2 = "CD8_stem_sc", group.by = "genotype")) + xlim(c(-0, 1)) + ylim(c(-1, 1.5)) + geom_hline(yintercept = 0) + geom_vline(xintercept = 0)
h <- plot(FeatureScatter(subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, idents = levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned))[13], subset = genotype == "KO"), feature1 = "Gueguen_CD8_FCGR3A", feature2 = "CD8_stem_sc", group.by = "genotype")) + xlim(c(-0, 1)) + ylim(c(-1, 1.5)) + geom_hline(yintercept = 0) + geom_vline(xintercept = 0)
pdf("FS_CD4_sign.pdf", width = 9, height = 20)
plot((c | d) / (e | f) + (g | h))
dev.off()
# Scatterplot cytotoxic / stem-like
c <- plot(FeatureScatter(subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, idents = levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned))[2], subset = genotype == "Mock"), feature1 = "GNLY", feature2 = "KLF2", group.by = "genotype"))
d <- plot(FeatureScatter(subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, idents = levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned))[2], subset = genotype == "KO"), feature1 = "GNLY", feature2 = "KLF2", group.by = "genotype"))
e <- plot(FeatureScatter(subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, idents = levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned))[3], subset = genotype == "Mock"), feature1 = "GNLY", feature2 = "KLF2", group.by = "genotype"))
f <- plot(FeatureScatter(subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, idents = levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned))[3], subset = genotype == "KO"), feature1 = "GNLY", feature2 = "KLF2", group.by = "genotype"))
g <- plot(FeatureScatter(subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, idents = levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned))[13], subset = genotype == "Mock"), feature1 = "GNLY", feature2 = "KLF2", group.by = "genotype"))
h <- plot(FeatureScatter(subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, idents = levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned))[13], subset = genotype == "KO"), feature1 = "GNLY", feature2 = "KLF2", group.by = "genotype"))
pdf("FS_CD4_genes.pdf", width = , height = 12)
plot((c | d) / (e | f) + (g | h))
dev.off()
# Computing the % of double positive cells (0.5 threshold)
length(WhichCells(subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, subset = genotype == "Mock"), idents = levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned))[3], expression = `Gueguen_CD8_FCGR3A_fix` > 0.5 & `CD8_stem_sc_fix` > 0.5)) / length(WhichCells(subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, subset = genotype == "Mock", idents = levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned))[3]))) * 100
length(WhichCells(subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, subset = genotype == "KO"), idents = levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned))[3], expression = `Gueguen_CD8_FCGR3A_fix` > 0.5 & `CD8_stem_sc_fix` > 0.5)) / length(WhichCells(subset(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned, subset = genotype == "KO", idents = levels(Idents(suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned))[3]))) * 100
suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned@meta.data$CD8_stem_sc_fix <- suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned@meta.data$CD8_stem_sc
suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned@meta.data$Gueguen_CD8_FCGR3A_fix <- suv.car.t.integrated.merged.timepoints.premerged.cd4s.cleaned@meta.data$Gueguen_CD8_FCGR3A
# Changing genotype
suv.car.t.integrated.merged.timepoints.premerged.cd4s[["genotype"]] <- [email protected]$orig.ident
[email protected]$genotype[which([email protected]$genotype == "Mock.1")] <- "Mock"
[email protected]$genotype[which([email protected]$genotype == "Mock.2")] <- "Mock"
[email protected]$genotype[which([email protected]$genotype == "Mock.3")] <- "Mock"
[email protected]$genotype[which([email protected]$genotype == "Suv.ko.1")] <- "KO"
[email protected]$genotype[which([email protected]$genotype == "Suv.ko.2")] <- "KO"
[email protected]$genotype[which([email protected]$genotype == "Suv.ko.3")] <- "KO"
suv.car.t.integrated.merged.timepoints.premerged.cd4s$genotype <- factor(x = suv.car.t.integrated.merged.timepoints.premerged.cd4s$genotype, levels = c("Mock", "KO"))
# tmp
# RoE heatmap