Dynamic transcriptomic remodeling in grafted human neural progenitor cells uncovers mechanisms for vision preservation in a rat model of retinitis pigmentosa.

1. # 20230922_rat_retina_with_GEO.R
3. #### SETUP ####
4. #### Samples and path variables ####
5. # Rat retina treated with human cells
6. # WT rat retina samples from GEO GSE209872
7. # ALIGNED to RAT ref rnor
8. # samples used:
9. #
10. # rnor_08182_P90_Treated_rat_human_cells_outs
11. # rnor_13933_P60_Treated_rat_cells_outs
12. # rnor_13933_P60_Untreated_rat_cells_outs
13. # rnor_16573_P90_Treated_rat_cells_outs
14. # rnor_16573_P90_Untreated_rat_cells_outs
15. # rnor_18482_P60_Treated_rat_cells_outs
16. # rnor_18482_P60_Untreated_rat_cells_outs
17. # rnor_19328_P60_Treated_rat_cells_outs
18. # rnor_19328_P60_Untreated_rat_cells_outs
19. # rnor_SRR20570827_outs
20. # rnor_SRR20570828_outs
22. date <- "20230922"
23. project <- "Rat_retina_with_GEO"
24. datadir <- "/Users/bells/Library/CloudStorage/Box-Box/20230922_rat_retina_with_GEO"
25. sourcedir1 <- "/Users/bells/Library/CloudStorage/Box-Box/20230317_rnor_outs_copy"
26. sourcedir2 <- "/Users/bells/Library/CloudStorage/Box-Box/20230926_rnor_wang_geo_outs_copy"
28. #### Load required packages ####
30. library(Seurat)
31. library(tidyverse)
32. library(Matrix)
33. library(cowplot)
34. library(viridis)
35. library(patchwork)
36. library(pheatmap)
37. library(irlba)
38. library(beepr)
39. library(SeuratDisk)
40. library(EnhancedVolcano)
41. library(Vennerable)
43. #### Set working dir and save session info ####
44. setwd(datadir)
46. sessionInfo()
48. sink(paste0(date,"_",project,"_devtools_sessionInfo3.txt"))
49. devtools::session_info()
50. sink()
52. sink(paste0(date,"_",project,"_sessionInfo3.txt"))
53. sessionInfo()
54. sink()
56. #### Functions for saving images in hires ####
58. hirestiff <- function(saveas){
59. ggsave(
60. saveas,
61. plot = last_plot(),
62. device = "tiff",
63. scale = 1,
64. dpi = 600,
65. limitsize = TRUE
66. )
67. }
69. lowrestiff <- function(saveas){
70. ggsave(
71. saveas,
72. plot = last_plot(),
73. device = "tiff",
74. scale = 1,
75. dpi = 100,
76. limitsize = TRUE
77. )
78. }
80. hirestiffsquare <- function(saveas){
81. ggsave(
82. saveas,
83. plot = last_plot(),
84. device = "tiff",
85. scale = 1,
86. dpi = 600,
87. limitsize = TRUE,
88. bg = "white",
89. width = 9,
90. height = 9,
91. units = c("in")
92. )
93. }
95. lowrestiffsquare <- function(saveas){
96. ggsave(
97. saveas,
98. plot = last_plot(),
99. device = "tiff",
100. scale = 1,
101. dpi = 100,
102. limitsize = TRUE,
103. bg = "white",
104. width = 9,
105. height = 9,
106. units = c("in")
107. )
108. }
110. savepdf <- function(saveas){
111. ggsave(
112. saveas,
113. plot = last_plot(),
114. device = "pdf",
115. scale = 1,
116. dpi = 600,
117. limitsize = TRUE
118. )
119. }
121. savepdfsquare <- function(saveas){
122. ggsave(
123. saveas,
124. plot = last_plot(),
125. device = "pdf",
126. scale = 1,
127. dpi = 600,
128. width = 9,
129. height = 9,
130. limitsize = TRUE
131. )
132. }
134. #### restore prefixes ####
136. # only save after determining cluster resolution
137. # must enter at minimum the date variable and datadir
139. setwd(datadir)
141. # for sct data:
143. # save(datadir, date, ident_res, meaningful.PCs,
144. # prefixPC, prefixPCres, project, res, sourcedir,
145. # hirestiff, hirestiffsquare,
146. # lowrestiff, lowrestiffsquare,
147. # savepdf, savepdfsquare,
148. # file = paste0(datadir,"/",date,"_",project,"_","prefixes.rdata"))
150. # load(paste0(datadir,"/",date,"_",project,"_","prefixes.rdata"))
152. # for integrated data:
154. # save(datadir, date, ident_res, meaningful.PCs, prefixPC, prefixPCres,
155. # project, res, sourcedir1, sourcedir2,
156. # hirestiff, hirestiffsquare,
157. # lowrestiff, lowrestiffsquare,
158. # savepdf, savepdfsquare,
159. # palopal, sorted_palopal,
160. # file = paste0(datadir,"/",date,"_",project,"_","prefixes.rdata"))
162. project <- "Rat_retina_with_GEO_int"
164. load(paste0(datadir,"/",date,"_",project,"_","prefixes.rdata"))
166. seurat_sct <- readRDS(file = paste0(prefixPCres,"_seurat_after_RNAnorm_CELLTYPES.rds"))
168. #### LOAD DATA ####
169. #### Load data from Cellranger output ####
171. sample_list <- c(
172. "rnor_08182_P90_Treated_rat_human_cells_outs",
173. "rnor_13933_P60_Treated_rat_cells_outs",
174. "rnor_13933_P60_Untreated_rat_cells_outs",
175. "rnor_16573_P90_Treated_rat_cells_outs",
176. "rnor_16573_P90_Untreated_rat_cells_outs",
177. "rnor_18482_P60_Treated_rat_cells_outs",
178. "rnor_18482_P60_Untreated_rat_cells_outs",
179. "rnor_19328_P60_Treated_rat_cells_outs",
180. "rnor_19328_P60_Untreated_rat_cells_outs",
181. "rnor_SRR20570827_outs",
182. "rnor_SRR20570828_outs"
183. )
185. # since we have data from two separate source directories
186. # we make two filelists and pathlists
187. filelist1 <- list.files(path = sourcedir1)
188. filelist1 <- as.data.frame(filelist1, files = list.files(path = sourcedir1))
189. path_list1 <- paste0(sourcedir1,"/",filelist1[,1],"/filtered_feature_bc_matrix")
191. filelist2 <- list.files(path = sourcedir2)
192. filelist2 <- as.data.frame(filelist2, files = list.files(path = sourcedir2))
193. path_list2 <- paste0(sourcedir2,"/",filelist2[,1],"/filtered_feature_bc_matrix")
195. path_list <- c(path_list1,path_list2)
196. path_list
198. # check that all files exist
199. all(file.exists(path_list))
201. # check that the sample list and path list are in the same order
202. for (i in 1:length(sample_list)) {
203. TF <- grepl(sample_list[i],path_list[i])
204. print(paste0(sample_list[i]," ",TF))
205. }
207. # make an empty list to store data into
208. data_list <- vector("list")
210. # make a list of all the 10X data
211. for (i in 1:length(sample_list)) {
212. # read in cellranger output
213. cellranger <- Read10X(data.dir = path_list[i])
214. # append data to the list
215. data_list <- c(data_list,cellranger)
216. }
218. updated_samps <- c("08182.p90.TR",
219. "13933.p60.TR",
220. "13933.p60.UT",
221. "16573.p90.TR",
222. "16573.p90.UT",
223. "18482.p60.TR",
224. "18482.p60.UT",
225. "19328.p60.TR",
226. "19328.p60.UT",
227. "GEO27.p00.WT",
228. "GEO28.p00.WT")
230. names(data_list) <- updated_samps
232. #### Seurat Object ####
234. # Initialize the Seurat object with the raw (non-normalized data).
235. # Keep all genes expressed in >= 1 cell
237. # make an empty list to store data into
238. seurat_list <- vector("list")
240. # make a list of all the 10X data
241. for (i in 1:length(data_list)) {
242. # create seurat object
243. s.data <- CreateSeuratObject(counts = data_list[[i]],
244. project = updated_samps[i],
245. min.cells = 1,
246. min.features = 0)
247. # append data to the list
248. seurat_list <- c(seurat_list,s.data)
249. }
251. names(seurat_list) <- updated_samps
253. seurat_list
255. #### Add metadata ####
257. # Create lists for metadata attributes
258. time_point <- c(
259. "08182.p90.TR" = "p90",
260. "13933.p60.TR" = "p60",
261. "13933.p60.UT" = "p60",
262. "16573.p90.TR" = "p90",
263. "16573.p90.UT" = "p90",
264. "18482.p60.TR" = "p60",
265. "18482.p60.UT" = "p60",
266. "19328.p60.TR" = "p60",
267. "19328.p60.UT" = "p60",
268. "GEO27.p00.WT" = "WT",
269. "GEO28.p00.WT" = "WT"
270. )
272. treatment <- c(
273. "08182.p90.TR" = "Cell_Treated",
274. "13933.p60.TR" = "Cell_Treated",
275. "13933.p60.UT" = "Untreated",
276. "16573.p90.TR" = "Cell_Treated",
277. "16573.p90.UT" = "Untreated",
278. "18482.p60.TR" = "Cell_Treated",
279. "18482.p60.UT" = "Untreated",
280. "19328.p60.TR" = "Cell_Treated",
281. "19328.p60.UT" = "Untreated",
282. "GEO27.p00.WT" = "Untreated",
283. "GEO28.p00.WT" = "Untreated"
284. )
286. batch <- c(
287. "08182.p90.TR" = "08182",
288. "13933.p60.TR" = "13933",
289. "13933.p60.UT" = "13933",
290. "16573.p90.TR" = "16573",
291. "16573.p90.UT" = "16573",
292. "18482.p60.TR" = "18482",
293. "18482.p60.UT" = "18482",
294. "19328.p60.TR" = "19328",
295. "19328.p60.UT" = "19328",
296. "GEO27.p00.WT" = "GSE209872",
297. "GEO28.p00.WT" = "GSE209872"
298. )
300. orig.ident <- c(
301. "08182.p90.TR" = "08182.p90.TR",
302. "13933.p60.TR" = "13933.p60.TR",
303. "13933.p60.UT" = "13933.p60.UT",
304. "16573.p90.TR" = "16573.p90.TR",
305. "16573.p90.UT" = "16573.p90.UT",
306. "18482.p60.TR" = "18482.p60.TR",
307. "18482.p60.UT" = "18482.p60.UT",
308. "19328.p60.TR" = "19328.p60.TR",
309. "19328.p60.UT" = "19328.p60.UT",
310. "GEO27.p00.WT" = "GEO27.WT.UT",
311. "GEO28.p00.WT" = "GEO28.WT.UT"
312. )
314. # Iterate through the Seurat objects and assign metadata
315. for (key in names(time_point)) {
316. seurat_list[[key]][["time_point"]] <- time_point[key]
317. seurat_list[[key]][["treatment"]] <- treatment[key]
318. seurat_list[[key]][["batch"]] <- batch[key]
319. seurat_list[[key]][["time_point_treatment"]] <- paste(time_point[key], treatment[key], sep = "_")
320. seurat_list[[key]][["time_point_batch"]] <- paste(time_point[key], batch[key], sep = "_")
321. seurat_list[[key]][["treatment_batch"]] <- paste(treatment[key], batch[key], sep = "_")
322. seurat_list[[key]][["time_point_treatment_batch"]] <- paste(time_point[key], treatment[key], batch[key], sep = "_")
323. seurat_list[[key]][["orig.ident"]] <- orig.ident[key]
324. }
326. #### Merge all data into a single seurat object ####
328. tenx1 <- merge(seurat_list$"08182.p90.TR",
329. y = c(
330. seurat_list$"13933.p60.TR",
331. seurat_list$"13933.p60.UT",
332. seurat_list$"16573.p90.TR",
333. seurat_list$"16573.p90.UT",
334. seurat_list$"18482.p60.TR",
335. seurat_list$"18482.p60.UT",
336. seurat_list$"19328.p60.TR",
337. seurat_list$"19328.p60.UT",
338. seurat_list$"GEO27.p00.WT",
339. seurat_list$"GEO28.p00.WT"
340. ),
341. add.cell.ids = c(
342. "08182.p90.TR",
343. "13933.p60.TR",
344. "13933.p60.UT",
345. "16573.p90.TR",
346. "16573.p90.UT",
347. "18482.p60.TR",
348. "18482.p60.UT",
349. "19328.p60.TR",
350. "19328.p60.UT",
351. "GEO27.p00.WT",
352. "GEO28.p00.WT"
353. ),
354. project = project)
356. tenx1
358. rm(cellranger,
359. data_list,
360. filelist1,
361. filelist2,
362. s.data,
363. seurat_list,
364. batch,
365. i,
366. key,
367. path_list,
368. path_list1,
369. path_list2,
370. sample_list,
371. TF,
372. time_point,
373. treatment,
374. updated_samps)
376. # save merged seurat object
377. # saveRDS(tenx1, file = paste0("./",date,"_",project,"_merged_seurat_prefilter.rds"))
379. # tenx1 <- read_rds(file = paste0("./",date,"_",project,"_merged_seurat_prefilter.rds"))
381. #### QC FILTERING ####
382. #### QC Setup ####
383. # copy tenx1 to another object data just in case
384. data <- tenx1
386. # check to make sure pattern is correct and not grepping other genes
387. grep("^Mt-",rownames(data),value = T) # mitochondrial genes
388. grep("^Rp[sl]",rownames(data),value = T) # ribosomal genes
389. grep("^Mrp[sl]",rownames(data),value = T) # mitochondrial-ribosomal genes
390. grep("^Hb[abq]",rownames(data),value = T) # hemoglobin genes
392. # Add PercentageFeatureSet of various factors to metadata
393. data[["percent.mt"]] <- PercentageFeatureSet(data, pattern = "^Mt-")
394. data[["percent.ribo"]] <- PercentageFeatureSet(data, pattern = "^Rp[sl]")
395. data[["percent.mt.ribo"]] <- PercentageFeatureSet(data, pattern = "^Mrp[sl]")
396. data[["percent.hb"]] <- PercentageFeatureSet(data, pattern = "^Hb[abq]")
398. #### Prefiltered Plots ####
399. vp <- VlnPlot(data, features = "nCount_RNA", pt.size = 0) +
400. NoLegend() +
401. ggtitle("nCount_RNA prefiltered")
402. vp
403. pdf(paste0("./",date,"_",project,"_prefilter_","nCount_RNA",".pdf"))
404. print(vp)
405. dev.off()
407. vp <- VlnPlot(data, features = "nFeature_RNA", pt.size = 0) +
408. NoLegend() +
409. ggtitle("nFeature_RNA prefiltered")
410. vp
411. pdf(paste0("./",date,"_",project,"_prefilter_","nFeature_RNA",".pdf"))
412. print(vp)
413. dev.off()
415. vp <- VlnPlot(data, features = "percent.mt", pt.size = 0) +
416. NoLegend() +
417. ggtitle("percent.mt prefiltered")
418. vp
419. pdf(paste0("./",date,"_",project,"_prefilter_","percent.mt",".pdf"))
420. print(vp)
421. dev.off()
423. vp <- VlnPlot(data, features = "percent.ribo", pt.size = 0) +
424. NoLegend() +
425. ggtitle("percent.ribo prefiltered")
426. vp
427. pdf(paste0("./",date,"_",project,"_prefilter_","percent.ribo",".pdf"))
428. print(vp)
429. dev.off()
431. vp <- VlnPlot(data, features = "percent.mt.ribo", pt.size = 0) +
432. NoLegend() +
433. ggtitle("percent.mt.ribo prefiltered")
434. vp
435. pdf(paste0("./",date,"_",project,"_prefilter_","percent.mt.ribo",".pdf"))
436. print(vp)
437. dev.off()
439. vp <- VlnPlot(data, features = "percent.hb", pt.size = 0) +
440. NoLegend() +
441. ggtitle("percent.hb prefiltered")
442. vp
443. pdf(paste0("./",date,"_",project,"_prefilter_","percent.hb",".pdf"))
444. print(vp)
445. dev.off()
447. p1 <- FeatureScatter(data, feature1 = "nCount_RNA", feature2 = "percent.mt")
448. p2 <- FeatureScatter(data, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
449. p1 + p2
451. pdf(paste0("./",date,"_",project,"_prefilter_","scatterplots",".pdf"))
452. print(p1 + p2)
453. dev.off()
454. #### Additional Prefiltered QC metrics and plots ####
455. qc.metrics <- as_tibble(data[[c("nCount_RNA",
456. "nFeature_RNA",
457. "percent.mt",
458. "percent.ribo",
459. "percent.mt.ribo",
460. "percent.hb")]],
461. rownames="Cell.Barcode")
463. sp <- qc.metrics %>%
464. arrange(percent.mt) %>%
465. ggplot(aes(nCount_RNA,nFeature_RNA,colour=percent.mt)) +
466. geom_point() +
467. scale_color_gradientn(colors=c("black","blue","green2","red","yellow")) +
468. ggtitle("QC metrics - percent Mitochondial DNA Prefiltered")
469. sp
470. pdf(paste0("./",date,"_",project,"_QC_Metrics_pct.mt_prefilter",".pdf"))
471. print(sp)
472. dev.off()
474. sp <- qc.metrics %>%
475. arrange(percent.ribo) %>%
476. ggplot(aes(nCount_RNA,nFeature_RNA,colour=percent.ribo)) +
477. geom_point() +
478. scale_color_gradientn(colors=c("black","blue","green2","red","yellow")) +
479. ggtitle("QC metrics - percent Ribosomal DNA Prefiltered")
480. sp
481. pdf(paste0("./",date,"_",project,"_QC_Metrics_pct.ribo_prefilter",".pdf"))
482. print(sp)
483. dev.off()
485. sp <- qc.metrics %>%
486. arrange(percent.mt.ribo) %>%
487. ggplot(aes(nCount_RNA,nFeature_RNA,colour=percent.mt.ribo)) +
488. geom_point() +
489. scale_color_gradientn(colors=c("black","blue","green2","red","yellow")) +
490. ggtitle("QC metrics - percent Mito-Ribosomal DNA Prefiltered")
491. sp
492. pdf(paste0("./",date,"_",project,"_QC_Metrics_pct.mt.ribo_prefilter",".pdf"))
493. print(sp)
494. dev.off()
496. sp <- qc.metrics %>%
497. arrange(percent.hb) %>%
498. ggplot(aes(nCount_RNA,nFeature_RNA,colour=percent.hb)) +
499. geom_point() +
500. scale_color_gradientn(colors=c("black","blue","green2","red","yellow")) +
501. ggtitle("QC metrics - percent Hemoglobin DNA Prefiltered")
502. sp
503. pdf(paste0("./",date,"_",project,"_QC_Metrics_pct.hb_prefilter",".pdf"))
504. print(sp)
505. dev.off()
507. sp <- qc.metrics %>%
508. ggplot(aes(percent.mt)) +
509. geom_histogram(binwidth = 0.5, fill="yellow", colour="black") +
510. ggtitle("Distribution of Percentage of reads from Mitochondria Prefiltered")
511. sp
512. pdf(paste0("./",date,"_",project,"_Dist_mito_prefilter",".pdf"))
513. print(sp)
514. dev.off()
516. sp <- qc.metrics %>%
517. ggplot(aes(percent.ribo)) +
518. geom_histogram(binwidth = 0.5, fill="yellow", colour="black") +
519. ggtitle("Distribution of Percentage of reads from Ribosomes Prefiltered")
520. sp
521. pdf(paste0("./",date,"_",project,"_Dist_ribo_prefilter",".pdf"))
522. print(sp)
523. dev.off()
525. sp <- qc.metrics %>%
526. ggplot(aes(percent.mt.ribo)) +
527. geom_histogram(binwidth = 0.05, fill="yellow", colour="black") +
528. ggtitle("Distribution of Percentage of reads from Mito-Ribosomes Prefiltered")
529. sp
530. pdf(paste0("./",date,"_",project,"_Dist_mt_ribo_prefilter",".pdf"))
531. print(sp)
532. dev.off()
534. sp <- qc.metrics %>%
535. ggplot(aes(percent.hb)) +
536. geom_histogram(binwidth = 0.5, fill="yellow", colour="black") +
537. ggtitle("Distribution of Percentage of reads from Hemoglobin Prefiltered")
538. sp
539. pdf(paste0("./",date,"_",project,"_Dist_hb_prefilter",".pdf"))
540. print(sp)
541. dev.off()
543. #### Add Z-scores ####
544. data[["nUMI.z"]] <- scale(data$nCount_RNA)
545. data[["nGene.z"]] <- scale(data$nFeature_RNA)
546. data[["percent.mt.z"]] <- scale(data$percent.mt)
547. data[["percent.ribo.z"]] <- scale(data$percent.ribo)
548. data[["percent.mt.ribo.z"]] <- scale(data$percent.mt.ribo)
549. data[["percent.hb.z"]] <- scale(data$percent.hb)
551. #### Prefiltered QC data ####
552. length([email hidden]$orig.ident)
553. mean([email hidden]$percent.mt)
554. mean([email hidden]$percent.ribo)
555. mean([email hidden]$percent.mt.ribo)
556. mean([email hidden]$percent.hb)
557. mean([email hidden]$nCount_RNA)
558. median([email hidden]$nCount_RNA)
559. mean([email hidden]$nFeature_RNA)
560. median([email hidden]$nFeature_RNA)
561. max([email hidden]$nCount_RNA)
563. # save prefilter QC data
564. sink(paste0("./",date,"_",project,"_prefilter_QC_metrics.txt"))
565. cat("length")
566. length([email hidden]$orig.ident)
567. cat("mean percent mito")
568. mean([email hidden]$percent.mt)
569. cat("mean percent ribo")
570. mean([email hidden]$percent.ribo)
571. cat("mean percent mito-ribo")
572. mean([email hidden]$percent.mt.ribo)
573. cat("mean hemoglobin")
574. mean([email hidden]$percent.hb)
575. cat("mean percent counts")
576. mean([email hidden]$nCount_RNA)
577. cat("median counts")
578. median([email hidden]$nCount_RNA)
579. cat("mean features")
580. mean([email hidden]$nFeature_RNA)
581. cat("median features")
582. median([email hidden]$nFeature_RNA)
583. cat("max counts")
584. max([email hidden]$nCount_RNA)
585. sink()
587. #### Filter cells based on Z-score ####
588. data <- subset(data, subset = percent.mt.z < 3)
589. data <- subset(data, subset = percent.ribo.z < 3)
590. data <- subset(data, subset = percent.mt.ribo.z < 3)
591. data <- subset(data, subset = percent.hb.z < 3)
593. # data <- subset(data, subset = nGene.z < 3)
594. # data <- subset(data, subset = nUMI.z < 3)
596. #### Postfiltered QC data ####
597. length([email hidden]$orig.ident)
598. mean([email hidden]$percent.mt)
599. mean([email hidden]$percent.ribo)
600. mean([email hidden]$percent.hb)
601. mean([email hidden]$nCount_RNA)
602. median([email hidden]$nCount_RNA)
603. mean([email hidden]$nFeature_RNA)
604. median([email hidden]$nFeature_RNA)
605. max([email hidden]$nCount_RNA)
607. # save postfilter QC data
608. sink(paste0("./",date,"_",project,"_postfilter_QC_metrics.txt"))
609. cat("length")
610. length([email hidden]$orig.ident)
611. cat("mean percent mito")
612. mean([email hidden]$percent.mt)
613. cat("mean percent ribo")
614. mean([email hidden]$percent.ribo)
615. cat("mean percent mito-ribo")
616. mean([email hidden]$percent.mt.ribo)
617. cat("mean hemoglobin")
618. mean([email hidden]$percent.hb)
619. cat("mean percent counts")
620. mean([email hidden]$nCount_RNA)
621. cat("median counts")
622. median([email hidden]$nCount_RNA)
623. cat("mean features")
624. mean([email hidden]$nFeature_RNA)
625. cat("median features")
626. median([email hidden]$nFeature_RNA)
627. cat("max counts")
628. max([email hidden]$nCount_RNA)
629. sink()
631. #### Postfiltered Plots ####
632. vp <- VlnPlot(data, features = "nCount_RNA", pt.size = 0, group.by = "orig.ident") +
633. NoLegend() +
634. ggtitle("nCount_RNA filtered")
635. vp
636. pdf(paste0("./",date,"_",project,"_filtered_","nCount_RNA",".pdf"))
637. print(vp)
638. dev.off()
640. vp <- VlnPlot(data, features = "nFeature_RNA", pt.size = 0, group.by = "orig.ident") +
641. NoLegend() +
642. ggtitle("nFeature_RNA filtered")
643. vp
644. pdf(paste0("./",date,"_",project,"_filtered_","nFeature_RNA",".pdf"))
645. print(vp)
646. dev.off()
648. vp <- VlnPlot(data, features = "percent.mt", pt.size = 0, group.by = "orig.ident") +
649. NoLegend() +
650. ggtitle("percent.mt filtered")
651. vp
652. pdf(paste0("./",date,"_",project,"_filtered_","percent.mt",".pdf"))
653. print(vp)
654. dev.off()
656. vp <- VlnPlot(data, features = "percent.ribo", pt.size = 0, group.by = "orig.ident") +
657. NoLegend() +
658. ggtitle("percent.ribo filtered")
659. vp
660. pdf(paste0("./",date,"_",project,"_filtered_","percent.ribo",".pdf"))
661. print(vp)
662. dev.off()
664. vp <- VlnPlot(data, features = "percent.mt.ribo", pt.size = 0) +
665. NoLegend() +
666. ggtitle("percent.mt.ribo filtered")
667. vp
668. pdf(paste0("./",date,"_",project,"_filtered_","percent.mt.ribo",".pdf"))
669. print(vp)
670. dev.off()
672. vp <- VlnPlot(data, features = "percent.hb", pt.size = 0) +
673. NoLegend() +
674. ggtitle("percent.hb filtered")
675. vp
676. pdf(paste0("./",date,"_",project,"_filtered_","percent.hb",".pdf"))
677. print(vp)
678. dev.off()
680. p1 <- FeatureScatter(data, feature1 = "nCount_RNA", feature2 = "percent.mt")
681. p2 <- FeatureScatter(data, feature1 = "nCount_RNA", feature2 = "nFeature_RNA")
682. p1 + p2
684. pdf(paste0("./",date,"_",project,"_filtered_","scatterplots",".pdf"))
685. print(p1 + p2)
686. dev.off()
688. #### Additional postfiltered QC metrics and plots ####
689. qc.metrics <- as_tibble(data[[c("nCount_RNA",
690. "nFeature_RNA",
691. "percent.mt",
692. "percent.ribo",
693. "percent.mt.ribo",
694. "percent.hb")]],
695. rownames="Cell.Barcode")
697. sp <- qc.metrics %>%
698. arrange(percent.mt) %>%
699. ggplot(aes(nCount_RNA,nFeature_RNA,colour=percent.mt)) +
700. geom_point() +
701. scale_color_gradientn(colors=c("black","blue","green2","red","yellow")) +
702. ggtitle("QC metrics - percent Mitochondial DNA Postfilter")
703. sp
704. pdf(paste0("./",date,"_",project,"_QC_Metrics_pct.mt_postfilter",".pdf"))
705. print(sp)
706. dev.off()
708. sp <- qc.metrics %>%
709. arrange(percent.ribo) %>%
710. ggplot(aes(nCount_RNA,nFeature_RNA,colour=percent.ribo)) +
711. geom_point() +
712. scale_color_gradientn(colors=c("black","blue","green2","red","yellow")) +
713. ggtitle("QC metrics - percent Ribosomal DNA Postfilter")
714. sp
715. pdf(paste0("./",date,"_",project,"_QC_Metrics_pct.ribo_postfilter",".pdf"))
716. print(sp)
717. dev.off()
719. sp <- qc.metrics %>%
720. arrange(percent.mt.ribo) %>%
721. ggplot(aes(nCount_RNA,nFeature_RNA,colour=percent.mt.ribo)) +
722. geom_point() +
723. scale_color_gradientn(colors=c("black","blue","green2","red","yellow")) +
724. ggtitle("QC metrics - percent Mito-Ribosomal DNA Postfilter")
725. sp
726. pdf(paste0("./",date,"_",project,"_QC_Metrics_pct.mt.ribo_postfilter",".pdf"))
727. print(sp)
728. dev.off()
730. sp <- qc.metrics %>%
731. arrange(percent.hb) %>%
732. ggplot(aes(nCount_RNA,nFeature_RNA,colour=percent.hb)) +
733. geom_point() +
734. scale_color_gradientn(colors=c("black","blue","green2","red","yellow")) +
735. ggtitle("QC metrics - percent Hemoglobin DNA Postfilter")
736. sp
737. pdf(paste0("./",date,"_",project,"_QC_Metrics_pct.hb_postfilter",".pdf"))
738. print(sp)
739. dev.off()
741. sp <- qc.metrics %>%
742. ggplot(aes(percent.mt)) +
743. geom_histogram(binwidth = 0.25, fill="yellow", colour="black") +
744. ggtitle("Distribution of Percentage of reads from Mitochondria Postfilter")
745. sp
746. pdf(paste0("./",date,"_",project,"_Dist_mito_postfilter",".pdf"))
747. print(sp)
748. dev.off()
750. sp <- qc.metrics %>%
751. ggplot(aes(percent.ribo)) +
752. geom_histogram(binwidth = 0.5, fill="yellow", colour="black") +
753. ggtitle("Distribution of Percentage of reads from Ribosomes Postfilter")
754. sp
755. pdf(paste0("./",date,"_",project,"_Dist_ribo_postfilter",".pdf"))
756. print(sp)
757. dev.off()
759. sp <- qc.metrics %>%
760. ggplot(aes(percent.mt.ribo)) +
761. geom_histogram(binwidth = 0.05, fill="yellow", colour="black") +
762. ggtitle("Distribution of Percentage of reads from Mito-Ribosomes Postfilter")
763. sp
764. pdf(paste0("./",date,"_",project,"_Dist_mt_ribo_postfilter",".pdf"))
765. print(sp)
766. dev.off()
768. sp <- qc.metrics %>%
769. ggplot(aes(percent.hb)) +
770. geom_histogram(binwidth = 0.05, fill="yellow", colour="black") +
771. ggtitle("Distribution of Percentage of reads from Hemoglobin Postfilter")
772. sp
773. pdf(paste0("./",date,"_",project,"_Dist_hb_postfilter",".pdf"))
774. print(sp)
775. dev.off()
777. #### post QC cleanup ####
779. tenx1
780. data
782. tenx1 <- data
784. rm(vp,sp,p1,p2,
785. qc.metrics,data)
787. #### set factor levels ####
788. seurat_sct <- tenx1
790. seurat_sct$orig.ident <- factor(
791. x = seurat_sct$orig.ident,
792. levels = c(
793. "GEO27.WT.UT",
794. "GEO28.WT.UT",
795. "08182.p90.TR",
796. "13933.p60.TR",
797. "13933.p60.UT",
798. "16573.p90.TR",
799. "16573.p90.UT",
800. "18482.p60.TR",
801. "18482.p60.UT",
802. "19328.p60.TR",
803. "19328.p60.UT"
804. )
805. )
807. seurat_sct$treatment <- factor(
808. x = seurat_sct$treatment,
809. levels = c(
810. "Untreated",
811. "Cell_Treated"
812. )
813. )
815. seurat_sct$batch <- factor(
816. x = seurat_sct$batch,
817. levels = c(
818. "GSE209872",
819. "08182",
820. "13933",
821. "16573",
822. "18482",
823. "19328"
824. )
825. )
827. seurat_sct$time_point <- factor(
828. x = seurat_sct$time_point,
829. levels = c(
830. "WT",
831. "p60",
832. "p90"
833. )
834. )
836. seurat_sct$time_point_treatment <- factor(
837. x = seurat_sct$time_point_treatment,
838. levels = c(
839. "WT_Untreated",
840. "p60_Cell_Treated",
841. "p60_Untreated",
842. "p90_Cell_Treated",
843. "p90_Untreated"
844. )
845. )
847. seurat_sct$time_point_batch <- factor(
848. x = seurat_sct$time_point_batch,
849. levels = c(
850. "WT_GSE209872",
851. "p60_13933",
852. "p60_18482",
853. "p60_19328",
854. "p90_08182",
855. "p90_16573"
856. )
857. )
859. seurat_sct$treatment_batch <- factor(
860. x = seurat_sct$treatment_batch,
861. levels = c(
862. "Untreated_GSE209872",
863. "Untreated_13933",
864. "Untreated_18482",
865. "Untreated_19328",
866. "Untreated_16573",
867. "Cell_Treated_13933",
868. "Cell_Treated_18482",
869. "Cell_Treated_19328",
870. "Cell_Treated_08182",
871. "Cell_Treated_16573"
872. )
873. )
875. seurat_sct$time_point_treatment_batch <- factor(
876. x = seurat_sct$time_point_treatment_batch,
877. levels = c(
878. "WT_Untreated_GSE209872",
879. "p60_Untreated_13933",
880. "p60_Untreated_18482",
881. "p60_Untreated_19328",
882. "p90_Untreated_16573",
883. "p60_Cell_Treated_13933",
884. "p60_Cell_Treated_18482",
885. "p60_Cell_Treated_19328",
886. "p90_Cell_Treated_08182",
887. "p90_Cell_Treated_16573"
888. )
889. )
891. rm(tenx1)
893. #### SCTransform ####
894. seurat_sct
896. seurat_sct <- SCTransform(seurat_sct,
897. ncells = 11848, #total cells / 10 (118478/10 = 11848)
898. vars.to.regress = c("percent.mt","percent.ribo"),
899. verbose = TRUE)
901. #### PCA ####
902. # this performs PCA on the seurat object
903. seurat_sct <- RunPCA(seurat_sct, npcs = 50, verbose = TRUE)
905. # make PC coordinate object a data frame
906. xx.coord <- as.data.frame(seurat_sct@reductions$[email hidden])
908. # make PC feature loadings object a data frame
909. xx.gload <- as.data.frame(seurat_sct@reductions$[email hidden])
911. # calculate eigenvalues for arrays
913. # generate squares of all sample coordinates
914. sq.xx.coord <- as.data.frame(xx.coord^2)
915. # create empty list for eigenvalues first
916. eig <- c()
917. # calculate the eigenvalue for each PC in sq.xx.coord by taking the sqrt of the sum of squares
918. for(i in 1:ncol(sq.xx.coord))
919. eig[i] = sqrt(sum(sq.xx.coord[,i]))
920. # calculate the total variance by adding up all the eigenvalues
921. sum.eig <- sum(eig)
922. # calculate the expected contribution of all PCs if they all contribute equally to the total variance
923. expected.contribution <- sum.eig/(length(xx.coord)-1)
924. # return the number of principal components with an eigenvalue greater than expected by equal variance
925. meaningful.PCs <- sum(eig > expected.contribution)
927. # create empty list for eigenvalue percentage
928. eig.percent <- c()
929. # calculate the percentage of the total variance by each PC eigenvalue
930. for(i in 1:length(eig))
931. eig.percent[i] = 100*eig[i]/sum.eig
932. # sum of all eig.percent should total to 100
933. sum(eig.percent)
934. # create empty list for scree values
935. scree <- c()
936. # calculate a running total of variance contribution
937. for(i in 1:length(eig))
938. if(i == 1) scree[i] = eig.percent[i] else scree[i] = scree[i-1] + eig.percent[i]
940. # create data frame for eigenvalue summaries
941. eigenvalues <- data.frame("PC" = colnames(xx.coord), "eig" = eig, "percent" = eig.percent, "scree" = scree)
943. # write csv for eigenvalues
944. # write.csv(eigenvalues, file = paste0("./",date,"_",project,"_PCA_eigenvalues.csv"), row.names = F)
946. # plot scree values
947. plot(eigenvalues$percent, ylim = c(0,100), type = "S", xlab = "PC", ylab = "Percent of variance",
948. main = paste0(date,"_",project," scree plot all samples PCA"))
949. points(eigenvalues$scree, ylim = c(0,100), type = "p", pch = 16)
950. lines(eigenvalues$scree)
951. # add red line to indicate cut-off
952. cut.off <- 100/(length(eig)-1)
953. abline(h = cut.off, col = "red")
954. # add blue line to indicate which PCs are meaningful and kept
955. abline(v = meaningful.PCs, col = "blue")
956. text(meaningful.PCs, cut.off, label = paste("cutoff PC",meaningful.PCs),
957. adj = c(-0.1, -0.5))
959. dev.copy(pdf, paste0("./",date,"_",project,"_scree_plot.pdf"))
960. dev.off()
962. rm(eigenvalues,sq.xx.coord,xx.coord,xx.gload,cut.off,
963. eig,eig.percent,expected.contribution,i,scree,sum.eig)
965. # meaningful.PCs <- 15
967. #### Run UMAP and look at UMAP plots ####
969. seurat_sct <- RunUMAP(seurat_sct, reduction = "pca", dims = 1:meaningful.PCs, verbose = TRUE)
971. # update prefixed variable
972. prefixPC <- paste0("./",date,"_",project,"_",meaningful.PCs,"PCs")
974. ## UMAP plot by sample name ("orig.ident")
975. DimPlot(seurat_sct,
976. reduction = "umap",
977. label = FALSE,
978. pt.size = .25,
979. group.by = "orig.ident",
980. raster = F
981. )
982. hirestiff(paste0(prefixPC,"_UMAP_by_sample_hires.tiff"))
983. lowrestiff(paste0(prefixPC,"_UMAP_by_sample_lowres.tiff"))
985. DimPlot(seurat_sct,
986. reduction = "umap",
987. label = FALSE,
988. pt.size = .25,
989. group.by = "treatment",
990. raster = F
991. )
992. hirestiff(paste0(prefixPC,"_UMAP_by_treatment_hires.tiff"))
993. lowrestiff(paste0(prefixPC,"_UMAP_by_treatment_lowres.tiff"))
995. DimPlot(seurat_sct,
996. reduction = "umap",
997. label = FALSE,
998. pt.size = .25,
999. split.by = "treatment",
1000. group.by = "treatment",
1001. raster = F
1002. )
1003. hirestiff(paste0(prefixPC,"_UMAP_split_by_treatment_hires.tiff"))
1004. lowrestiff(paste0(prefixPC,"_UMAP_split_by_treatment_lowres.tiff"))
1006. DimPlot(seurat_sct,
1007. reduction = "umap",
1008. label = FALSE,
1009. pt.size = .25,
1010. group.by = "batch",
1011. raster = F
1012. )
1013. hirestiff(paste0(prefixPC,"_UMAP_by_batch_hires.tiff"))
1014. lowrestiff(paste0(prefixPC,"_UMAP_by_batch_lowres.tiff"))
1016. DimPlot(seurat_sct,
1017. reduction = "umap",
1018. label = FALSE,
1019. pt.size = .25,
1020. group.by = "batch",
1021. split.by = "batch",
1022. raster = F
1023. )
1024. hirestiff(paste0(prefixPC,"_UMAP_split_by_batch_hires.tiff"))
1025. lowrestiff(paste0(prefixPC,"_UMAP_split_by_batch_lowres.tiff"))
1027. DimPlot(seurat_sct,
1028. reduction = "umap",
1029. label = FALSE,
1030. pt.size = .25,
1031. group.by = "time_point",
1032. raster = F
1033. )
1034. hirestiff(paste0(prefixPC,"_UMAP_by_time_point_hires.tiff"))
1035. lowrestiff(paste0(prefixPC,"_UMAP_by_time_point_lowres.tiff"))
1037. DimPlot(seurat_sct,
1038. reduction = "umap",
1039. label = FALSE,
1040. pt.size = .25,
1041. split.by = "time_point",
1042. group.by = "time_point",
1043. raster = F
1044. )
1045. hirestiff(paste0(prefixPC,"_UMAP_split_by_time_point_hires.tiff"))
1046. lowrestiff(paste0(prefixPC,"_UMAP_split_by_time_point_lowres.tiff"))
1048. # saveRDS(seurat_sct, file = paste0(prefixPC,"_seurat_integrated_preclustering.rds"))
1050. # seurat_sct <- read_rds(file = paste0(prefixPC,"_seurat_integrated_preclustering.rds"))
1051. #### (skip to integration) ####
1052. #### SCT Clustering and Resolution ####
1053. # DefaultAssay(seurat_sct) <- "SCT"
1055. # Determine the K-nearest neighbor graph
1056. seurat_sct <- FindNeighbors(object = seurat_sct, reduction = "pca", dims = 1:meaningful.PCs)
1058. # Determine the clusters
1060. seurat_sct <- FindClusters(object = seurat_sct,
1061. resolution = c(0.1,0.2,0.3,0.4,0.5))
1063. res <- "_res.0.1"
1064. ident_res <- paste0("SCT_snn",res)
1066. Idents(seurat_sct) <- ident_res
1067. DimPlot(seurat_sct, reduction = "umap", label = TRUE, label.size = 5, pt.size = 0.8) +
1068. NoLegend() +
1069. ggtitle(paste0(ident_res))
1071. #update prefix
1072. prefixPCres <- paste0(prefixPC,res)
1074. # Plot the UMAP
1075. DimPlot(seurat_sct, reduction = "umap", label = TRUE, label.size = 5, pt.size = 0.8) +
1076. NoLegend() +
1077. ggtitle(paste0(ident_res))
1078. hirestiff(paste0(prefixPCres,"_UMAP","_by_","cluster","_hires.tiff"))
1079. lowrestiff(paste0(prefixPCres,"_UMAP","_by_","cluster","_lowres.tiff"))
1081. # UMAP of cells in each cluster by treatment without cluster labels
1082. DimPlot(seurat_sct, reduction = "umap", label = FALSE, split.by = "treatment", pt.size = 0.8) + NoLegend()
1083. hirestiff(paste0(prefixPCres,"_UMAP","_by_","treatment","_no_labels","_hires.tiff"))
1084. lowrestiff(paste0(prefixPCres,"_UMAP","_by_","treatment","_no_labels","_lowres.tiff"))
1086. # UMAP of cells in each cluster by treatment with cluster labels
1087. DimPlot(seurat_sct, reduction = "umap", label = TRUE,
1088. split.by = "treatment", pt.size = 0.8, label.size = 5) + NoLegend()
1089. hirestiff(paste0(prefixPCres,"_UMAP","_by_","treatment","_with_clusters_labels","_hires.tiff"))
1090. lowrestiff(paste0(prefixPCres,"_UMAP","_by_","treatment","_with_clusters_labels","_lowres.tiff"))
1092. # UMAP of cells in each cluster by batch without cluster labels
1093. DimPlot(seurat_sct, reduction = "umap", label = FALSE, split.by = "batch", pt.size = 0.8) + NoLegend()
1094. hirestiff(paste0(prefixPCres,"_UMAP","_by_","batch","_no_labels","_hires.tiff"))
1095. lowrestiff(paste0(prefixPCres,"_UMAP","_by_","batch","_no_labels","_lowres.tiff"))
1097. # UMAP of cells in each cluster by batch with cluster labels
1098. DimPlot(seurat_sct, reduction = "umap", label = TRUE,
1099. split.by = "batch", pt.size = 0.8, label.size = 5) + NoLegend()
1100. hirestiff(paste0(prefixPCres,"_UMAP","_by_","batch","_with_clusters_labels","_hires.tiff"))
1101. lowrestiff(paste0(prefixPCres,"_UMAP","_by_","batch","_with_clusters_labels","_lowres.tiff"))
1103. # UMAP of cells in each cluster by time_point without cluster labels
1104. DimPlot(seurat_sct, reduction = "umap", label = FALSE, split.by = "time_point", pt.size = 0.8) + NoLegend()
1105. hirestiff(paste0(prefixPCres,"_UMAP","_by_","time_point","_no_labels","_hires.tiff"))
1106. lowrestiff(paste0(prefixPCres,"_UMAP","_by_","time_point","_no_labels","_lowres.tiff"))
1108. # UMAP of cells in each cluster by time_point with cluster labels
1109. DimPlot(seurat_sct, reduction = "umap", label = TRUE,
1110. split.by = "time_point", pt.size = 0.8, label.size = 5) + NoLegend()
1111. hirestiff(paste0(prefixPCres,"_UMAP","_by_","time_point","_with_clusters_labels","_hires.tiff"))
1112. lowrestiff(paste0(prefixPCres,"_UMAP","_by_","time_point","_with_clusters_labels","_lowres.tiff"))
1114. #### Extract number of cells per cluster per orig.ident ####
1115. n_cells <- FetchData(seurat_sct, vars = c("ident", "orig.ident")) %>%
1116. dplyr::count(ident, orig.ident) %>%
1117. tidyr::spread(ident, n)
1119. write.csv(n_cells, file = paste0(prefixPCres,"_cells_per_cluster.csv"))
1121. rm(n_cells)
1123. #### save RDS containing reduction and cluster idents ####
1124. saveRDS(seurat_sct, paste0(prefixPCres,"_seurat_after_clustering.rds"))
1126. # seurat_sct <- readRDS(file = paste0(prefixPCres,"_seurat_after_clustering.rds"))
1128. #### normalize rna slot ####
1130. # Select the RNA counts slot to be the default assay for visualization purposes
1131. DefaultAssay(seurat_sct) <- "RNA"
1133. # Normalize, find variable features, scale data
1134. seurat_sct <- NormalizeData(seurat_sct)
1135. seurat_sct <- FindVariableFeatures(seurat_sct)
1136. all.genes <- rownames(seurat_sct)
1137. seurat_sct <- ScaleData(seurat_sct, features = all.genes)
1139. # save object containing RNA normalized data
1140. saveRDS(seurat_sct, paste0(prefixPCres,"_seurat_after_RNAnorm.rds"))
1142. seurat_sct <- readRDS(file = paste0(prefixPCres,"_seurat_after_RNAnorm.rds"))
1144. # SaveH5Seurat(seurat_sct,
1145. # filename = paste0(prefixPCres,"_seurat_after_RNAnorm.h5Seurat"))
1146. #
1147. # seurat_sct_test <- LoadH5Seurat(file = paste0(prefixPCres,"_seurat_after_RNAnorm.h5Seurat"))
1149. #### ####
1150. #### Run integration (RPCA method) ####
1151. DefaultAssay(seurat_sct) <- "RNA"
1153. s.list <- SplitObject(seurat_sct, split.by = "orig.ident")
1155. # normalize data and find variable features
1156. s.list <- lapply(X = s.list, FUN = function(x) {
1157. x <- NormalizeData(x, verbose = FALSE)
1158. x <- FindVariableFeatures(x, verbose = FALSE)
1159. })
1161. # Next, select features for downstream integration, and run PCA on each
1162. # object in the list, which is required for running the alternative
1163. # reciprocal PCA workflow.
1165. features <- SelectIntegrationFeatures(object.list = s.list)
1167. s.list <- lapply(X = s.list, FUN = function(x) {
1168. x <- ScaleData(x, features = features, verbose = FALSE)
1169. x <- RunPCA(x, features = features, verbose = FALSE)
1170. })
1172. anchors <- FindIntegrationAnchors(object.list = s.list,
1173. reduction = "rpca",
1174. k.anchor = 10) # default is 5
1176. s.integrated <- IntegrateData(anchorset = anchors, dims = 1:50)
1177. s.integrated <- ScaleData(s.integrated, verbose = FALSE)
1179. # saveRDS(s.integrated, file = paste0("./",date,"_",project,"_seurat_integrated.rds"))
1181. # s.integrated <- read_rds(file = paste0("./",date,"_",project,"_seurat_integrated.rds"))
1183. seurat_sct <- s.integrated
1185. rm(s.list,anchors,s.integrated)
1187. #### set factor levels ####
1188. seurat_sct$orig.ident <- factor(
1189. x = seurat_sct$orig.ident,
1190. levels = c(
1191. "GEO27.WT.UT",
1192. "GEO28.WT.UT",
1193. "08182.p90.TR",
1194. "13933.p60.TR",
1195. "13933.p60.UT",
1196. "16573.p90.TR",
1197. "16573.p90.UT",
1198. "18482.p60.TR",
1199. "18482.p60.UT",
1200. "19328.p60.TR",
1201. "19328.p60.UT"
1202. )
1203. )
1205. seurat_sct$treatment <- factor(
1206. x = seurat_sct$treatment,
1207. levels = c(
1208. "Untreated",
1209. "Cell_Treated"
1210. )
1211. )
1213. seurat_sct$batch <- factor(
1214. x = seurat_sct$batch,
1215. levels = c(
1216. "GSE209872",
1217. "08182",
1218. "13933",
1219. "16573",
1220. "18482",
1221. "19328"
1222. )
1223. )
1225. seurat_sct$time_point <- factor(
1226. x = seurat_sct$time_point,
1227. levels = c(
1228. "WT",
1229. "p60",
1230. "p90"
1231. )
1232. )
1234. seurat_sct$time_point_treatment <- factor(
1235. x = seurat_sct$time_point_treatment,
1236. levels = c(
1237. "WT_Untreated",
1238. "p60_Cell_Treated",
1239. "p60_Untreated",
1240. "p90_Cell_Treated",
1241. "p90_Untreated"
1242. )
1243. )
1245. seurat_sct$time_point_batch <- factor(
1246. x = seurat_sct$time_point_batch,
1247. levels = c(
1248. "WT_GSE209872",
1249. "p60_13933",
1250. "p60_18482",
1251. "p60_19328",
1252. "p90_08182",
1253. "p90_16573"
1254. )
1255. )
1257. seurat_sct$treatment_batch <- factor(
1258. x = seurat_sct$treatment_batch,
1259. levels = c(
1260. "Untreated_GSE209872",
1261. "Untreated_13933",
1262. "Untreated_18482",
1263. "Untreated_19328",
1264. "Untreated_16573",
1265. "Cell_Treated_13933",
1266. "Cell_Treated_18482",
1267. "Cell_Treated_19328",
1268. "Cell_Treated_08182",
1269. "Cell_Treated_16573"
1270. )
1271. )
1273. seurat_sct$time_point_treatment_batch <- factor(
1274. x = seurat_sct$time_point_treatment_batch,
1275. levels = c(
1276. "WT_Untreated_GSE209872",
1277. "p60_Untreated_13933",
1278. "p60_Untreated_18482",
1279. "p60_Untreated_19328",
1280. "p90_Untreated_16573",
1281. "p60_Cell_Treated_13933",
1282. "p60_Cell_Treated_18482",
1283. "p60_Cell_Treated_19328",
1284. "p90_Cell_Treated_08182",
1285. "p90_Cell_Treated_16573"
1286. )
1287. )
1289. #### update project variable to include integration ####
1290. project <- paste0(project,"_int")
1292. #### INT PCA ####
1293. # this performs PCA on the seurat object
1294. seurat_sct <- RunPCA(seurat_sct, npcs = 50, verbose = TRUE)
1296. # make PC coordinate object a data frame
1297. xx.coord <- as.data.frame(seurat_sct@reductions$[email hidden])
1299. # make PC feature loadings object a data frame
1300. xx.gload <- as.data.frame(seurat_sct@reductions$[email hidden])
1302. # calculate eigenvalues for arrays
1304. # generate squares of all sample coordinates
1305. sq.xx.coord <- as.data.frame(xx.coord^2)
1306. # create empty list for eigenvalues first
1307. eig <- c()
1308. # calculate the eigenvalue for each PC in sq.xx.coord by taking the sqrt of the sum of squares
1309. for(i in 1:ncol(sq.xx.coord))
1310. eig[i] = sqrt(sum(sq.xx.coord[,i]))
1311. # calculate the total variance by adding up all the eigenvalues
1312. sum.eig <- sum(eig)
1313. # calculate the expected contribution of all PCs if they all contribute equally to the total variance
1314. expected.contribution <- sum.eig/(length(xx.coord)-1)
1315. # return the number of principal components with an eigenvalue greater than expected by equal variance
1316. meaningful.PCs <- sum(eig > expected.contribution)
1318. # create empty list for eigenvalue percentage
1319. eig.percent <- c()
1320. # calculate the percentage of the total variance by each PC eigenvalue
1321. for(i in 1:length(eig))
1322. eig.percent[i] = 100*eig[i]/sum.eig
1323. # sum of all eig.percent should total to 100
1324. sum(eig.percent)
1325. # create empty list for scree values
1326. scree <- c()
1327. # calculate a running total of variance contribution
1328. for(i in 1:length(eig))
1329. if(i == 1) scree[i] = eig.percent[i] else scree[i] = scree[i-1] + eig.percent[i]
1331. # create data frame for eigenvalue summaries
1332. eigenvalues <- data.frame("PC" = colnames(xx.coord), "eig" = eig, "percent" = eig.percent, "scree" = scree)
1334. # write csv for eigenvalues
1335. # write.csv(eigenvalues, file = paste0("./",date,"_",project,"_PCA_eigenvalues.csv"), row.names = F)
1337. # plot scree values
1338. plot(eigenvalues$percent, ylim = c(0,100), type = "S", xlab = "PC", ylab = "Percent of variance",
1339. main = paste0(date,"_",project," scree plot all samples PCA"))
1340. points(eigenvalues$scree, ylim = c(0,100), type = "p", pch = 16)
1341. lines(eigenvalues$scree)
1342. # add red line to indicate cut-off
1343. cut.off <- 100/(length(eig)-1)
1344. abline(h = cut.off, col = "red")
1345. # add blue line to indicate which PCs are meaningful and kept
1346. abline(v = meaningful.PCs, col = "blue")
1347. text(meaningful.PCs, cut.off, label = paste("cutoff PC",meaningful.PCs),
1348. adj = c(-0.1, -0.5))
1350. dev.copy(pdf, paste0("./",date,"_",project,"_scree_plot.pdf"))
1351. dev.off()
1353. rm(eigenvalues,sq.xx.coord,xx.coord,xx.gload,cut.off,
1354. eig,eig.percent,expected.contribution,i,scree,sum.eig)
1356. # meaningful.PCs <- 12
1358. #### INT Run UMAP and look at UMAP plots ####
1360. seurat_sct <- RunUMAP(seurat_sct,
1361. reduction = "pca",
1362. dims = 1:meaningful.PCs,
1363. verbose = TRUE)
1365. # update prefixed variable
1366. prefixPC <- paste0("./",date,"_",project,"_",meaningful.PCs,"PCs")
1368. ## UMAP plot by sample name ("orig.ident")
1369. DimPlot(seurat_sct,
1370. reduction = "umap",
1371. label = FALSE,
1372. pt.size = .25,
1373. group.by = "orig.ident",
1374. raster = F
1375. )
1376. hirestiff(paste0(prefixPC,"_UMAP_by_sample_hires.tiff"))
1377. lowrestiff(paste0(prefixPC,"_UMAP_by_sample_lowres.tiff"))
1379. DimPlot(seurat_sct,
1380. reduction = "umap",
1381. label = FALSE,
1382. pt.size = .25,
1383. group.by = "treatment",
1384. raster = F
1385. )
1386. hirestiff(paste0(prefixPC,"_UMAP_by_treatment_hires.tiff"))
1387. lowrestiff(paste0(prefixPC,"_UMAP_by_treatment_lowres.tiff"))
1389. DimPlot(seurat_sct,
1390. reduction = "umap",
1391. label = FALSE,
1392. pt.size = .25,
1393. split.by = "treatment",
1394. group.by = "treatment",
1395. raster = F
1396. )
1397. hirestiff(paste0(prefixPC,"_UMAP_split_by_treatment_hires.tiff"))
1398. lowrestiff(paste0(prefixPC,"_UMAP_split_by_treatment_lowres.tiff"))
1400. DimPlot(seurat_sct,
1401. reduction = "umap",
1402. label = FALSE,
1403. pt.size = .25,
1404. group.by = "batch",
1405. raster = F
1406. )
1407. hirestiff(paste0(prefixPC,"_UMAP_by_batch_hires.tiff"))
1408. lowrestiff(paste0(prefixPC,"_UMAP_by_batch_lowres.tiff"))
1410. DimPlot(seurat_sct,
1411. reduction = "umap",
1412. label = FALSE,
1413. pt.size = .25,
1414. group.by = "batch",
1415. split.by = "batch",
1416. raster = F
1417. )
1418. hirestiff(paste0(prefixPC,"_UMAP_split_by_batch_hires.tiff"))
1419. lowrestiff(paste0(prefixPC,"_UMAP_split_by_batch_lowres.tiff"))
1421. DimPlot(seurat_sct,
1422. reduction = "umap",
1423. label = FALSE,
1424. pt.size = .25,
1425. group.by = "time_point",
1426. raster = F
1427. )
1428. hirestiff(paste0(prefixPC,"_UMAP_by_time_point_hires.tiff"))
1429. lowrestiff(paste0(prefixPC,"_UMAP_by_time_point_lowres.tiff"))
1431. DimPlot(seurat_sct,
1432. reduction = "umap",
1433. label = FALSE,
1434. pt.size = .25,
1435. split.by = "time_point",
1436. group.by = "time_point",
1437. raster = F
1438. )
1439. hirestiff(paste0(prefixPC,"_UMAP_split_by_time_point_hires.tiff"))
1440. lowrestiff(paste0(prefixPC,"_UMAP_split_by_time_point_lowres.tiff"))
1442. DimPlot(seurat_sct,
1443. reduction = "umap",
1444. label = FALSE,
1445. pt.size = .25,
1446. group.by = "time_point_treatment",
1447. raster = F
1448. )
1449. hirestiff(paste0(prefixPC,"_UMAP_by_time_point_treatment_hires.tiff"))
1450. lowrestiff(paste0(prefixPC,"_UMAP_by_time_point_treatment_lowres.tiff"))
1452. DimPlot(seurat_sct,
1453. reduction = "umap",
1454. label = FALSE,
1455. pt.size = .25,
1456. split.by = "time_point_treatment",
1457. group.by = "time_point_treatment",
1458. raster = F
1459. )
1460. hirestiff(paste0(prefixPC,"_UMAP_split_by_time_point_treatment_hires.tiff"))
1461. lowrestiff(paste0(prefixPC,"_UMAP_split_by_time_point_treatment_lowres.tiff"))
1463. DimPlot(seurat_sct,
1464. reduction = "umap",
1465. label = FALSE,
1466. pt.size = .25,
1467. group.by = "time_point_batch",
1468. raster = F
1469. )
1470. hirestiff(paste0(prefixPC,"_UMAP_by_time_point_batch_hires.tiff"))
1471. lowrestiff(paste0(prefixPC,"_UMAP_by_time_point_batch_lowres.tiff"))
1473. DimPlot(seurat_sct,
1474. reduction = "umap",
1475. label = FALSE,
1476. pt.size = .25,
1477. split.by = "time_point_batch",
1478. group.by = "time_point_batch",
1479. raster = F
1480. )
1481. hirestiff(paste0(prefixPC,"_UMAP_split_by_time_point_batch_hires.tiff"))
1482. lowrestiff(paste0(prefixPC,"_UMAP_split_by_time_point_batch_lowres.tiff"))
1484. DimPlot(seurat_sct,
1485. reduction = "umap",
1486. label = FALSE,
1487. pt.size = .25,
1488. group.by = "treatment_batch",
1489. raster = F
1490. )
1491. hirestiff(paste0(prefixPC,"_UMAP_by_treatment_batch_hires.tiff"))
1492. lowrestiff(paste0(prefixPC,"_UMAP_by_treatment_batch_lowres.tiff"))
1494. DimPlot(seurat_sct,
1495. reduction = "umap",
1496. label = FALSE,
1497. pt.size = .25,
1498. split.by = "treatment_batch",
1499. group.by = "treatment_batch",
1500. raster = F
1501. )
1502. hirestiff(paste0(prefixPC,"_UMAP_split_by_treatment_batch_hires.tiff"))
1503. lowrestiff(paste0(prefixPC,"_UMAP_split_by_treatment_batch_lowres.tiff"))
1505. DimPlot(seurat_sct,
1506. reduction = "umap",
1507. label = FALSE,
1508. pt.size = .25,
1509. group.by = "time_point_treatment_batch",
1510. raster = F
1511. )
1512. hirestiff(paste0(prefixPC,"_UMAP_by_time_point_treatment_batch_hires.tiff"))
1513. lowrestiff(paste0(prefixPC,"_UMAP_by_time_point_treatment_batch_lowres.tiff"))
1515. DimPlot(seurat_sct,
1516. reduction = "umap",
1517. label = FALSE,
1518. pt.size = .25,
1519. split.by = "time_point_treatment_batch",
1520. group.by = "time_point_treatment_batch",
1521. raster = F
1522. )
1523. hirestiff(paste0(prefixPC,"_UMAP_split_by_time_point_treatment_batch_hires.tiff"))
1524. lowrestiff(paste0(prefixPC,"_UMAP_split_by_time_point_treatment_batch_lowres.tiff"))
1526. # saveRDS(seurat_sct, file = paste0(prefixPC,"_seurat_integrated_preclustering.rds"))
1528. # seurat_sct <- read_rds(file = paste0(prefixPC,"_seurat_integrated_preclustering.rds"))
1530. #### INT Clustering and Resolution ####
1531. # DefaultAssay(seurat_sct) <- "integrated"
1533. # Determine the K-nearest neighbor graph
1534. seurat_sct <- FindNeighbors(object = seurat_sct,
1535. reduction = "pca",
1536. dims = 1:meaningful.PCs)
1538. # Determine the clusters
1539. # seurat_sct <- FindClusters(object = seurat_sct,
1540. # resolution = c(0.1,0.2,0.3,0.4,0.5))
1542. seurat_sct <- FindClusters(object = seurat_sct,
1543. resolution = c(0.5))
1545. res <- "_res.0.5"
1546. ident_res <- paste0("integrated_snn",res)
1548. Idents(seurat_sct) <- ident_res
1550. DimPlot(seurat_sct,
1551. reduction = "umap",
1552. label = TRUE,
1553. label.size = 5,
1554. pt.size = 0.8,
1555. raster = F
1556. ) +
1557. NoLegend() +
1558. ggtitle(paste0(ident_res))
1560. # set cluster colors with Palo (makes sure adjacent colors are different)
1561. palopal <- Palo::Palo(seurat_sct[["umap"]]@cell.embeddings,
1562. as.character(Idents(seurat_sct)),
1563. scales::hue_pal()(length(levels(seurat_sct))))
1565. # sort these so they are in numeric order and can be transferred to new idents
1566. sorted_palopal <- palopal[order(as.numeric(names(palopal)))]
1568. # test resolution with Palo colors
1569. DimPlot(seurat_sct,
1570. reduction = "umap",
1571. label = TRUE,
1572. label.size = 5,
1573. pt.size = 0.8,
1574. cols = sorted_palopal,
1575. raster = F
1576. ) +
1577. NoLegend() +
1578. ggtitle(paste0(ident_res))
1580. #update prefix
1581. prefixPCres <- paste0(prefixPC,res)
1583. # Plot the UMAP
1584. DimPlot(seurat_sct,
1585. reduction = "umap",
1586. label = TRUE,
1587. label.size = 5,
1588. pt.size = 0.8,
1589. cols = sorted_palopal,
1590. raster = F
1591. ) +
1592. NoLegend() +
1593. ggtitle(paste0(ident_res))
1594. hirestiff(paste0(prefixPCres,"_UMAP","_by_","cluster","_hires.tiff"))
1595. lowrestiff(paste0(prefixPCres,"_UMAP","_by_","cluster","_lowres.tiff"))
1597. # UMAP of cells in each cluster by treatment without cluster labels
1598. DimPlot(seurat_sct,
1599. reduction = "umap",
1600. label = FALSE,
1601. split.by = "treatment",
1602. pt.size = 0.8,
1603. cols = sorted_palopal,
1604. raster = F
1605. ) +
1606. NoLegend()
1607. hirestiff(paste0(prefixPCres,"_UMAP","_by_","treatment",
1608. "_no_labels","_hires.tiff"))
1609. lowrestiff(paste0(prefixPCres,"_UMAP","_by_","treatment",
1610. "_no_labels","_lowres.tiff"))
1612. # UMAP of cells in each cluster by treatment with cluster labels
1613. DimPlot(seurat_sct,
1614. reduction = "umap",
1615. label = TRUE,
1616. split.by = "treatment",
1617. pt.size = 0.8,
1618. label.size = 5,
1619. cols = sorted_palopal,
1620. raster = F
1621. ) +
1622. NoLegend()
1623. hirestiff(paste0(prefixPCres,"_UMAP","_by_","treatment",
1624. "_with_clusters_labels","_hires.tiff"))
1625. lowrestiff(paste0(prefixPCres,"_UMAP","_by_","treatment",
1626. "_with_clusters_labels","_lowres.tiff"))
1628. # UMAP of cells in each cluster by batch without cluster labels
1629. DimPlot(seurat_sct,
1630. reduction = "umap",
1631. label = FALSE,
1632. split.by = "batch",
1633. pt.size = 0.8,
1634. cols = sorted_palopal,
1635. raster = F
1636. ) +
1637. NoLegend()
1638. hirestiff(paste0(prefixPCres,"_UMAP","_by_","batch",
1639. "_no_labels","_hires.tiff"))
1640. lowrestiff(paste0(prefixPCres,"_UMAP","_by_","batch",
1641. "_no_labels","_lowres.tiff"))
1643. # UMAP of cells in each cluster by batch with cluster labels
1644. DimPlot(seurat_sct,
1645. reduction = "umap",
1646. label = TRUE,
1647. split.by = "batch",
1648. pt.size = 0.8,
1649. label.size = 3,
1650. cols = sorted_palopal,
1651. raster = F
1652. ) +
1653. NoLegend()
1654. hirestiff(paste0(prefixPCres,"_UMAP","_by_","batch",
1655. "_with_clusters_labels","_hires.tiff"))
1656. lowrestiff(paste0(prefixPCres,"_UMAP","_by_","batch",
1657. "_with_clusters_labels","_lowres.tiff"))
1659. # UMAP of cells in each cluster by time_point without cluster labels
1660. DimPlot(seurat_sct,
1661. reduction = "umap",
1662. label = FALSE,
1663. split.by = "time_point",
1664. pt.size = 0.8,
1665. cols = sorted_palopal,
1666. raster = F
1667. ) +
1668. NoLegend()
1669. hirestiff(paste0(prefixPCres,"_UMAP","_by_","time_point",
1670. "_no_labels","_hires.tiff"))
1671. lowrestiff(paste0(prefixPCres,"_UMAP","_by_","time_point",
1672. "_no_labels","_lowres.tiff"))
1674. # UMAP of cells in each cluster by time_point with cluster labels
1675. DimPlot(seurat_sct,
1676. reduction = "umap",
1677. label = TRUE,
1678. split.by = "time_point",
1679. pt.size = 0.8,
1680. label.size = 5,
1681. cols = sorted_palopal,
1682. raster = F
1683. ) +
1684. NoLegend()
1685. hirestiff(paste0(prefixPCres,"_UMAP","_by_","time_point",
1686. "_with_clusters_labels","_hires.tiff"))
1687. lowrestiff(paste0(prefixPCres,"_UMAP","_by_","time_point",
1688. "_with_clusters_labels","_lowres.tiff"))
1690. # UMAP of cells in each cluster by time_point without cluster labels
1691. DimPlot(seurat_sct,
1692. reduction = "umap",
1693. label = FALSE,
1694. split.by = "time_point_treatment",
1695. pt.size = 0.8,
1696. cols = sorted_palopal,
1697. raster = F
1698. ) +
1699. NoLegend()
1700. hirestiff(paste0(prefixPCres,"_UMAP","_by_","time_point_treatment",
1701. "_no_labels","_hires.tiff"))
1702. lowrestiff(paste0(prefixPCres,"_UMAP","_by_","time_point_treatment",
1703. "_no_labels","_lowres.tiff"))
1705. # UMAP of cells in each cluster by time_point with cluster labels
1706. DimPlot(seurat_sct,
1707. reduction = "umap",
1708. label = TRUE,
1709. split.by = "time_point_treatment",
1710. pt.size = 0.8,
1711. label.size = 5,
1712. cols = sorted_palopal,
1713. raster = F
1714. ) +
1715. NoLegend()
1716. hirestiff(paste0(prefixPCres,"_UMAP","_by_","time_point_treatment",
1717. "_with_clusters_labels","_hires.tiff"))
1718. lowrestiff(paste0(prefixPCres,"_UMAP","_by_","time_point_treatment",
1719. "_with_clusters_labels","_lowres.tiff"))
1721. #### INT Extract number of cells per cluster per orig.ident ####
1722. n_cells <- FetchData(seurat_sct, vars = c("ident", "orig.ident")) %>%
1723. dplyr::count(ident, orig.ident) %>%
1724. tidyr::spread(ident, n)
1726. write.csv(n_cells, file = paste0(prefixPCres,"_cells_per_cluster.csv"))
1728. rm(n_cells)
1730. #### INT save RDS containing reduction and cluster idents ####
1731. # saveRDS(seurat_sct, paste0(prefixPCres,"_seurat_after_clustering.rds"))
1733. # seurat_sct <- readRDS(file = paste0(prefixPCres,"_seurat_after_clustering.rds"))
1735. #### INT normalize rna slot ####
1737. # Select the RNA counts slot to be the default assay for visualization purposes
1738. DefaultAssay(seurat_sct) <- "RNA"
1740. # Normalize, find variable features, scale data
1741. seurat_sct <- NormalizeData(seurat_sct)
1742. seurat_sct <- FindVariableFeatures(seurat_sct)
1743. all.genes <- rownames(seurat_sct)
1744. seurat_sct <- ScaleData(seurat_sct, features = all.genes)
1746. # Export normalized counts
1747. norm_counts <- GetAssayData(seurat_sct, slot = "data")
1749. save(norm_counts, file = paste0(prefixPCres,"_normalized_counts_RNA_sparse_matrix.rdata"))
1751. rm(norm_counts)
1753. # save RDS containing RNA normalized data
1754. # saveRDS(seurat_sct, paste0(prefixPCres,"_seurat_after_RNAnorm.rds"))
1755. # seurat_sct <- readRDS(file = paste0(prefixPCres,"_seurat_after_RNAnorm.rds"))
1757. # save RDS containing RNA normalized data AND CELLTYPE LABELS
1758. # saveRDS(seurat_sct, paste0(prefixPCres,"_seurat_after_RNAnorm_CELLTYPES.rds"))
1759. # seurat_sct <- readRDS(file = paste0(prefixPCres,"_seurat_after_RNAnorm_CELLTYPES.rds"))
1761. #### INT Cluster QC ####
1763. # Look at QC metrics for clustering by cell quality
1765. VlnPlot(seurat_sct,
1766. features = "nCount_RNA",
1767. raster = F) +
1768. NoLegend()
1770. hirestiff(paste(prefixPCres,"vln","plot","QC","nCount_RNA","hires.tiff", sep = "_"))
1771. lowrestiff(paste(prefixPCres,"vln","plot","QC","nCount_RNA","lowres.tiff", sep = "_"))
1773. # Low count clusters: nothing obvious
1775. VlnPlot(seurat_sct,
1776. features = "nFeature_RNA",
1777. raster = F) +
1778. NoLegend()
1780. hirestiff(paste(prefixPCres,"vln","plot","QC","nFeature_RNA","hires.tiff", sep = "_"))
1781. lowrestiff(paste(prefixPCres,"vln","plot","QC","nFeature_RNA","lowres.tiff", sep = "_"))
1783. # Low feature clusters: nothing obvious
1785. VlnPlot(seurat_sct,
1786. features = "percent.mt",
1787. raster = F) +
1788. NoLegend()
1790. hirestiff(paste(prefixPCres,"vln","plot","QC","percent.mt","hires.tiff", sep = "_"))
1791. lowrestiff(paste(prefixPCres,"vln","plot","QC","percent.mt","lowres.tiff", sep = "_"))
1793. # High percent mito clusters: various
1795. VlnPlot(seurat_sct,
1796. features = "percent.ribo",
1797. raster = F) +
1798. NoLegend()
1800. hirestiff(paste(prefixPCres,"vln","plot","QC","percent.ribo","hires.tiff", sep = "_"))
1801. lowrestiff(paste(prefixPCres,"vln","plot","QC","percent.ribo","lowres.tiff", sep = "_"))
1803. # High percent ribo clusters: varies
1805. #### PLOTS ####
1806. #### Test Plots ####
1807. # apoptosis genes (bax and casp3 up, other two down)
1808. FeaturePlot(seurat_sct,
1809. reduction = "umap",
1810. features = c("Bax","Casp3","Bnip3","Bcl2"),
1811. order = TRUE,
1812. min.cutoff = 'q10',
1813. label = TRUE,
1814. label.size = 3,
1815. pt.size = 0.5,
1816. raster = F)
1818. # Plot the UMAP
1819. DimPlot(seurat_sct,
1820. reduction = "umap",
1821. label = TRUE,
1822. label.size = 5,
1823. pt.size = 0.8) +
1824. NoLegend()
1826. # Plot the UMAP split by Treatment
1827. DimPlot(seurat_sct,
1828. reduction = "umap",
1829. label = TRUE,
1830. split.by = "treatment",
1831. pt.size = 0.8) +
1832. NoLegend()
1835. Rods <- c("Rho","Gnat1","Nrl","Pde6a","Pde6b","Cngb1")
1836. Rods1 <- c("Rho","Gnat1","Nrl")
1837. Rods2 <- c("Pde6a","Pde6b","Cngb1")
1838. Cones <- c("Arr3","Opn1sw","Gnat2","Pde6h")
1839. Mueller_Glia <- c("Glul","Rlbp1")
1840. Astrocytes <- c("Gfap","Aqp4","Gpr37")
1841. Microglia <- c("Tyrobp","Cx3cr1","Tmem119","Clec7a")
1842. Bipolar <- c("Scgn","Vsx1","Vsx2")
1843. Retinal_Ganglion <- c("Rbpms","Pou4f1","Pou4f2")
1844. Amacrine <- c("Calb1","Gad1","Gad2")
1845. Horizontal <- c("Onecut2")
1846. Pericytes <- c("Pecam1","Acta2")
1847. Reticulocytes <- c("Hbb-bs")
1849. default_FP <- function(obj,features){
1850. FeaturePlot(obj,
1851. reduction = "umap",
1852. features = features,
1853. order = TRUE,
1854. # min.cutoff = 'q10',
1855. label = TRUE,
1856. label.size = 3,
1857. pt.size = 0.5)
1858. }
1860. #### Genes from Dr. Wang ####
1862. FeaturePlot(seurat_sct,
1863. reduction = "umap",
1864. features = c("Ifitm3","Socs3","Slpi","Ch25h"),
1865. order = TRUE,
1866. min.cutoff = 'q10',
1867. label = TRUE,
1868. label.size = 3,
1869. pt.size = 0.5)
1871. FeaturePlot(seurat_sct,
1872. reduction = "umap",
1873. features = c("Rlbp1", # Rlbp1 = Crlbp1
1874. "Anxa1","Anxa2","Nfkbia"),
1875. order = TRUE,
1876. min.cutoff = 'q10',
1877. label = TRUE,
1878. label.size = 3,
1879. pt.size = 0.5)
1881. FeaturePlot(seurat_sct,
1882. reduction = "umap",
1883. features = c("Aif1", # Aif1 = Iba1
1884. "Adgre1", # Adgre1 = F4/80
1885. "Cd68"),
1886. order = TRUE,
1887. min.cutoff = 'q10',
1888. label = TRUE,
1889. label.size = 3,
1890. pt.size = 0.5)
1892. FeaturePlot(seurat_sct,
1893. reduction = "umap",
1894. features = c("Aifm1","Aif1"),
1895. order = TRUE,
1896. min.cutoff = 'q10',
1897. label = TRUE,
1898. label.size = 3,
1899. pt.size = 0.5)
1901. FeaturePlot(seurat_sct,
1902. reduction = "umap",
1903. features = c("Mydgf","Kdelr1","Kdelr2","Ern1"),
1904. order = TRUE,
1905. min.cutoff = 'q10',
1906. label = TRUE,
1907. label.size = 3,
1908. pt.size = 0.5)
1910. FeaturePlot(seurat_sct,
1911. reduction = "umap",
1912. features = c("Mydgf"),
1913. order = TRUE,
1914. min.cutoff = 'q10',
1915. label = TRUE,
1916. label.size = 3,
1917. pt.size = 0.5,
1918. split.by = "time_point_treatment")
1920. FeaturePlot(seurat_sct,
1921. reduction = "umap",
1922. features = c("Kdelr1"),
1923. order = TRUE,
1924. min.cutoff = 'q10',
1925. label = TRUE,
1926. label.size = 3,
1927. pt.size = 0.5,
1928. split.by = "time_point_treatment")
1930. FeaturePlot(seurat_sct,
1931. reduction = "umap",
1932. features = c("Kdelr2"),
1933. order = TRUE,
1934. min.cutoff = 'q10',
1935. label = TRUE,
1936. label.size = 3,
1937. pt.size = 0.5,
1938. split.by = "time_point_treatment")
1940. FeaturePlot(seurat_sct,
1941. reduction = "umap",
1942. features = c("Ern1"),
1943. order = TRUE,
1944. min.cutoff = 'q10',
1945. label = TRUE,
1946. label.size = 3,
1947. pt.size = 0.5,
1948. split.by = "time_point_treatment")
1950. FeaturePlot(seurat_sct,
1951. reduction = "umap",
1952. features = c("Manf"),
1953. order = TRUE,
1954. min.cutoff = 'q10',
1955. label = TRUE,
1956. label.size = 3,
1957. pt.size = 0.5,
1958. split.by = "time_point_treatment")
1960. VlnPlot(seurat_sct,
1961. idents = c("ROD","CONE"),
1962. features = c("Mydgf","Kdelr1","Kdelr2","Ern1","Manf"),
1963. split.by = "time_point_treatment",
1964. stack = T,
1965. flip = T)
1967. FeaturePlot(seurat_sct,
1968. reduction = "umap",
1969. features = c("Creb1",
1970. "Crem"),
1971. order = TRUE,
1972. min.cutoff = 'q10',
1973. label = TRUE,
1974. label.size = 3,
1975. pt.size = 0.5)
1977. FeaturePlot(seurat_sct,
1978. reduction = "umap",
1979. features = c("Tmem119","P2ry13","P2ry12","Siglech"),
1980. order = TRUE,
1981. min.cutoff = 'q10',
1982. label = TRUE,
1983. label.size = 3,
1984. pt.size = 0.5)
1986. #### retina genes from Saba SCT ####
1988. # Rod Photoreceptors
1989. # Cluster 7
1990. VlnPlot(seurat_sct,
1991. features = c("Rho","Gnat1","Nrl",
1992. "Pde6a","Pde6b","Cngb1"),
1993. stack = TRUE,
1994. flip = TRUE) +
1995. NoLegend() +
1996. ggtitle("Rod Photoreceptors")
1998. hirestiff(paste(prefixPCres,"VLN","Rod","Photoreceptors","hires.tiff", sep = "_"))
1999. lowrestiff(paste(prefixPCres,"VLN","Rod","Photoreceptors","lowres.tiff", sep = "_"))
2001. # Cone Photoreceptors
2002. # Cluster 8
2003. VlnPlot(seurat_sct,
2004. features = c("Arr3","Opn1sw","Gnat2","Pde6h"),
2005. stack = TRUE,
2006. flip = TRUE) +
2007. NoLegend() +
2008. ggtitle("Cone Photoreceptors")
2010. hirestiff(paste(prefixPCres,"VLN","Cone","Photoreceptors","hires.tiff", sep = "_"))
2011. lowrestiff(paste(prefixPCres,"VLN","Cone","Photoreceptors","lowres.tiff", sep = "_"))
2013. # Muller Glia
2014. # Clusters 4, 10, 12
2015. VlnPlot(seurat_sct,
2016. features = c("Glul","Rlbp1","Nes","Vim"),
2017. stack = TRUE,
2018. flip = TRUE) +
2019. NoLegend() +
2020. ggtitle("Muller Glia")
2022. hirestiff(paste(prefixPCres,"VLN","Muller","Glia","hires.tiff", sep = "_"))
2023. lowrestiff(paste(prefixPCres,"VLN","Muller","Glia","lowres.tiff", sep = "_"))
2025. # Retinal Astrocytes
2026. # Clusters 4, 10, 12
2027. VlnPlot(seurat_sct,
2028. features = c("Gfap","Aqp4","Gpr37"),
2029. stack = TRUE,
2030. flip = TRUE) +
2031. NoLegend() +
2032. ggtitle("Retinal Astrocytes")
2034. hirestiff(paste(prefixPCres,"VLN","Retinal","Astrocytes","hires.tiff", sep = "_"))
2035. lowrestiff(paste(prefixPCres,"VLN","Retinal","Astrocytes","lowres.tiff", sep = "_"))
2037. # Microglia
2038. # Clusters 3, 13, 16
2039. VlnPlot(seurat_sct,
2040. features = c("Tyrobp","Cx3cr1","Tmem119"),
2041. stack = TRUE,
2042. flip = TRUE) +
2043. NoLegend() +
2044. ggtitle("Microglia")
2046. hirestiff(paste(prefixPCres,"VLN","Microglia","hires.tiff", sep = "_"))
2047. lowrestiff(paste(prefixPCres,"VLN","Microglia","lowres.tiff", sep = "_"))
2049. # Bipolar Cells
2050. # Clusters 0, 1, (4), 6, 9, (10), ?11 ,(12,14)
2051. VlnPlot(seurat_sct,
2052. features = c("Scgn","Vsx1","Vsx2"),
2053. stack = TRUE,
2054. flip = TRUE) +
2055. NoLegend() +
2056. ggtitle("Bipolar Cells")
2058. hirestiff(paste(prefixPCres,"VLN","Bipolar","Cells","hires.tiff", sep = "_"))
2059. lowrestiff(paste(prefixPCres,"VLN","Bipolar","Cells","lowres.tiff", sep = "_"))
2061. # Retinal Ganglion Cells
2062. # Clusters 2, 5, 13, 14, 16
2063. VlnPlot(seurat_sct,
2064. features = c("Rbpms","Pou4f1","Pou4f2"),
2065. stack = TRUE,
2066. flip = TRUE) +
2067. NoLegend() +
2068. ggtitle("Retinal Ganglion Cells")
2070. hirestiff(paste(prefixPCres,"VLN","Retinal","Ganglion","Cells","hires.tiff", sep = "_"))
2071. lowrestiff(paste(prefixPCres,"VLN","Retinal","Ganglion","Cells","lowres.tiff", sep = "_"))
2073. # Amacrine Cells
2074. # Part of Cluster 6
2075. VlnPlot(seurat_sct,
2076. features = c("Calb1","Gad1","Gad2"),
2077. stack = TRUE,
2078. flip = TRUE) +
2079. NoLegend() +
2080. ggtitle("Amacrine Cells")
2082. hirestiff(paste(prefixPCres,"VLN","Amacrine","Cells","hires.tiff", sep = "_"))
2083. lowrestiff(paste(prefixPCres,"VLN","Amacrine","Cells","lowres.tiff", sep = "_"))
2085. # Horizontal Cells
2086. # Cluster
2087. VlnPlot(seurat_sct,
2088. features = c("Onecut2"),
2089. # stack = TRUE,
2090. flip = TRUE,
2091. pt.size = 0) +
2092. NoLegend() +
2093. ggtitle("Horizontal Cells")
2095. hirestiff(paste(prefixPCres,"VLN","Horizontal","Cells","hires.tiff", sep = "_"))
2096. lowrestiff(paste(prefixPCres,"VLN","Horizontal","Cells","lowres.tiff", sep = "_"))
2098. # Pericytes
2099. # Clusters 2,16
2100. VlnPlot(seurat_sct,
2101. features = c(#"Pecam1",
2102. "Acta2"),
2103. # stack = TRUE,
2104. flip = TRUE,
2105. pt.size = 0) +
2106. NoLegend() #+ ggtitle("Pericytes")
2108. hirestiff(paste(prefixPCres,"VLN","Pericytes","hires.tiff", sep = "_"))
2109. lowrestiff(paste(prefixPCres,"VLN","Pericytes","lowres.tiff", sep = "_"))
2111. # Reticulocytes
2112. # Cluster 15
2113. VlnPlot(seurat_sct,
2114. features = c(#"Pecam1",
2115. "Hbb-bs"),
2116. # stack = TRUE,
2117. flip = TRUE,
2118. pt.size = 0) +
2119. NoLegend() #+ ggtitle("Pericytes")
2121. hirestiff(paste(prefixPCres,"VLN","Pericytes","hires.tiff", sep = "_"))
2122. lowrestiff(paste(prefixPCres,"VLN","Pericytes","lowres.tiff", sep = "_"))
2124. #### retina genes from Saba Integrated ####
2126. VlnPlot(seurat_sct,
2127. features = c("Nfe2l2"),
2128. # stack = TRUE,
2129. flip = TRUE) +
2130. NoLegend()
2132. default_FP(seurat_sct,c("Nfe2l2"))
2134. VlnPlot(seurat_sct,
2135. features = c(RPEpi,"Mertk"),
2136. stack = TRUE,
2137. flip = TRUE) +
2138. NoLegend()
2140. # Rod Photoreceptors
2141. # Cluster 8
2142. VlnPlot(seurat_sct,
2143. features = Rods,
2144. stack = TRUE,
2145. flip = TRUE) +
2146. NoLegend() +
2147. ggtitle("Rod Photoreceptors")
2149. hirestiff(paste(prefixPCres,"VLN","Rod","Photoreceptors","hires.tiff", sep = "_"))
2150. lowrestiff(paste(prefixPCres,"VLN","Rod","Photoreceptors","lowres.tiff", sep = "_"))
2152. default_FP(seurat_sct,Rods1)
2154. hirestiff(paste(prefixPCres,"UMAP","Rod","genes1","hires.tiff", sep = "_"))
2155. lowrestiff(paste(prefixPCres,"UMAP","Rod","genes1","lowres.tiff", sep = "_"))
2157. default_FP(seurat_sct,Rods2)
2159. hirestiff(paste(prefixPCres,"UMAP","Rod","genes2","hires.tiff", sep = "_"))
2160. lowrestiff(paste(prefixPCres,"UMAP","Rod","genes2","lowres.tiff", sep = "_"))
2163. default_FP(seurat_sct,"Grm6")
2165. default_FP(seurat_sct,c("Pou4f1","Pou4f2","Pou4f3"))
2167. VlnPlot(seurat_sct,
2168. features = c("Onecut1","Onecut2"),
2169. stack = TRUE,
2170. flip = TRUE) +
2171. NoLegend()
2173. # Cone Photoreceptors
2174. # Cluster 9
2175. VlnPlot(seurat_sct,
2176. features = Cones,
2177. stack = TRUE,
2178. flip = TRUE) +
2179. NoLegend() +
2180. ggtitle("Cone Photoreceptors")
2182. hirestiff(paste(prefixPCres,"VLN","Cone","Photoreceptors","hires.tiff", sep = "_"))
2183. lowrestiff(paste(prefixPCres,"VLN","Cone","Photoreceptors","lowres.tiff", sep = "_"))
2185. default_FP(seurat_sct,Cones)
2187. hirestiff(paste(prefixPCres,"UMAP","Cone","genes","hires.tiff", sep = "_"))
2188. lowrestiff(paste(prefixPCres,"UMAP","Cone","genes","lowres.tiff", sep = "_"))
2190. FeaturePlot(seurat_sct,
2191. reduction = "umap",
2192. features = Cones,
2193. order = TRUE,
2194. min.cutoff = 'q10',
2195. label = TRUE,
2196. label.size = 3,
2197. pt.size = 0.5)
2199. # Muller Glia
2200. # Clusters 10,11,14,15,16,20,24,27,29
2201. VlnPlot(seurat_sct,
2202. features = c(Mueller_Glia,"Gfap","Aqp4","S100b"),
2203. stack = TRUE,
2204. flip = TRUE) +
2205. NoLegend() +
2206. ggtitle("Muller Glia")
2208. hirestiff(paste(prefixPCres,"VLN","Muller","Glia","hires.tiff", sep = "_"))
2209. lowrestiff(paste(prefixPCres,"VLN","Muller","Glia","lowres.tiff", sep = "_"))
2211. default_FP(seurat_sct,Mueller_Glia)
2213. hirestiff(paste(prefixPCres,"UMAP","Mueller_Glia","genes","hires.tiff", sep = "_"))
2214. lowrestiff(paste(prefixPCres,"UMAP","Mueller_Glia","genes","lowres.tiff", sep = "_"))
2216. Macroglia <- c("Apoe","Slc1a3","Trpm3","Glul","Clu" )
2219. FeaturePlot(seurat_sct,
2220. reduction = "umap",
2221. features = Macroglia,
2222. order = TRUE,
2223. min.cutoff = 'q10',
2224. label = TRUE,
2225. label.size = 3,
2226. pt.size = 0.5)
2228. FeaturePlot(seurat_sct,
2229. reduction = "umap",
2230. features = Mueller_Glia,
2231. order = TRUE,
2232. min.cutoff = 'q10',
2233. label = TRUE,
2234. label.size = 3,
2235. pt.size = 0.5)
2237. FeaturePlot(seurat_sct,
2238. reduction = "umap",
2239. features = c("Glul"),
2240. order = TRUE,
2241. min.cutoff = 'q10',
2242. label = TRUE,
2243. label.size = 3,
2244. pt.size = 0.5,
2245. split.by = "time_point")
2247. FeaturePlot(TRseurat,
2248. reduction = "umap",
2249. features = c("Glul"),
2250. order = TRUE,
2251. min.cutoff = 'q10',
2252. label = TRUE,
2253. label.size = 3,
2254. pt.size = 0.5,
2255. split.by = "time_point")
2257. FeaturePlot(UTseurat,
2258. reduction = "umap",
2259. features = c("Glul"),
2260. order = TRUE,
2261. min.cutoff = 'q10',
2262. label = TRUE,
2263. label.size = 3,
2264. pt.size = 0.5,
2265. split.by = "time_point")
2267. # Retinal Astrocytes
2268. # Clusters 10,11,14,15,16,20,24,27,29
2269. VlnPlot(seurat_sct,
2270. features = Astrocytes,
2271. stack = TRUE,
2272. flip = TRUE) +
2273. NoLegend() +
2274. ggtitle("Retinal Astrocytes")
2276. hirestiff(paste(prefixPCres,"VLN","Retinal","Astrocytes","hires.tiff", sep = "_"))
2277. lowrestiff(paste(prefixPCres,"VLN","Retinal","Astrocytes","lowres.tiff", sep = "_"))
2279. default_FP(seurat_sct,Astrocytes)
2281. hirestiff(paste(prefixPCres,"UMAP","Astrocytes","genes","hires.tiff", sep = "_"))
2282. lowrestiff(paste(prefixPCres,"UMAP","Astrocytes","genes","lowres.tiff", sep = "_"))
2284. # Microglia
2285. # Clusters 5,7,13,21,22,26,28,29
2286. VlnPlot(seurat_sct,
2287. features = Microglia,
2288. stack = TRUE,
2289. flip = TRUE) +
2290. NoLegend() +
2291. ggtitle("Microglia")
2293. hirestiff(paste(prefixPCres,"VLN","Microglia","hires.tiff", sep = "_"))
2294. lowrestiff(paste(prefixPCres,"VLN","Microglia","lowres.tiff", sep = "_"))
2296. default_FP(seurat_sct,Microglia)
2298. hirestiff(paste(prefixPCres,"UMAP","Microglia","genes","hires.tiff", sep = "_"))
2299. lowrestiff(paste(prefixPCres,"UMAP","Microglia","genes","lowres.tiff", sep = "_"))
2301. # Bipolar Cells
2302. # Clusters (0),1,2,3
2303. VlnPlot(seurat_sct,
2304. features = c(Bipolar,"Snap25","Neto1"),
2305. stack = TRUE,
2306. flip = TRUE) +
2307. NoLegend() +
2308. ggtitle("Bipolar Cells")
2310. hirestiff(paste(prefixPCres,"VLN","Bipolar","Cells","hires.tiff", sep = "_"))
2311. lowrestiff(paste(prefixPCres,"VLN","Bipolar","Cells","lowres.tiff", sep = "_"))
2313. default_FP(seurat_sct,Bipolar)
2315. hirestiff(paste(prefixPCres,"UMAP","Bipolar","genes","hires.tiff", sep = "_"))
2316. lowrestiff(paste(prefixPCres,"UMAP","Bipolar","genes","lowres.tiff", sep = "_"))
2318. # Retinal Ganglion Cells
2319. # Clusters 4,6,12,18,19,21,24,25,26,28
2320. VlnPlot(seurat_sct,
2321. features = c(Retinal_Ganglion,
2322. "Slc17a6","Pou4f3","Cacna1c"),
2323. stack = TRUE,
2324. flip = TRUE,
2325. raster = F) +
2326. NoLegend() +
2327. ggtitle("Retinal Ganglion Cells")
2329. hirestiff(paste(prefixPCres,"VLN","Retinal","Ganglion","Cells","hires.tiff", sep = "_"))
2330. lowrestiff(paste(prefixPCres,"VLN","Retinal","Ganglion","Cells","lowres.tiff", sep = "_"))
2332. default_FP(seurat_sct,Retinal_Ganglion)
2334. hirestiff(paste(prefixPCres,"UMAP","Retinal_Ganglion","genes","hires.tiff", sep = "_"))
2335. lowrestiff(paste(prefixPCres,"UMAP","Retinal_Ganglion","genes","lowres.tiff", sep = "_"))
2337. FeaturePlot(seurat_sct,
2338. reduction = "umap",
2339. features = c(Retinal_Ganglion,
2340. "Slc17a6","Pou4f3"),
2341. order = TRUE,
2342. min.cutoff = 'q10',
2343. label = TRUE,
2344. label.size = 3,
2345. pt.size = 0.5,
2346. raster = F)
2348. # Amacrine Cells
2349. # Cluster 17
2350. VlnPlot(seurat_sct,
2351. features = Amacrine,
2352. stack = TRUE,
2353. flip = TRUE) +
2354. NoLegend() +
2355. ggtitle("Amacrine Cells")
2357. hirestiff(paste(prefixPCres,"VLN","Amacrine","Cells","hires.tiff", sep = "_"))
2358. lowrestiff(paste(prefixPCres,"VLN","Amacrine","Cells","lowres.tiff", sep = "_"))
2360. default_FP(seurat_sct,Amacrine)
2362. hirestiff(paste(prefixPCres,"UMAP","Amacrine","genes","hires.tiff", sep = "_"))
2363. lowrestiff(paste(prefixPCres,"UMAP","Amacrine","genes","lowres.tiff", sep = "_"))
2365. # Horizontal Cells
2366. # Clusters ?
2367. VlnPlot(seurat_sct,
2368. features = c(Horizontal,"Onecut1","Gfap"),
2369. stack = TRUE,
2370. flip = TRUE,
2371. pt.size = 0) +
2372. NoLegend() +
2373. ggtitle("Horizontal Cells")
2375. hirestiff(paste(prefixPCres,"VLN","Horizontal","Cells","hires.tiff", sep = "_"))
2376. lowrestiff(paste(prefixPCres,"VLN","Horizontal","Cells","lowres.tiff", sep = "_"))
2378. default_FP(seurat_sct,Horizontal)
2380. hirestiff(paste(prefixPCres,"UMAP","Horizontal","genes","hires.tiff", sep = "_"))
2381. lowrestiff(paste(prefixPCres,"UMAP","Horizontal","genes","lowres.tiff", sep = "_"))
2383. # Pericytes
2384. # Cluster 4,12,18,24,25,26
2385. VlnPlot(seurat_sct,
2386. features = c(Pericytes,"Cacna1c"),
2387. stack = TRUE,
2388. flip = TRUE,
2389. pt.size = 0) +
2390. NoLegend() +
2391. ggtitle("Pericytes")
2393. hirestiff(paste(prefixPCres,"VLN","Pericytes","hires.tiff", sep = "_"))
2394. lowrestiff(paste(prefixPCres,"VLN","Pericytes","lowres.tiff", sep = "_"))
2396. default_FP(seurat_sct,Pericytes)
2398. hirestiff(paste(prefixPCres,"UMAP","Pericytes","genes","hires.tiff", sep = "_"))
2399. lowrestiff(paste(prefixPCres,"UMAP","Pericytes","genes","lowres.tiff", sep = "_"))
2401. # Reticulocytes 23
2402. VlnPlot(seurat_sct,
2403. features = c(Reticulocytes,"Gfap"),
2404. stack = TRUE,
2405. flip = TRUE,
2406. pt.size = 0) +
2407. NoLegend() +
2408. ggtitle("Reticulocytes")
2410. hirestiff(paste(prefixPCres,"VLN","Reticulocytes","hires.tiff", sep = "_"))
2411. lowrestiff(paste(prefixPCres,"VLN","Reticulocytes","lowres.tiff", sep = "_"))
2413. default_FP(seurat_sct,Reticulocytes)
2415. hirestiff(paste(prefixPCres,"UMAP","Reticulocytes","genes","hires.tiff", sep = "_"))
2416. lowrestiff(paste(prefixPCres,"UMAP","Reticulocytes","genes","lowres.tiff", sep = "_"))
2418. #### retinal genes ####
2420. Rods <- c("Rho","Gnat1","Nrl","Pde6a","Pde6b","Cngb1")
2421. Rods1 <- c("Rho","Gnat1","Nrl")
2422. Rods2 <- c("Pde6a","Pde6b","Cngb1")
2423. Cones <- c("Arr3","Opn1sw","Gnat2","Pde6h")
2424. Mueller_Glia <- c("Glul","Rlbp1")
2425. Astrocytes <- c("Gfap","Aqp4","Gpr37")
2426. Microglia <- c("Tyrobp","Cx3cr1","Tmem119")
2427. Bipolar <- c("Scgn","Vsx1","Vsx2")
2428. Retinal_Ganglion <- c("Rbpms","Pou4f1","Pou4f2")
2429. Amacrine <- c("Calb1","Gad1","Gad2")
2430. Horizontal <- c("Onecut2")
2431. Pericytes <- c("Pecam1","Acta2")
2432. Reticulocytes <- c("Hbb-bs")
2433. RPEpi <- c("Mitf","Tjp1","Rpe65","Rlbp1","Best1")
2435. retina_genes <- c(Rods,Cones,Mueller_Glia,Astrocytes,
2436. Microglia,Bipolar,Retinal_Ganglion,
2437. Amacrine,Horizontal,Reticulocytes)
2439. FeaturePlot(seurat_sct,
2440. reduction = "umap",
2441. features = c(Horizontal,"Lhx1","Isl1","Prox1"),
2442. order = TRUE,
2443. min.cutoff = 'q10',
2444. label = TRUE,
2445. label.size = 3,
2446. pt.size = 0.5)
2448. FeaturePlot(seurat_sct,
2449. reduction = "umap",
2450. features = c("Ntrk2"),
2451. order = TRUE,
2452. min.cutoff = 'q10',
2453. label = TRUE,
2454. label.size = 3,
2455. pt.size = 0.5)
2457. #### retinal genes from Chen et al 2022 ####
2459. Rods <- c("Rho","Nr2e3","Nrl","Pdc","Rp1")
2460. Cones <- c("Opn1sw","Opn1mw","Pde6h","Arr3","Gnat2")
2461. Mueller_Glia <- c("Zfp36l1","Dbi","Apoe","Slc1a3","Sparc")
2462. Bipolar <- c("Pcp2","Isl1","Grm6","Trnp1")
2463. Amacrine <- c("C1ql1","Snhg11","Tfap2b","Pcsk1n")
2464. Horizontal <- c("Slc4a3","Calb1","Septin4","Tpm3")
2465. Microglia <- c("Ctsd","Ccl4","C1qb","C1qc")
2466. Vascular <- c("Trpm1","Igfbp7")
2468. # Astrocytes <- c("Gfap","Aqp4","Gpr37")
2469. # Retinal_Ganglion <- c("Rbpms","Pou4f1","Pou4f2")
2470. # Pericytes <- c("Pecam1","Acta2")
2471. # Reticulocytes <- c("Hbb-bs")
2473. retina_genes <- c(Rods,Cones,Mueller_Glia,Bipolar,
2474. Amacrine,Horizontal,Microglia,Vascular)
2476. Idents(seurat_sct) <- ident_res
2478. Idents(seurat_sct) <- "Cell.Type"
2480. Idents(seurat_sct) <- "Cluster.Type"
2482. VlnPlot(seurat_sct,
2483. features = c("Olig2","Apc","Cnp",
2484. "Nes","Fabp7","Pax6"),
2485. stack = T,
2486. flip = T) + NoLegend()
2488. VlnPlot(seurat_sct,
2489. features = c("Casp3","Casp6","Bax",
2490. "Bcl2","Fas","Faslg","Tnfrsf1a"),
2491. stack = T,
2492. flip = T) + NoLegend()
2494. VlnPlot(seurat_sct,
2495. features = Rods,
2496. stack = T,
2497. flip = T) + NoLegend()
2499. VlnPlot(seurat_sct,
2500. features = Cones,
2501. stack = T,
2502. flip = T) + NoLegend()
2504. VlnPlot(seurat_sct,
2505. features = Mueller_Glia,
2506. stack = T,
2507. flip = T) + NoLegend()
2509. VlnPlot(seurat_sct,
2510. features = Retinal_Ganglion,
2511. stack = T,
2512. flip = T) + NoLegend()
2514. VlnPlot(seurat_sct,
2515. features = c("Glul","Rlbp1","Gfap","S100b"),
2516. stack = T,
2517. flip = T) + NoLegend()
2519. VlnPlot(seurat_sct,
2520. features = Bipolar,
2521. stack = T,
2522. flip = T) + NoLegend()
2524. VlnPlot(seurat_sct,
2525. features = c("Isl1","Grm6","Scgn","Vsx2","Prdm8",
2526. "Prkca","Pcp2","Bhlhe23","Cabp5","Car8",
2527. "Trpm1","Slc5a8","Cacna2d3","Adamts5"),
2528. stack = T,
2529. flip = T) + NoLegend()
2531. VlnPlot(seurat_sct,
2532. features = c("Gsg1","Tmem215","Trnp1"),
2533. stack = T,
2534. flip = T) + NoLegend()
2536. VlnPlot(seurat_sct,
2537. features = Amacrine,
2538. stack = T,
2539. flip = T) + NoLegend()
2541. VlnPlot(seurat_sct,
2542. features = c("Snhg11","Tfap2b","Gad1","Gad2"),
2543. stack = T,
2544. flip = T) + NoLegend()
2546. VlnPlot(seurat_sct,
2547. features = Horizontal,
2548. stack = T,
2549. flip = T) + NoLegend()
2551. VlnPlot(seurat_sct,
2552. features = Microglia,
2553. stack = T,
2554. flip = T) + NoLegend()
2556. VlnPlot(seurat_sct,
2557. features = c("Ccl4","C1qb","C1qc","Tyrobp","Cx3cr1","Tmem119"),
2558. stack = T,
2559. flip = T) + NoLegend()
2561. VlnPlot(seurat_sct,
2562. features = Vascular,
2563. stack = T,
2564. flip = T) + NoLegend()
2566. VlnPlot(seurat_sct,
2567. features = Retinal_Ganglion,
2568. stack = T,
2569. flip = T) + NoLegend()
2571. VlnPlot(seurat_sct,
2572. features = Reticulocytes,
2573. # stack = T,
2574. flip = T) + NoLegend()
2576. VlnPlot(seurat_sct,
2577. features = c("Ruvbl1","Tfrc","Ucp3","Folr1",
2578. "Ube2o","Cd36","Itga4","Slc11a2"),
2579. stack = T,
2580. flip = T) + NoLegend() # reticulocyte markers
2582. VlnPlot(seurat_sct,
2583. features = c("Tfrc","Cga","Wrn","Kit","Piezo1",
2584. "Kmt5a","Tal1","Gata1","Lmo2","Eng",
2585. "Cfp","Itgb1","Use1","Grsf1","Lrp8",
2586. "Lmnb1","Gypa","Tfr2","Igbp1","Epor",
2587. "Ets1","Bpgm","Cd44","Cdh1"),
2588. stack = T,
2589. flip = T) + NoLegend() # Erythroblast markers
2591. DotPlot(seurat_sct,
2592. features = c("Tfrc","Cga","Wrn","Kit","Piezo1",
2593. "Kmt5a","Tal1","Gata1","Lmo2","Eng",
2594. "Cfp","Itgb1","Use1","Grsf1","Lrp8",
2595. "Lmnb1","Gypa","Tfr2","Igbp1","Epor",
2596. "Ets1","Bpgm","Cd44","Cdh1"),
2597. cluster.idents = T)
2601. # Erythroid-like and erythroid precursor cells
2602. DotPlot(seurat_sct,
2603. features = c("Ahsp","Hba-a1","Hba-a2","Hbb-bs"),
2604. cluster.idents = T)
2606. # Muller
2607. DotPlot(seurat_sct,
2608. features = c("S100a16","Dbi","Gnai2","Crb1",
2609. "Rdh10","Abca8a","Crot","Dapl1",
2610. "Itm2b"),
2611. cluster.idents = T)
2613. default_FP(seurat_sct,Rods)
2614. default_FP(seurat_sct,Cones)
2615. default_FP(seurat_sct,Mueller_Glia)
2616. default_FP(seurat_sct,Bipolar)
2617. default_FP(seurat_sct,Amacrine)
2618. default_FP(seurat_sct,Horizontal)
2619. default_FP(seurat_sct,Microglia)
2620. default_FP(seurat_sct,Vascular)
2623. default_FP(seurat_sct,c("Trpm1","Prkca","Pcp2"))
2625. default_FP(seurat_sct,c("Tp53"))
2627. VlnPlot(seurat_sct,
2628. features = c("Tp53")) + NoLegend()
2630. FeaturePlot(seurat_sct,
2631. reduction = "umap",
2632. features = c("Trpm1","Prkca","Pcp2","Bhlhe23"),
2633. order = TRUE,
2634. min.cutoff = 'q10',
2635. label = TRUE,
2636. label.size = 3,
2637. pt.size = 0.5)
2639. FeaturePlot(seurat_sct,
2640. reduction = "umap",
2641. features = c("Prdm8","Prkca","Pcp2","Bhlhe23"),
2642. order = TRUE,
2643. min.cutoff = 'q10',
2644. label = TRUE,
2645. label.size = 3,
2646. pt.size = 0.5)
2648. #### dotplots ####
2650. DotPlot(seurat_sct,
2651. features = c("Rho", # 7 rods
2652. "Arr3", # 8 cones
2653. "Rlbp1", # 4,10,12 muller glia
2654. "Scgn", # 6 BC
2655. "Sebox", # 0,1 RBCs
2656. "Slc6a9", # 8,7,2,6 Gly-AC
2657. "Gad1", # 2,6 GABA-AC
2658. "Thy1", # 2,6 GABA-AC / GLY -AC
2659. "Onecut2",# 6 HC
2660. "Cx3cr1", # 16,3 MG
2661. # "Pecam1", # NONE Pericytes
2662. "Acta2", # 2,16 Pericytes
2663. "Hbb-bs"),# 15 Reticulocytes
2664. cluster.idents = TRUE)
2666. DotPlot(seurat_sct,
2667. features = c("Rho","Gnat1","Nrl","Pde6a","Pde6b","Cngb1", # rods
2668. "Arr3","Opn1sw","Gnat2","Pde6h", # cones
2669. "Glul","Rlbp1", # mueller glia
2670. "Gfap","Aqp4","Gpr37", # astrocytes
2671. "Tyrobp","Cx3cr1","Tmem119", # Microglia
2672. "Scgn","Vsx1","Vsx2", # Bipolar
2673. "Rbpms","Pou4f1","Pou4f2", # Retinal Ganglion
2674. "Calb1","Gad1","Gad2", # Amacrine
2675. "Onecut2", # Horizontal
2676. "Pecam1","Acta2", # Pericytes
2677. "Hbb-bs"), # Reticulocytes
2678. cluster.idents = TRUE) +
2679. theme(axis.text.x = element_text(angle = 45, vjust = 0.5))
2681. Idents(p60seurat) <- ident_res
2683. DotPlot(p60seurat,
2684. features = c("Rho","Gnat1","Nrl","Pde6a","Pde6b","Cngb1", # rods
2685. "Arr3","Opn1sw","Gnat2","Pde6h", # cones
2686. "Glul","Rlbp1", # mueller glia
2687. "Gfap","Aqp4","Gpr37", # astrocytes
2688. "Tyrobp","Cx3cr1","Tmem119", # Microglia
2689. "Scgn","Vsx1","Vsx2", # Bipolar
2690. "Rbpms","Pou4f1","Pou4f2", # Retinal Ganglion
2691. "Calb1","Gad1","Gad2", # Amacrine
2692. "Onecut2", # Horizontal
2693. "Pecam1","Acta2", # Pericytes
2694. "Hbb-bs"), # Reticulocytes
2695. cluster.idents = TRUE) +
2696. theme(axis.text.x = element_text(angle = 45, vjust = 0.5))
2698. Idents(p90seurat) <- ident_res
2700. DotPlot(p90seurat,
2701. features = c("Rho","Gnat1","Nrl","Pde6a","Pde6b","Cngb1", # rods
2702. "Arr3","Opn1sw","Gnat2","Pde6h", # cones
2703. "Glul","Rlbp1", # mueller glia
2704. "Gfap","Aqp4","Gpr37", # astrocytes
2705. "Tyrobp","Cx3cr1","Tmem119", # Microglia
2706. "Scgn","Vsx1","Vsx2", # Bipolar
2707. "Rbpms","Pou4f1","Pou4f2", # Retinal Ganglion
2708. "Calb1","Gad1","Gad2", # Amacrine
2709. "Onecut2", # Horizontal
2710. "Pecam1","Acta2", # Pericytes
2711. "Hbb-bs"), # Reticulocytes
2712. cluster.idents = TRUE) +
2713. theme(axis.text.x = element_text(angle = 45, vjust = 0.5))
2715. DotPlot(seurat_sct,
2716. features = c(Rods,Cones,Mueller_Glia,Astrocytes,
2717. Microglia,Bipolar,Retinal_Ganglion,
2718. Amacrine,Horizontal,Pericytes,Reticulocytes),
2719. cluster.idents = TRUE) +
2720. theme(axis.text.x = element_text(angle = 45, vjust = 0.5))
2722. VlnPlot(seurat_sct,
2723. features = Retinal_Ganglion,
2724. stack = TRUE,
2725. flip = TRUE,
2726. sort = TRUE) + NoLegend()
2728. # rods 7
2729. VlnPlot(seurat_sct,
2730. features = Rods,
2731. stack = TRUE,
2732. flip = TRUE,
2733. sort = FALSE) + NoLegend()
2735. FeaturePlot()
2737. # cones 8
2738. VlnPlot(seurat_sct,
2739. features = Cones,
2740. stack = TRUE,
2741. flip = TRUE,
2742. sort = TRUE) + NoLegend()
2744. # Mueller_Glia
2745. VlnPlot(seurat_sct,
2746. features = Mueller_Glia,
2747. stack = TRUE,
2748. flip = TRUE,
2749. sort = TRUE) + NoLegend()
2754. #### Identify cells with non-zero expression for each gene ####
2755. cells_with_expression <- WhichCells(seurat_sct, expression = )
2756. all.genes <- rownames(seurat_sct)
2758. # Create a list of all genes with expression greater than zero
2759. genes_with_expression <- names(which(Sums(cells_with_expression) > 0))
2761. which(rowSums(seurat_sct) > 0)
2763. rowSums(seurat_sct) > 0
2764. namesgenes <- names(rowSums(seurat_sct) > 0)
2765. FeaturePlot(seurat_sct,
2766. reduction = "umap",
2767. features = c("Opn1sw"),
2768. order = TRUE,
2769. min.cutoff = 'q10',
2770. label = TRUE,
2771. label.size = 3,
2772. pt.size = 0.5)
2774. FeaturePlot(seurat_sct,
2775. reduction = "umap",
2776. features = c("Arr3"),
2777. order = TRUE,
2778. min.cutoff = 'q10',
2779. label = TRUE,
2780. label.size = 3,
2781. pt.size = 0.5)
2783. FeaturePlot(seurat_sct,
2784. reduction = "umap",
2785. features = c("Pax6","Tfap2b","Gad1","Slc6a9"),
2786. order = TRUE,
2787. min.cutoff = 'q10',
2788. label = TRUE,
2789. label.size = 3,
2790. pt.size = 0.5)
2792. FeaturePlot(seurat_sct,
2793. reduction = "umap",
2794. features = c("Tfap2b"),
2795. order = TRUE,
2796. min.cutoff = 'q10',
2797. label = FALSE,
2798. label.size = 3,
2799. pt.size = 0.5,
2800. split.by = "treatment")
2802. FeaturePlot(seurat_sct,
2803. reduction = "umap",
2804. features = c("Tmem119","Adgre1","Tyrobp"),
2805. order = TRUE,
2806. min.cutoff = 'q10',
2807. label = TRUE,
2808. label.size = 3,
2809. pt.size = 0.5)
2811. FeaturePlot(seurat_sct,
2812. reduction = "umap",
2813. features = c("Glul"),
2814. order = TRUE,
2815. min.cutoff = 'q10',
2816. label = TRUE,
2817. label.size = 3,
2818. pt.size = 0.5)
2820. FeaturePlot(seurat_sct,
2821. reduction = "umap",
2822. features = c("Gfap"),
2823. order = TRUE,
2824. min.cutoff = 'q10',
2825. label = TRUE,
2826. label.size = 3,
2827. pt.size = 0.5)
2829. FeaturePlot(seurat_sct,
2830. reduction = "umap",
2831. features = c("Gsg1","Vsx1"),
2832. order = TRUE,
2833. min.cutoff = 'q10',
2834. label = TRUE,
2835. label.size = 3,
2836. pt.size = 0.5)
2838. FeaturePlot(seurat_sct,
2839. reduction = "umap",
2840. features = c("Onecut1"),
2841. order = TRUE,
2842. min.cutoff = 'q10',
2843. label = TRUE,
2844. label.size = 3,
2845. pt.size = 0.5)
2847. FeaturePlot(seurat_sct,
2848. reduction = "umap",
2849. features = c("Grm6"),
2850. order = TRUE,
2851. min.cutoff = 'q10',
2852. label = TRUE,
2853. label.size = 3,
2854. pt.size = 0.5)
2856. FeaturePlot(seurat_sct,
2857. reduction = "umap",
2858. features = c("Rbpms"),
2859. order = TRUE,
2860. min.cutoff = 'q10',
2861. label = TRUE,
2862. label.size = 3,
2863. pt.size = 0.5)
2865. VlnPlot(seurat_sct,
2866. features = "Grm6",
2867. flip = TRUE,
2868. pt.size = 0)
2872. #### Wang et al 2022 ####
2874. DimPlot(seurat_sct,
2875. reduction = "umap",
2876. label = TRUE,
2877. label.size = 5,
2878. pt.size = 1,
2879. raster = F,
2880. cols = sorted_palopal) +
2881. NoLegend()
2883. DimPlot(seurat_sct,
2884. reduction = "umap",
2885. label = TRUE,
2886. label.size = 5,
2887. pt.size = 1,
2888. raster = F,
2889. cols = sorted_palopal,
2890. split.by = "time_point") +
2891. NoLegend()
2893. # Rod - Pde6a, Rho, Sag, Gnat1, Nrl
2894. VlnPlot(seurat_sct,
2895. features = c("Pde6a", "Rho", "Sag", "Gnat1", "Nrl",
2896. "Map2","Rom1"),
2897. stack = T,
2898. flip = T,
2899. raster = F) + NoLegend()
2901. FeaturePlot(seurat_sct,
2902. reduction = "umap",
2903. features = c("Pde6a", "Rho", "Sag", "Gnat1", "Nrl",
2904. "Map2"),
2905. order = TRUE,
2906. min.cutoff = 'q10',
2907. label = TRUE,
2908. label.size = 3,
2909. pt.size = 0.5,
2910. raster = F)
2912. # S Cone - Pde6h, Arr3, Gnat2, Opn1sw, Ccdc136, Ttr
2913. VlnPlot(seurat_sct,
2914. features = c("Pde6h", "Arr3", "Gnat2", "Opn1sw", "Ccdc136" ,"Ttr"),
2915. stack = T,
2916. flip = T,
2917. raster = F) + NoLegend()
2919. FeaturePlot(seurat_sct,
2920. reduction = "umap",
2921. features = c("Pde6h", "Arr3", "Gnat2", "Opn1mw"),
2922. order = TRUE,
2923. min.cutoff = 'q10',
2924. label = TRUE,
2925. label.size = 3,
2926. pt.size = 0.5,
2927. raster = F)
2929. FeaturePlot(seurat_sct,
2930. reduction = "umap",
2931. features = c("Ccdc136" ,"Ttr", "Vopp1", "Lmo4"),
2932. order = TRUE,
2933. min.cutoff = 'q10',
2934. label = TRUE,
2935. label.size = 3,
2936. pt.size = 0.5,
2937. raster = F)
2939. # M/L Cone - Pde6h, Arr3, Gnat2, Opn1mw/Opn1lw, Vopp1, Lmo4
2940. VlnPlot(seurat_sct,
2941. features = c("Pde6h", "Arr3", "Gnat2", "Opn1mw",
2942. "Opn1lw", "Vopp1", "Lmo4"),
2943. stack = T,
2944. flip = T,
2945. raster = F) + NoLegend()
2947. # HC - Pvalb, Lhx1, Pax6, Onecut1
2948. VlnPlot(seurat_sct,
2949. features = c("Pvalb", "Lhx1", "Pax6", "Onecut1"),
2950. stack = T,
2951. flip = T,
2952. raster = F) + NoLegend()
2954. FeaturePlot(seurat_sct,
2955. reduction = "umap",
2956. features = c("Pvalb", "Lhx1", "Pax6", "Onecut1"),
2957. order = TRUE,
2958. min.cutoff = 'q10',
2959. label = TRUE,
2960. label.size = 3,
2961. pt.size = 0.5,
2962. raster = F)
2964. # AC - Snhg11, Pax6, Slc32a1, Crabp1, Nrxn2, Gad1/Gad2, Maf, Tfap2a, Slc6a9, Ebf3
2965. VlnPlot(seurat_sct,
2966. features = c("Snhg11", "Pax6", "Slc32a1", "Crabp1",
2967. "Nrxn2", "Gad1", "Gad2", "Maf",
2968. "Tfap2a", "Slc6a9", "Ebf3"),
2969. stack = T,
2970. flip = T) + NoLegend()
2972. hirestiff(paste(prefixPCres,"UMAP","Amacrine","markers","hires.tiff", sep = "_"))
2973. lowrestiff(paste(prefixPCres,"UMAP","Amacrine","markers","lowres.tiff", sep = "_"))
2975. # RBC - Gsg1, Pax6, Slc6a9, Vsx2, Otx2, Prkca, Isl1, Grm6, Cabp5, Vstm2b, Casp7, Rpa1
2976. VlnPlot(seurat_sct,
2977. features = c("Gsg1", "Pax6", "Slc6a9", "Vsx2", "Otx2",
2978. "Prkca", "Isl1", "Grm6", "Cabp5", "Vstm2b",
2979. "Casp7", "Rpa1"),
2980. stack = T,
2981. flip = T) + NoLegend()
2983. FeaturePlot(seurat_sct,
2984. reduction = "umap",
2985. features = c("Gsg1", "Pax6", "Slc6a9", "Vsx2", "Otx2",
2986. "Prkca", "Isl1", "Grm6", "Cabp5", "Vstm2b",
2987. "Casp7", "Rpa1"),
2988. order = TRUE,
2989. min.cutoff = 'q10',
2990. label = TRUE,
2991. label.size = 3,
2992. pt.size = 0.5,
2993. raster = F)
2995. # ON CBC - Gsg1, Pax6, Slc6a9, Vsx2, Otx2, App, Scgn, Vsx1, Isl1, Grm6
2996. VlnPlot(seurat_sct,
2997. features = c("Gsg1", "Pax6", "Slc6a9", "Vsx2", "Otx2",
2998. "App", "Scgn", "Vsx1", "Isl1", "Grm6"),
2999. stack = T,
3000. flip = T) + NoLegend()
3002. FeaturePlot(seurat_sct,
3003. reduction = "umap",
3004. features = c("App", "Scgn", "Vsx1", "Isl1", "Grm6"),
3005. order = TRUE,
3006. min.cutoff = 'q10',
3007. label = TRUE,
3008. label.size = 3,
3009. pt.size = 0.5,
3010. raster = F)
3012. # OFF CBC - Gsg1, Pax6, Slc6a9, Vsx2, Otx2, App, Scgn, Vsx1, Grin2b, Grik1
3013. VlnPlot(seurat_sct,
3014. features = c("Gsg1", "Pax6", "Slc6a9", "Vsx2", "Otx2",
3015. "App", "Scgn", "Vsx1", "Grin2b", "Grik1"),
3016. stack = T,
3017. flip = T) + NoLegend()
3019. FeaturePlot(seurat_sct,
3020. reduction = "umap",
3021. features = c("App", "Scgn", "Vsx1", "Grin2b", "Grik1"),
3022. order = TRUE,
3023. min.cutoff = 'q10',
3024. label = TRUE,
3025. label.size = 3,
3026. pt.size = 0.5,
3027. raster = F)
3029. # Müller - Gsta1, Pax6, Rlbp1, Aqp4, Slc1a3, Apoe, Dkk3, Gpr37, Rax, Hes1, Notch1, Glul
3030. VlnPlot(seurat_sct,
3031. features = c("Gsta1", "Pax6", "Rlbp1", "Aqp4", "Slc1a3",
3032. "Apoe", "Dkk3", "Gpr37", "Rax", "Hes1",
3033. "Notch1", "Glul"),
3034. stack = T,
3035. flip = T,
3036. raster = F) + NoLegend()
3038. FeaturePlot(seurat_sct,
3039. reduction = "umap",
3040. features = c("Gfap", "Aqp4", "S100b"),
3041. order = TRUE,
3042. min.cutoff = 'q10',
3043. label = TRUE,
3044. label.size = 3,
3045. pt.size = 0.5,
3046. raster = F)
3048. # Microglia - C1qa, Tmem119, Cx3cr1, Apoe, Pax6, P2ry12, Aif1
3049. VlnPlot(seurat_sct,
3050. features = c("C1qa", "Tmem119", "Cx3cr1", "Apoe", "Pax6",
3051. "P2ry12", "Aif1"),
3052. stack = T,
3053. flip = T,
3054. raster = F) + NoLegend()
3056. FeaturePlot(seurat_sct,
3057. reduction = "umap",
3058. features = c("C1qa", "Tmem119", "Cx3cr1",
3059. "P2ry12", "Aif1"),
3060. order = TRUE,
3061. min.cutoff = 'q10',
3062. label = TRUE,
3063. label.size = 3,
3064. pt.size = 0.5,
3065. raster = F)
3067. # EC - Vwf, Pecam1, Cldn5, Cdh5, Tek, Kdr, Flt1
3068. VlnPlot(seurat_sct,
3069. features = c("Vwf", "Pecam1", "Cldn5", "Cdh5", "Tek",
3070. "Kdr", "Flt1", "Rbpms"),
3071. stack = T,
3072. flip = T,
3073. raster = F) + NoLegend()
3075. # Pericyte - Rgs5, Kcnj8, Myl9, Cspg4, Pdgfrb, Myh11, Acta2
3076. VlnPlot(seurat_sct,
3077. features = c("Rgs5","Kcnj8","Myl9","Cspg4","Pdgfrb",
3078. "Myh11","Acta2","Mgp", "Crip1", "Tpm2",
3079. "Rbpms","Pgf","Igf1","Igf2r","Vegfa",
3080. "Axl"),
3081. stack = T,
3082. flip = T,
3083. raster = F) + NoLegend()
3085. VlnPlot(seurat_sct,
3086. features = c("Dnmbp","Tubb3","Nefl","Nefm","Nefh","Rbfox3","Map2","Syp","Gap43"),
3087. stack = T,
3088. flip = T,
3089. raster = F) + NoLegend()
3091. FeaturePlot(seurat_sct,
3092. reduction = "umap",
3093. features = c("Dnmbp","Tubb3","Nefl","Nefm","Nefh",
3094. "Rbfox3","Map2","Syp","Gap43"),
3095. order = TRUE,
3096. min.cutoff = 'q10',
3097. label = TRUE,
3098. label.size = 3,
3099. pt.size = 0.5,
3100. raster = F)
3102. VlnPlot(seurat_sct,
3103. features = c("Mgp", "Crip1", "Tpm2"),
3104. stack = T,
3105. flip = T,
3106. raster = F) + NoLegend()
3109. # Macrophage - Cxcr4, Cd53, Ptprc
3110. VlnPlot(seurat_sct,
3111. features = c("Cxcr4", "Cd53", "Ptprc"),
3112. stack = T,
3113. flip = T,
3114. raster = F) + NoLegend()
3116. FeaturePlot(seurat_sct,
3117. reduction = "umap",
3118. features = c("Cxcr4", "Cd53", "Ptprc"),
3119. order = TRUE,
3120. min.cutoff = 'q10',
3121. label = TRUE,
3122. label.size = 3,
3123. pt.size = 0.5,
3124. raster = F)
3126. # RGC - Rbpms
3127. VlnPlot(seurat_sct,
3128. features = c("Rbpms","Pou4f1","Pou4f2","Pou4f3","Thy1",
3129. "Map2","Opn4","Slc17a6","Jam2","Tbr1"),
3130. stack = T,
3131. flip = T,
3132. raster = F) + NoLegend()
3134. VlnPlot(seurat_sct,
3135. features = c("Rbpms","Pou4f1","Pou4f2",
3136. "Jam2","Foxp1","Foxp2"),
3137. stack = T,
3138. flip = T,
3139. raster = F) + NoLegend()
3141. # Ctxn3+Müller - Cd9, Penk, Ctxn3
3142. VlnPlot(seurat_sct,
3143. features = c("Cd9", "Penk", "Ctxn3"),
3144. stack = T,
3145. flip = T,
3146. raster = F) + NoLegend()
3148. # Ctxn3-Müller - Cd9, Penk
3149. VlnPlot(seurat_sct,
3150. features = c("Cd9", "Penk"),
3151. stack = T,
3152. flip = T,
3153. raster = F) + NoLegend()
3155. # EC - Pltp, Csrp2, Cxcl12, Id1, Ramp2, Slc2a1
3156. VlnPlot(seurat_sct,
3157. features = c("Pltp", "Csrp2", "Cxcl12", "Id1",
3158. "Ramp2", "Slc2a1"),
3159. stack = T,
3160. flip = T,
3161. raster = F) + NoLegend()
3163. # Pericyte - Mgp, Crip1, Tpm2
3164. VlnPlot(seurat_sct,
3165. features = c("Mgp", "Crip1", "Tpm2"),
3166. stack = T,
3167. flip = T,
3168. raster = F) + NoLegend()
3170. #### Genes from 9/18 ####
3172. # Egr11 in rods and cones - probably Egr1
3173. # Amacrine: Snhg11, Pax6, Slc32a1, Crabp1, Nrxn2,
3174. # Gad1/Gad2, Maf, Tfap2a, Slc6a9, Ebf3 - done above
3175. # Muller glia and astrocytes:
3176. # glutamine synthetase (GLUL), clusterin (CLU),
3177. # and apolipoprotein E (APOE).
3178. # NxnL1: metabolic dysfunction-glucosse metabolism.
3179. # Hk2: a key enzyme for aerobic glycolysis and required for survival of photoreceptors during aging.
3180. # Cebpd: upregulation is related to neuroprotection.
3182. VlnPlot(seurat_sct,
3183. features = c("Egr1","Nxnl1","Hk2","Cebpd"),
3184. stack = T,
3185. flip = T) + NoLegend()
3187. hirestiff(paste(prefixPCres,"VLN","Egr1","Nxnl1","Hk2","Cebpd","hires.tiff", sep = "_"))
3188. lowrestiff(paste(prefixPCres,"VLN","Egr1","Nxnl1","Hk2","Cebpd","lowres.tiff", sep = "_"))
3190. VlnPlot(seurat_sct,
3191. features = c("Egr1"),
3192. flip = T,
3193. idents = c(8,9),
3194. split.by = "time_point_treatment")
3196. VlnPlot(seurat_sct,
3197. features = c("Egr1"),
3198. flip = T,
3199. idents = c("ROD","CONE"),
3200. split.by = "time_point_treatment")
3202. hirestiff(paste(prefixPCres,"VLN","Egr1","rod-cone","hires.tiff", sep = "_"))
3203. lowrestiff(paste(prefixPCres,"VLN","Egr1","rod-cone","lowres.tiff", sep = "_"))
3205. VlnPlot(seurat_sct,
3206. features = c("Nxnl1"),
3207. flip = T,
3208. idents = c("ROD","CONE"),
3209. split.by = "time_point_treatment")
3211. hirestiff(paste(prefixPCres,"VLN","Nxnl1","rod-cone","hires.tiff", sep = "_"))
3212. lowrestiff(paste(prefixPCres,"VLN","Nxnl1","rod-cone","lowres.tiff", sep = "_"))
3214. VlnPlot(seurat_sct,
3215. features = c("Hk2"),
3216. flip = T,
3217. idents = c("ROD","CONE"),
3218. split.by = "time_point_treatment")
3220. hirestiff(paste(prefixPCres,"VLN","Hk2","rod-cone","hires.tiff", sep = "_"))
3221. lowrestiff(paste(prefixPCres,"VLN","Hk2","rod-cone","lowres.tiff", sep = "_"))
3223. VlnPlot(seurat_sct,
3224. features = c("Cebpd"),
3225. flip = T,
3226. idents = c("GLN1","GLN2","MULG"),
3227. split.by = "time_point_treatment")
3229. hirestiff(paste(prefixPCres,"VLN","Cebpd","gln_mulg","hires.tiff", sep = "_"))
3230. lowrestiff(paste(prefixPCres,"VLN","Cebpd","gln_mulg","lowres.tiff", sep = "_"))
3232. VlnPlot(seurat_sct,
3233. features = c("Glul", "Clu", "Apoe"),
3234. stack = T,
3235. flip = T) + NoLegend()
3237. hirestiff(paste(prefixPCres,"VLN","Muller","markers","09_18","hires.tiff", sep = "_"))
3238. lowrestiff(paste(prefixPCres,"VLN","Muller","markers","09_18","lowres.tiff", sep = "_"))
3241. #### Genes from 9/28 ####
3242. # VlnPlot(seurat_sct,
3243. # features = c("Pou4f1", "Pou4f2", "Pou4f3", "Sncg", "Isl1", "Nefl", "Nefm",
3244. # "Elavl4", "Eomes", "Calb1", "Calb2", "Opn4", "Fes", "Tpbg", "Irx3",
3245. # "Irx4", "Trb1", "Neurod2", "Penk", "Jam2", "Calca", "Spp1", "Plpp4",
3246. # "Coch", "Foxp1", "Foxp2", "Tusc5", "Satb2", "Cdk15", "Anxa3", "Slc17a6",
3247. # "Tbx20", "Serpine2", "Nmb", "Adcyap1", "Nmb", "Cartpt", "Mmp17", "Gpr88",
3248. # "Fam19a4", "Ebf3", "Col25a1", "A230065H16Rik", "Dcx"),
3249. # stack = T,
3250. # flip = T) + NoLegend()
3252. VlnPlot(seurat_sct,
3253. features = c("Pou4f1","Pou4f2","Pou4f3","Sncg","Isl1",
3254. "Nefl","Nefm","Elavl4",#"Eomes",
3255. "Calb1",
3256. "Calb2","Opn4","Fes","Tpbg","Irx3"),
3257. stack = T,
3258. flip = T,
3259. idents = c(6,15,20)) + NoLegend() + ggtitle("RGC set 1")
3261. hirestiff(paste(prefixPCres,"VLN","RGC","set","1","sept_28","hires.tiff", sep = "_"))
3262. lowrestiff(paste(prefixPCres,"VLN","RGC","set","1","sept_28","lowres.tiff", sep = "_"))
3264. VlnPlot(seurat_sct,
3265. features = c("Irx4","Trib1","Neurod2","Penk","Jam2",
3266. "Calca","Spp1","Plpp4","Coch","Foxp1",
3267. "Foxp2","Trarg1","Satb2","Cdk15","Anxa3"),
3268. stack = T,
3269. flip = T,
3270. idents = c(6,15,20)) + NoLegend() + ggtitle("RGC set 2")
3272. hirestiff(paste(prefixPCres,"VLN","RGC","set","2","sept_28","hires.tiff", sep = "_"))
3273. lowrestiff(paste(prefixPCres,"VLN","RGC","set","2","sept_28","lowres.tiff", sep = "_"))
3275. VlnPlot(seurat_sct,
3276. features = c("Slc17a6","Tbx20","Serpine2","Nmb","Adcyap1",
3277. "Nmb","Cartpt","Mmp17","Gpr88","Tafa4",
3278. "Ebf3","Col25a1","Lbhd2","Dcx"),
3279. stack = T,
3280. flip = T,
3281. idents = c(6,15,20)) + NoLegend() + ggtitle("RGC set 3")
3283. hirestiff(paste(prefixPCres,"VLN","RGC","set","3","sept_28","hires.tiff", sep = "_"))
3284. lowrestiff(paste(prefixPCres,"VLN","RGC","set","3","sept_28","lowres.tiff", sep = "_"))
3286. # VlnPlot(seurat_sct,
3287. # features = c("Gldn","Chl1","Qpct","Dgkg","Timp2",
3288. # "Junb","Ifitm10","Fgf1","Sema5a","Ramp1",
3289. # "Ostf1","Esrrg","Etl4","Kbtbd11","Ctxn3",
3290. # "Igf11","Sdk1","Sdk2","Man1a","Ifi27"),
3291. # stack = T,
3292. # flip = T) + NoLegend()
3294. VlnPlot(seurat_sct,
3295. features = c("Gldn","Chl1","Qpct","Dgkg","Timp2",
3296. "Junb","Ifitm10","Fgf1","Sema5a","Ramp1"),
3297. stack = T,
3298. flip = T,
3299. idents = c(6,15,20)) + NoLegend() + ggtitle("RGC set 4")
3301. hirestiff(paste(prefixPCres,"VLN","RGC","set","4","sept_28","hires.tiff", sep = "_"))
3302. lowrestiff(paste(prefixPCres,"VLN","RGC","set","4","sept_28","lowres.tiff", sep = "_"))
3304. VlnPlot(seurat_sct,
3305. features = c("Ostf1","Esrrg","Etl4","Kbtbd11","Ctxn3",
3306. "Igf2","Sdk1","Sdk2","Man1a1","Ifi27"),
3307. stack = T,
3308. flip = T,
3309. idents = c(6,15,20)) + NoLegend() + ggtitle("RGC set 5")
3311. hirestiff(paste(prefixPCres,"VLN","RGC","set","5","sept_28","hires.tiff", sep = "_"))
3312. lowrestiff(paste(prefixPCres,"VLN","RGC","set","5","sept_28","lowres.tiff", sep = "_"))
3314. VlnPlot(seurat_sct,
3315. features = c("Pdlim1","Col1a1","Pdgfrb","Rbpms"),
3316. stack = T,
3317. flip = T,
3318. idents = c(6,15,20)) + NoLegend() + ggtitle("Pericytes set 1")
3320. hirestiff(paste(prefixPCres,"VLN","Pericytes","set","1","sept_28","hires.tiff", sep = "_"))
3321. lowrestiff(paste(prefixPCres,"VLN","Pericytes","set","1","sept_28","lowres.tiff", sep = "_"))
3323. VlnPlot(seurat_sct,
3324. features = c("Rgs5","Kcnj8","Myl9","Cspg4","Pdgfrb",
3325. "Myh11","Acta2","Mgp","Crip1","Tpm2",
3326. "Pdlim1","Col1a1"),
3327. stack = T,
3328. flip = T,
3329. idents = c(6,15,20)) + NoLegend() + ggtitle("Pericytes set 1 with Wang et al markers")
3331. hirestiff(paste(prefixPCres,"VLN","Pericytes","with","wang","markers","sept_28","hires.tiff", sep = "_"))
3332. lowrestiff(paste(prefixPCres,"VLN","Pericytes","with","wang","markers","sept_28","lowres.tiff", sep = "_"))
3335. VlnPlot(seurat_sct,
3336. features = c("Prph","Ebf3","Atp13a5",
3337. "Cend1","Esrrb","Uchl1"),
3338. stack = T,
3339. flip = T,
3340. idents = c(6,15,20)) +
3341. NoLegend() +
3342. ggtitle("New markers 10-6-23")
3344. hirestiff(paste(prefixPCres,"VLN","new","markers","Oct6","hires.tiff", sep = "_"))
3345. lowrestiff(paste(prefixPCres,"VLN","new","markers","Oct6","lowres.tiff", sep = "_"))
3347. Prph, Ebf3
3348. Pericyte: Atp13a5
3350. BM88(Cend1), ERRbeta (Esrrb), PGP9.5 (Uchl1)
3353. FeaturePlot(seurat_sct,
3354. reduction = "umap",
3355. features = c("Rbpms"),
3356. order = TRUE,
3357. min.cutoff = 'q10',
3358. label = TRUE,
3359. label.size = 3,
3360. pt.size = 0.5,
3361. raster = F,
3362. split.by = "time_point")
3364. FeaturePlot(seurat_sct,
3365. reduction = "umap",
3366. features = c("Sncg","Rbpms","Pou4f1","Kcnj8"),
3367. order = TRUE,
3368. min.cutoff = 'q10',
3369. label = TRUE,
3370. label.size = 3,
3371. pt.size = 0.5,
3372. raster = F)
3374. FeaturePlot(seurat_sct,
3375. reduction = "umap",
3376. features = c("Sncg","Rbpms","Isl1","Rbfox3"),
3377. order = TRUE,
3378. min.cutoff = 'q10',
3379. label = TRUE,
3380. label.size = 3,
3381. pt.size = 0.5,
3382. raster = F)
3384. FeaturePlot(seurat_sct,
3385. reduction = "umap",
3386. features = c("Sncg","Rbpms","Pou4f1","Pou4f2"),
3387. order = TRUE,
3388. min.cutoff = 'q10',
3389. label = TRUE,
3390. label.size = 3,
3391. pt.size = 0.5,
3392. raster = F)
3394. #### Find Cluster Markers for all cells ####
3396. # find all markers
3397. cluster_markers <- FindAllMarkers(seurat_sct, only.pos = TRUE, min.pct = 0.50, logfc.threshold = 0.4)
3399. write.csv(cluster_markers,
3400. file = paste(prefixPCres,"clustermarkers_pos_only_minpct0.50_logfcthresh0.4.csv", sep = "_"))
3402. #### find DEGs between treatments ####
3404. Idents(object = seurat_sct) <- "treatment"
3406. cluster_markers <- FindAllMarkers(seurat_sct, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.4)
3408. write.csv(cluster_markers, file =
3409. paste(prefixPCres,"degs_treatment_all_cells_minpct0.25_logfcthresh0.4.csv", sep = "_"))
3411. Idents(object = seurat_sct) <- ident_res
3413. #### find DEGs between time points ####
3415. Idents(object = seurat_sct) <- "time_point"
3417. cluster_markers <- FindAllMarkers(seurat_sct, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.4)
3419. write.csv(cluster_markers, file =
3420. paste(prefixPCres,"degs_timepoint_all_cells_minpct0.25_logfcthresh0.4.csv", sep = "_"))
3422. Idents(object = seurat_sct) <- ident_res
3424. #### subsets for further DEG analysis ####
3426. # subset seurat obj by time point for treated vs untreated DEGs
3427. Idents(object = seurat_sct) <- "time_point"
3429. p60seurat <- subset(seurat_sct, idents = "p60")
3430. p90seurat <- subset(seurat_sct, idents = "p90")
3432. Idents(object = p60seurat) <- "treatment"
3433. cluster_markers <- FindAllMarkers(p60seurat, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.4)
3435. write.csv(cluster_markers, file =
3436. paste(prefixPCres,"degs_treatment_p60_minpct0.25_logfcthresh0.4.csv", sep = "_"))
3438. Idents(object = p90seurat) <- "treatment"
3439. cluster_markers <- FindAllMarkers(p90seurat, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.4)
3441. write.csv(cluster_markers, file =
3442. paste(prefixPCres,"degs_treatment_p90_minpct0.25_logfcthresh0.4.csv", sep = "_"))
3444. # subset seurat obj by treatment for p60 vs p90 DEGs
3445. Idents(object = seurat_sct) <- "treatment"
3447. TRseurat <- subset(seurat_sct, idents = "Cell_Treated")
3448. UTseurat <- subset(seurat_sct, idents = "Untreated")
3450. Idents(object = TRseurat) <- "time_point"
3451. cluster_markers <- FindAllMarkers(TRseurat, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.4)
3453. write.csv(cluster_markers, file =
3454. paste(prefixPCres,"degs_treated_p60_vs_p90_minpct0.25_logfcthresh0.4.csv", sep = "_"))
3456. Idents(object = UTseurat) <- "time_point"
3457. cluster_markers <- FindAllMarkers(UTseurat, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.4)
3459. write.csv(cluster_markers, file =
3460. paste(prefixPCres,"degs_untreated_p60_vs_p90_minpct0.25_logfcthresh0.4.csv", sep = "_"))
3462. Idents(object = seurat_sct) <- ident_res
3463. Idents(object = p60seurat) <- ident_res
3464. Idents(object = p90seurat) <- ident_res
3465. Idents(object = TRseurat) <- ident_res
3466. Idents(object = UTseurat) <- ident_res
3468. #### save subset objects ####
3469. saveRDS(p60seurat, paste0(prefixPCres,"_p60_seurat.rds"))
3471. # p60seurat <- readRDS(file = paste0(prefixPCres,"_p60_seurat.rds"))
3473. saveRDS(p90seurat, paste0(prefixPCres,"_p90_seurat.rds"))
3475. p90seurat <- readRDS(file = paste0(prefixPCres,"_p90_seurat.rds"))
3477. saveRDS(TRseurat, paste0(prefixPCres,"_cell_treated_seurat.rds"))
3479. # TRseurat <- readRDS(file = paste0(prefixPCres,"_cell_treated_seurat.rds"))
3481. saveRDS(UTseurat, paste0(prefixPCres,"_untreated_seurat.rds"))
3483. # UTseurat <- readRDS(file = paste0(prefixPCres,"_untreated_seurat.rds"))
3485. #### for loop for DEGs within clusters ####
3487. clusterlist <- as.list(0:(max(as.integer([email hidden]$integrated_snn_res.0.2))-1))
3488. clusterlist
3490. for (k in 1:length(clusterlist)){
3492. nam <- paste("clusterDEG_treatment", clusterlist[k], sep = "_")
3493. assign(nam, FindMarkers(seurat_sct, ident.1 = "Cell_Treated", group.by = "treatment", subset.ident = clusterlist[k]))
3495. write.csv(get(nam), file = paste(prefixPCres, nam, "list.csv", sep = "_"), row.names = TRUE)
3497. }
3499. # positive LFC means upregulation in Cell_Treated cells
3501. for (k in 1:length(clusterlist)){
3503. nam <- paste("clusterDEG_timepoint", clusterlist[k], sep = "_")
3504. assign(nam, FindMarkers(seurat_sct, ident.1 = "p60", group.by = "time_point", subset.ident = clusterlist[k]))
3506. write.csv(get(nam), file = paste(prefixPCres, nam, "list.csv", sep = "_"), row.names = TRUE)
3508. }
3510. # positive LFC means upregulation in p60 cells
3512. #### LABEL CELL TYPES ####
3513. #### Label cell types (short name) ####
3514. seurat_sct <- RenameIdents(object = seurat_sct,
3515. '0' = "RBC",
3516. '1' = "RBC",
3517. '2' = "ROD",
3518. '3' = "MG",
3519. '4' = "CBC",
3520. '5' = "RBC",
3521. '6' = "RGC",
3522. '7' = "EC",
3523. '8' = "ROD",
3524. '9' = "MULG",
3525. '10' = "MULG",
3526. '11' = "CONE",
3527. '12' = "MG",
3528. '13' = "MG",
3529. '14' = "CBC",
3530. '15' = "RGC",
3531. '16' = "AC",
3532. '17' = "MULG",
3533. '18' = "MULG",
3534. '19' = "MULG",
3535. '20' = "RGC",
3536. '21' = "MG",
3537. '22' = "EC",
3538. '23' = "M0",
3539. '24' = "MULG",
3540. '25' = "ROD",
3541. '26' = "EC",
3542. '27' = "MULG",
3543. '28' = "MG")
3545. #### save cell type (short name) as a metadata column ####
3546. seurat_sct[["Cell.Type"]] <- [email hidden]
3548. #### dimplots of cell type (short name) labeled clusters ####
3550. DimPlot(seurat_sct,
3551. reduction = "umap",
3552. label = TRUE,
3553. label.size = 6,
3554. pt.size = 1,
3555. raster = F) +
3556. NoLegend()
3557. hirestiff(paste(prefixPCres,"UMAP","celltypes","hires.tiff",sep = "_"))
3558. lowrestiff(paste(prefixPCres,"UMAP","celltypes","lowres.tiff",sep = "_"))
3560. DimPlot(seurat_sct,
3561. reduction = "umap",
3562. label = TRUE,
3563. label.size = 5,
3564. pt.size = 1,
3565. raster = F,
3566. split.by = "time_point") +
3567. NoLegend()
3568. hirestiff(paste(prefixPCres,"UMAP","celltypes","by","time_point","hires.tiff",sep = "_"))
3569. lowrestiff(paste(prefixPCres,"UMAP","celltypes","by","time_point","lowres.tiff",sep = "_"))
3571. DimPlot(seurat_sct,
3572. reduction = "umap",
3573. label = TRUE,
3574. label.size = 3,
3575. pt.size = 1,
3576. raster = F,
3577. split.by = "time_point_treatment") +
3578. NoLegend()
3579. hirestiff(paste(prefixPCres,"UMAP","celltypes","by","time_point_treatment","hires.tiff",sep = "_"))
3580. lowrestiff(paste(prefixPCres,"UMAP","celltypes","by","time_point_treatment","lowres.tiff",sep = "_"))
3582. for (i in 1:length(levels(seurat_sct$time_point_treatment))){
3583. p <- (DimPlot(seurat_sct,
3584. reduction = "umap",
3585. label = TRUE,
3586. label.size = 8,
3587. pt.size = 1,
3588. cells = c(WhichCells(seurat_sct,
3589. expression =
3590. time_point_treatment ==
3591. levels(seurat_sct$time_point_treatment)[i]))) +
3592. NoLegend() +
3593. ggtitle(paste0(levels(seurat_sct$time_point_treatment)[i])))
3595. print(p)
3597. hirestiff(paste(prefixPCres,"UMAP","split_by","time_point_treatment","with_cluster_labels",
3598. levels(seurat_sct$time_point_treatment)[i],"hires.tiff",sep = "_"))
3599. lowrestiff(paste(prefixPCres,"UMAP","split_by","time_point_treatment","with_cluster_labels",
3600. levels(seurat_sct$time_point_treatment)[i],"lowres.tiff",sep = "_"))
3601. hirestiffsquare(paste(prefixPCres,"UMAP","split_by","time_point_treatment","with_cluster_labels",
3602. levels(seurat_sct$time_point_treatment)[i],"hires_square.tiff",sep = "_"))
3603. lowrestiffsquare(paste(prefixPCres,"UMAP","split_by","time_point_treatment","with_cluster_labels",
3604. levels(seurat_sct$time_point_treatment)[i],"lowres_square.tiff",sep = "_"))
3605. }
3607. rm(p)
3608. #x#### Label cell types (long name) ####
3610. seurat_sct <- RenameIdents(object = seurat_sct,
3611. 'BIP1' = "Bipolar_1",
3612. 'BIP2' = "Bipolar_2",
3613. 'GLN1' = "Ganglion_1",
3614. 'GLN2' = "Ganglion_2",
3615. 'MG' = "Microglia",
3616. 'ROD' = "Rods",
3617. 'CONE' = "Cones",
3618. 'MULG' = "Müller_Glia",
3619. 'AMCR' = "Amacrine",
3620. 'RTIC' = "Reticulocytes")
3622. #x#### save cell type (long name) as a metadata column ####
3623. seurat_sct[["Cell.Type.Long"]] <- [email hidden]
3625. #x#### Label cell types plus cluster number ####
3626. seurat_sct <- RenameIdents(object = seurat_sct,
3627. '0' = "RGC0",
3628. '1' = "BP1",
3629. '2' = "BP2",
3630. '3' = "BP_ROD3",
3631. '4' = "FIB4",
3632. '5' = "MG5",
3633. '6' = "FIB6",
3634. '7' = "MG7",
3635. '8' = "ROD8",
3636. '9' = "CONE9",
3637. '10' = "MULG10",
3638. '11' = "FIB_MULG_PGMT11",
3639. '12' = "FIB12",
3640. '13' = "MG13",
3641. '14' = "MULG14",
3642. '15' = "FIB_ASTR15",
3643. '16' = "FIB_PGMT16",
3644. '17' = "RGC_AMCR17",
3645. '18' = "FIB18",
3646. '19' = "FIB19",
3647. '20' = "FIB_MEL20",
3648. '21' = "MG21",
3649. '22' = "MG22",
3650. '23' = "BLOOD23",
3651. '24' = "FIB24",
3652. '25' = "FIB25",
3653. '26' = "MG26",
3654. '27' = "MULG27",
3655. '28' = "MG28",
3656. '29' = "MG29")
3658. #### Label cell types
3659. seurat_sct <- RenameIdents(object = seurat_sct,
3660. '0' = "BPC_0",
3661. '1' = "BPC_1",
3662. '2' = "BPC_2",
3663. '3' = "BPC_3",
3664. '4' = "RGC_4",
3665. '5' = "MG_5",
3666. '6' = "RGC_6",
3667. '7' = "MG_7",
3668. '8' = "ROD_8",
3669. '9' = "CONE_9",
3670. '10' = "MULG_10",
3671. '11' = "MULG_11",
3672. '12' = "RGC_12",
3673. '13' = "MG_13",
3674. '14' = "MULG_14",
3675. '15' = "MULG_15",
3676. '16' = "PGMT_16",
3677. '17' = "AMCR_17",
3678. '18' = "RGC_18",
3679. '19' = "RGC_19",
3680. '20' = "MEL_20",
3681. '21' = "MG_21",
3682. '22' = "MG_22",
3683. '23' = "RTIC_23",
3684. '24' = "RGC_24",
3685. '25' = "RGC_25",
3686. '26' = "MG_26",
3687. '27' = "MULG_27",
3688. '28' = "MG_28",
3689. '29' = "MG_29")
3691. #x#### save cell type plus cluster numbers as a metadata column ####
3692. seurat_sct[["Cluster.Type"]] <- [email hidden]
3693. seurat_sct[["Type.Cluster"]] <- [email hidden]
3695. saveRDS(seurat_sct, paste0(prefixPCres,"_seurat_after_RNAnorm_NEW_CELLTYPES.rds"))
3696. seurat_sct <- readRDS(file = paste0(prefixPCres,"_seurat_after_RNAnorm_NEW_CELLTYPES.rds"))
3698. #x#### dimplots of cell type (long name) labeled clusters ####
3700. DimPlot(seurat_sct,
3701. reduction = "umap",
3702. label = TRUE,
3703. label.size = 8,
3704. pt.size = 1) +
3705. NoLegend()
3706. hirestiff(paste(prefixPCres,"UMAP","celltypesLN","hires.tiff",sep = "_"))
3707. lowrestiff(paste(prefixPCres,"UMAP","celltypesLN","lowres.tiff",sep = "_"))
3709. DimPlot(seurat_sct,
3710. reduction = "umap",
3711. label = TRUE,
3712. label.size = 6,
3713. pt.size = 1,
3714. split.by = "time_point") +
3715. NoLegend()
3716. hirestiff(paste(prefixPCres,"UMAP","celltypesLN","by","time_point","hires.tiff",sep = "_"))
3717. lowrestiff(paste(prefixPCres,"UMAP","celltypesLN","by","time_point","lowres.tiff",sep = "_"))
3719. DimPlot(seurat_sct,
3720. reduction = "umap",
3721. label = TRUE,
3722. label.size = 6,
3723. pt.size = 1,
3724. split.by = "treatment") +
3725. NoLegend()
3726. hirestiff(paste(prefixPCres,"UMAP","celltypesLN","by","treatment","hires.tiff",sep = "_"))
3727. lowrestiff(paste(prefixPCres,"UMAP","celltypesLN","by","treatment","lowres.tiff",sep = "_"))
3729. DimPlot(seurat_sct,
3730. reduction = "umap",
3731. label = TRUE,
3732. label.size = 4,
3733. pt.size = 1,
3734. split.by = "time_point_treatment") +
3735. NoLegend()
3736. hirestiff(paste(prefixPCres,"UMAP","celltypesLN","by","time_point_treatment","hires.tiff",sep = "_"))
3737. lowrestiff(paste(prefixPCres,"UMAP","celltypesLN","by","time_point_treatment","lowres.tiff",sep = "_"))
3739. # set cluster colors
3740. clust_cols <- scales::hue_pal()(length(levels(seurat_sct)))
3741. names(clust_cols) <- levels(seurat_sct)
3743. levels(seurat_sct)
3745. # UMAP of cells in each cluster by time_point_treatment without cluster labels
3746. DimPlot(seurat_sct,
3747. reduction = "umap",
3748. label = FALSE,
3749. cols = clust_cols,
3750. split.by = "time_point_treatment",
3751. pt.size = 1) +
3752. NoLegend() +
3753. ggtitle(paste0(ident_res))
3754. hirestiff(paste(prefixPCres,"UMAP","split_by","time_point_treatment","without_cluster_labels","hires.tiff",sep = "_"))
3755. lowrestiff(paste(prefixPCres,"UMAP","split_by","time_point_treatment","without_cluster_labels","lowres.tiff",sep = "_"))
3757. for (i in 1:length(levels(seurat_sct$time_point_treatment))){
3758. p <- (DimPlot(seurat_sct,
3759. reduction = "umap",
3760. label = FALSE,
3761. pt.size = 1,
3762. cols = clust_cols,
3763. cells = c(WhichCells(seurat_sct, expression = time_point_treatment == levels(seurat_sct$time_point_treatment)[i]))) +
3764. NoLegend() +
3765. ggtitle(paste0(levels(seurat_sct$time_point_treatment)[i])))
3767. print(p)
3769. hirestiff(paste(prefixPCres,"UMAP","split_by","time_point_treatment","no_cluster_labels",
3770. levels(seurat_sct$time_point_treatment)[i],"hires.tiff",sep = "_"))
3771. lowrestiff(paste(prefixPCres,"UMAP","split_by","time_point_treatment","no_cluster_labels",
3772. levels(seurat_sct$time_point_treatment)[i],"lowres.tiff",sep = "_"))
3773. hirestiffsquare(paste(prefixPCres,"UMAP","split_by","time_point_treatment","no_cluster_labels",
3774. levels(seurat_sct$time_point_treatment)[i],"hires_square.tiff",sep = "_"))
3775. lowrestiffsquare(paste(prefixPCres,"UMAP","split_by","time_point_treatment","no_cluster_labels",
3776. levels(seurat_sct$time_point_treatment)[i],"lowres_square.tiff",sep = "_"))
3777. }
3779. # UMAP of cells in each cluster by time_point_treatment with cluster labels
3780. DimPlot(seurat_sct,
3781. reduction = "umap",
3782. label = TRUE,
3783. cols = clust_cols,
3784. split.by = "time_point_treatment",
3785. label.size = 5,
3786. pt.size = 1) +
3787. NoLegend() +
3788. ggtitle(paste0(ident_res))
3789. hirestiff(paste(prefixPCres,"UMAP","split_by","time_point_treatment","with_cluster_labels","hires.tiff",sep = "_"))
3790. lowrestiff(paste(prefixPCres,"UMAP","split_by","time_point_treatment","with_cluster_labels","lowres.tiff",sep = "_"))
3792. for (i in 1:length(levels(seurat_sct$time_point_treatment))){
3793. p <- (DimPlot(seurat_sct,
3794. reduction = "umap",
3795. label = TRUE,
3796. label.size = 8,
3797. pt.size = 1,
3798. cols = clust_cols,
3799. cells = c(WhichCells(seurat_sct, expression = time_point_treatment == levels(seurat_sct$time_point_treatment)[i]))) +
3800. NoLegend() +
3801. ggtitle(paste0(levels(seurat_sct$time_point_treatment)[i])))
3803. print(p)
3805. hirestiff(paste(prefixPCres,"UMAP","split_by","time_point_treatment","with_cluster_labels",
3806. levels(seurat_sct$time_point_treatment)[i],"hires.tiff",sep = "_"))
3807. lowrestiff(paste(prefixPCres,"UMAP","split_by","time_point_treatment","with_cluster_labels",
3808. levels(seurat_sct$time_point_treatment)[i],"lowres.tiff",sep = "_"))
3809. hirestiffsquare(paste(prefixPCres,"UMAP","split_by","time_point_treatment","with_cluster_labels",
3810. levels(seurat_sct$time_point_treatment)[i],"hires_square.tiff",sep = "_"))
3811. lowrestiffsquare(paste(prefixPCres,"UMAP","split_by","time_point_treatment","with_cluster_labels",
3812. levels(seurat_sct$time_point_treatment)[i],"lowres_square.tiff",sep = "_"))
3813. }
3815. # UMAP of cells in each cluster by time_point_treatment with no cluster labels plus legend
3816. DimPlot(seurat_sct,
3817. reduction = "umap",
3818. label = FALSE,
3819. cols = clust_cols,
3820. label.size = 5,
3821. pt.size = 1) +
3822. ggtitle(paste0(ident_res))
3823. hirestiff(paste(prefixPCres,"UMAP","with_no_labels_plus_legend","hires.tiff",sep = "_"))
3824. lowrestiff(paste(prefixPCres,"UMAP","with_no_labels_plus_legend","lowres.tiff",sep = "_"))
3826. DimPlot(seurat_sct,
3827. reduction = "umap",
3828. label = FALSE,
3829. cols = clust_cols,
3830. split.by = "time_point_treatment",
3831. label.size = 5,
3832. pt.size = 1) +
3833. ggtitle(paste0(ident_res))
3834. hirestiff(paste(prefixPCres,"UMAP","split_by","time_point_treatment","with_no_labels_plus_legend","hires.tiff",sep = "_"))
3835. lowrestiff(paste(prefixPCres,"UMAP","split_by","time_point_treatment","with_no_labels_plus_legend","lowres.tiff",sep = "_"))
3837. for (i in 1:length(levels(seurat_sct$time_point_treatment))){
3838. p <- (DimPlot(seurat_sct,
3839. reduction = "umap",
3840. label = FALSE,
3841. label.size = 8,
3842. pt.size = 1,
3843. cols = clust_cols,
3844. cells = c(WhichCells(seurat_sct, expression = time_point_treatment == levels(seurat_sct$time_point_treatment)[i]))) +
3845. ggtitle(paste0(levels(seurat_sct$time_point_treatment)[i])))
3847. print(p)
3849. hirestiff(paste(prefixPCres,"UMAP","split_by","time_point_treatment","with_no_labels_plus_legend",
3850. levels(seurat_sct$time_point_treatment)[i],"hires.tiff",sep = "_"))
3851. lowrestiff(paste(prefixPCres,"UMAP","split_by","time_point_treatment","with_no_labels_plus_legend",
3852. levels(seurat_sct$time_point_treatment)[i],"lowres.tiff",sep = "_"))
3853. hirestiffsquare(paste(prefixPCres,"UMAP","split_by","time_point_treatment","with_no_labels_plus_legend",
3854. levels(seurat_sct$time_point_treatment)[i],"hires_square.tiff",sep = "_"))
3855. lowrestiffsquare(paste(prefixPCres,"UMAP","split_by","time_point_treatment","with_no_labels_plus_legend",
3856. levels(seurat_sct$time_point_treatment)[i],"lowres_square.tiff",sep = "_"))
3857. }
3859. rm(p,i)
3863. #### Extract number of cells per cluster per orig.ident ####
3864. n_cells <- FetchData(seurat_sct, vars = c("ident", "orig.ident")) %>%
3865. dplyr::count(ident, orig.ident) %>%
3866. tidyr::spread(ident, n)
3867. n_cells
3868. write.csv(n_cells, file = paste0(prefixPCres,"_celltypes_per_sample.csv"))
3870. rm(n_cells)
3872. #### switch between cluster numbers and cell type idents ####
3874. Idents(seurat_sct) <- ident_res
3876. Idents(seurat_sct) <- "Cell.Type"
3878. # Idents(seurat_sct) <- "Cell.Type.Long"
3879. #
3880. # Idents(seurat_sct) <- "Cluster.Type"
3881. #
3882. # Idents(seurat_sct) <- "Type.Cluster"
3884. #### check which ident is active by dimplot ####
3886. unique(sort([email hidden]))
3888. view([email hidden])
3890. DimPlot(seurat_sct,
3891. reduction = "umap",
3892. label = TRUE,
3893. label.size = 5,
3894. pt.size = 1,
3895. raster = F) +
3896. NoLegend()
3898. DimPlot(seurat_sct,
3899. reduction = "umap",
3900. label = TRUE,
3901. label.size = 5,
3902. pt.size = 1,
3903. cells.highlight = WhichCells(seurat_sct, idents = "Rods")) +
3904. NoLegend()
3906. DimPlot(seurat_sct,
3907. reduction = "umap",
3908. label = TRUE,
3909. label.size = 5,
3910. pt.size = 1,
3911. cells.highlight = WhichCells(seurat_sct, idents = "Cones")) +
3912. NoLegend()
3914. DimPlot(seurat_sct,
3915. reduction = "umap",
3916. label = TRUE,
3917. label.size = 5,
3918. pt.size = 1,
3919. cols = palopal) +
3920. NoLegend()
3922. #### save RDS containing RNA normalized data AND CELLTYPE LABELS ####
3924. # saveRDS(seurat_sct, paste0(prefixPCres,"_seurat_after_RNAnorm_CELLTYPES.rds"))
3926. seurat_sct <- readRDS(file = paste0(prefixPCres,"_seurat_after_RNAnorm_CELLTYPES.rds"))
3929. seurat_sct <- readRDS("~/Desktop/DEG_and_Pathway_Analysis/20230922_rat_retina_with_GEO/20230922_Rat_retina_with_GEO_int_12PCs_res.0.5_seurat_after_RNAnorm_CELLTYPES.rds")
3931. View([email hidden])
3933. DimPlot(seurat_sct, reduction = 'umap', group.by = "Cell.Type", split.by = "treatment",label = T, pt.size = 0.5, raster=FALSE)
3934. #### save seurat object as an hdf5 file for import into python ####
3936. SaveH5Seurat(seurat_sct,
3937. filename = paste(prefixPCres,"H5seurat_obj.H5seurat",sep = "_"),
3938. overwrite = FALSE,
3939. verbose = TRUE)
3941. #### Find Cluster Markers for all cells as celltypes ####
3943. # find all markers
3944. cluster_markers <- FindAllMarkers(seurat_sct,
3945. only.pos = TRUE,
3946. min.pct = 0.50,
3947. logfc.threshold = 0.4)
3949. write.csv(cluster_markers,
3950. file = paste(prefixPCres,"clustermarkers_CELLTYPES_pos_only_minpct0.50_logfcthresh0.4.csv", sep = "_"))
3952. #### find DEGs between treatments regardless of cell type ####
3954. Idents(object = seurat_sct) <- "treatment"
3956. cluster_markers <- FindAllMarkers(seurat_sct,
3957. only.pos = TRUE,
3958. min.pct = 0.25,
3959. logfc.threshold = 0.4)
3961. write.csv(cluster_markers, file =
3962. paste(prefixPCres,"degs_treatment_all_cells_minpct0.25_logfcthresh0.4.csv", sep = "_"))
3964. Idents(object = seurat_sct) <- ident_res
3966. #### find DEGs between time points regardless of cell type ####
3968. Idents(object = seurat_sct) <- "time_point"
3970. cluster_markers <- FindAllMarkers(seurat_sct,
3971. only.pos = TRUE,
3972. min.pct = 0.25,
3973. logfc.threshold = 0.4)
3975. write.csv(cluster_markers, file =
3976. paste(prefixPCres,"degs_timepoint_all_cells_minpct0.25_logfcthresh0.4.csv", sep = "_"))
3978. Idents(object = seurat_sct) <- ident_res
3980. #### find DEGs between time points and treatment regardless of cell type ####
3982. Idents(object = seurat_sct) <- "time_point_treatment"
3984. cluster_markers <- FindAllMarkers(seurat_sct,
3985. only.pos = TRUE,
3986. min.pct = 0.25,
3987. logfc.threshold = 0.4)
3989. write.csv(cluster_markers, file =
3990. paste(prefixPCres,"degs_timepoint_treatment_all_cells_minpct0.25_logfcthresh0.4.csv", sep = "_"))
3992. Idents(object = seurat_sct) <- ident_res
3994. Idents(object = seurat_sct) <- "Cell.Type"
3995. #### for loop for DEGs within clusters TR vs UT ####
3997. clusterlist <- c(levels([email hidden]$Cell.Type))
3999. clusterlist
4000. length(clusterlist)
4002. for (k in 1:length(clusterlist)){
4003. message(crayon::red(paste("Calculating DEGs within",
4004. clusterlist[k],
4005. sep = " ")))
4007. nam <- paste("clusterDEG_TR_vs_UT", clusterlist[k], sep = "_")
4009. assign(nam, FindMarkers(seurat_sct,
4010. ident.1 = "Cell_Treated",
4011. group.by = "treatment",
4012. subset.ident = clusterlist[k]))
4014. write.csv(get(nam), file = paste(prefixPCres,
4015. nam,
4016. "list.csv",
4017. sep = "_"),
4018. row.names = TRUE)
4020. }
4022. # positive LFC means upregulation in Cell_Treated cells
4024. #### for loop for DEGs within clusters by time_point ####
4026. # get the name of clusters
4027. clusterlist <- c(levels([email hidden]$Cell.Type))
4028. clusterlist
4030. # set the groups to compare
4031. grouplist <- c(levels([email hidden]$time_point))
4032. grouplist <- c(grouplist, grouplist[1])
4033. grouplist
4035. for (i in 1:(length(grouplist)-1)){
4036. for (k in 1:length(clusterlist)){
4037. message(crayon::red(paste("Calculating DEGs between",
4038. grouplist[i],
4039. "and",
4040. grouplist[i+1],
4041. "within the",
4042. clusterlist[k],
4043. "cluster",
4044. sep = " ")))
4046. nam <- paste("clusterDEG",
4047. grouplist[i],
4048. "vs",
4049. grouplist[i+1],
4050. "in",
4051. clusterlist[k],
4052. sep = "_")
4054. assign(nam, FindMarkers(seurat_sct,
4055. ident.1 = grouplist[i],
4056. ident.2 = grouplist[i+1],
4057. group.by = "time_point",
4058. subset.ident = clusterlist[k]))
4060. write.csv(get(nam), file = paste(prefixPCres,
4061. nam,
4062. "list.csv",
4063. sep = "_"),
4064. row.names = TRUE)
4066. }
4067. }
4071. #### for loop for DEGs within clusters by time_point ####
4073. # get the name of clusters
4074. clusterlist <- c(levels([email hidden]$Cell.Type))
4075. clusterlist
4077. # set the groups to compare
4078. grouplist <- c(levels([email hidden]$time_point_treatment))
4079. grouplist
4081. # Loop through each pair of groups (i vs j)
4082. for (i in 1:(length(grouplist)-1)){
4083. for (j in (i+1):length(grouplist)){
4084. group1 <- grouplist[i]
4085. group2 <- grouplist[j]
4087. # Run find markers on group1 and group2 for each element of clusterlist
4088. # print the comparison to the terminal
4089. for (k in 1:length(clusterlist)){
4090. message(crayon::red(paste("Calculating DEGs between",
4091. group1,
4092. "and",
4093. group2,
4094. "within the",
4095. clusterlist[k],
4096. "cluster",
4097. sep = " ")))
4099. # create the results variable unique to each pair of groups and
4100. # element of clusterlist
4101. nam <- paste("clusterDEG",
4102. group1,
4103. "vs",
4104. group2,
4105. "in",
4106. clusterlist[k],
4107. sep = "_")
4109. # run FindMarkers and save results in the "nam" variable
4110. assign(nam, FindMarkers(seurat_sct,
4111. ident.1 = group1,
4112. ident.2 = group2,
4113. group.by = "time_point_treatment",
4114. subset.ident = clusterlist[k]))
4116. # write out the results
4117. write.csv(get(nam), file = paste(prefixPCres,
4118. nam,
4119. "list.csv",
4120. sep = "_"),
4121. row.names = TRUE)
4122. }
4123. }
4124. }
4128. #### subsets for further DEG analysis ####
4130. # subset seurat obj by time point for treated vs untreated DEGs
4131. Idents(seurat_sct) <- "time_point"
4133. p60seurat <- subset(seurat_sct, idents = "p60")
4134. rm(seurat_sct)
4135. p90seurat <- subset(seurat_sct, idents = "p90")
4136. Idents(p90seurat) <- "Cell.Type"
4138. Cone_p90seurat <- subset(p90seurat, idents = "CONE")
4139. Cone <- subset(Cone_p90seurat, subset = Arr3 > .5)
4141. CONE <- WhichCells(object=Cone_p90seurat, expression = Arr3 > .5)
4142. Idents(Cone) <- "treatment"
4143. Cone_DEG <- FindMarkers(Cone,
4144. logfc.threshold = 0.2,
4145. min.pct = 0.2,
4146. ident.1 = "Cell_Treated",
4147. ident.2 = "Untreated")
4148. setwd("~/Desktop/DEG_and_Pathway_Analysis/Rat_DEGs/RAT_Pathway_analysis_V5")
4149. write.csv(Cone_DEG,"Cone_p90seurat1.csv")
4151. Cone_DEG1 <- FindMarkers(Cone,
4152. ident.1 = "Cell_Treated",
4153. ident.2 = "Untreated")
4154. write.csv(Cone_DEG1,"Cone_p90seurat.csv")
4159. VlnPlot(Cone, features = "Tyrobp", pt.size=0)
4161. FeaturePlot(Cone, features = "Apoe", min.cutoff = "q9", split.by = "treatment", cols = c("grey", "blue"), pt.size = .8)
4163. DimPlot(p90seurat, reduction = 'umap', group.by = "Cell.Type", split.by = "treatment", label = T, pt.size = 0.5, raster=FALSE)
4164. DimPlot(Cone_p90seurat, reduction = 'umap', split.by = "treatment", label = T, pt.size = 0.5, raster=FALSE)
4165. View([email hidden])
4166. View([email hidden])
4167. Idents(p60seurat) <- "treatment"
4168. cluster_markers <- FindAllMarkers(p60seurat, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.4)
4170. write.csv(cluster_markers, file =
4171. paste(prefixPCres,"degs_treatment_p60_minpct0.25_logfcthresh0.4.csv", sep = "_"))
4173. Idents(p90seurat) <- "treatment"
4174. cluster_markers <- FindAllMarkers(p90seurat, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.4)
4176. write.csv(cluster_markers, file =
4177. paste(prefixPCres,"degs_treatment_p90_minpct0.25_logfcthresh0.4.csv", sep = "_"))
4179. # subset seurat obj by treatment for p60 vs p90 DEGs
4180. Idents(seurat_sct) <- "treatment"
4182. TRseurat <- subset(seurat_sct, idents = "Cell_Treated")
4183. [email hidden] <- droplevels([email hidden])
4185. UTseurat <- subset(seurat_sct, idents = "Untreated")
4186. [email hidden] <- droplevels([email hidden])
4188. Idents(UTseurat) <- "time_point_treatment"
4190. UTseuratp60p90 <- subset(UTseurat,
4191. idents = "WT_Untreated",
4192. invert = T)
4193. [email hidden] <- droplevels([email hidden])
4195. Idents(object = seurat_sct) <- "Cell.Type"
4196. # Idents(object = p60seurat) <- "Cell.Type"
4197. # Idents(object = p90seurat) <- "Cell.Type"
4198. Idents(object = TRseurat) <- "Cell.Type"
4199. Idents(object = UTseurat) <- "Cell.Type"
4200. Idents(object = UTseuratp60p90) <- "Cell.Type"
4202. DimPlot(seurat_sct,
4203. reduction = "umap",
4204. label = TRUE,
4205. label.size = 5,
4206. pt.size = 1,
4207. raster = F) +
4208. NoLegend()
4210. DimPlot(TRseurat,
4211. reduction = "umap",
4212. label = TRUE,
4213. label.size = 5,
4214. pt.size = 1,
4215. raster = F) +
4216. NoLegend()
4218. DimPlot(UTseurat,
4219. reduction = "umap",
4220. label = TRUE,
4221. label.size = 5,
4222. pt.size = 1,
4223. raster = F) +
4224. NoLegend()
4226. DimPlot(UTseuratp60p90,
4227. reduction = "umap",
4228. label = TRUE,
4229. label.size = 5,
4230. pt.size = 1,
4231. raster = F) +
4232. NoLegend()
4234. #### DEGs in RODS Treated p60 vs p90 ####
4236. rod_tr_tp_markers <- FindMarkers(TRseurat,
4237. logfc.threshold = 0.4,
4238. min.pct = 0.25,
4239. only.pos = F,
4240. ident.1 = "p60",
4241. group.by = "time_point",
4242. subset.ident = "ROD")
4244. write.csv(rod_tr_tp_markers, file =
4245. paste(prefixPCres,"degs_ROD_p60_TR_vs_p90_TR",
4246. "minpct0.25_logfcthresh0.4.csv", sep = "_"))
4248. #### DEGs in RODS Untreated p60 vs p90 ####
4250. rod_ut_tp_markers <- FindMarkers(UTseuratp60p90,
4251. logfc.threshold = 0.4,
4252. min.pct = 0.25,
4253. only.pos = F,
4254. ident.1 = "p60",
4255. group.by = "time_point",
4256. subset.ident = "ROD")
4258. write.csv(rod_tr_tp_markers, file =
4259. paste(prefixPCres,"degs_ROD_p60_UT_vs_p90_UT",
4260. "minpct0.25_logfcthresh0.4.csv", sep = "_"))
4262. #### DEGs in CONES Treated p60 vs p90 ####
4264. cone_tr_tp_markers <- FindMarkers(TRseurat,
4265. logfc.threshold = 0.4,
4266. min.pct = 0.25,
4267. only.pos = F,
4268. ident.1 = "p60",
4269. group.by = "time_point",
4270. subset.ident = "CONE")
4272. write.csv(cone_tr_tp_markers, file =
4273. paste(prefixPCres,"degs_CONE_p60_TR_vs_p90_TR",
4274. "minpct0.25_logfcthresh0.4.csv", sep = "_"))
4276. #### DEGs in CONES Untreated p60 vs p90 ####
4278. cone_ut_tp_markers <- FindMarkers(UTseuratp60p90,
4279. logfc.threshold = 0.4,
4280. min.pct = 0.25,
4281. only.pos = F,
4282. ident.1 = "p60",
4283. group.by = "time_point",
4284. subset.ident = "CONE")
4286. write.csv(cone_tr_tp_markers, file =
4287. paste(prefixPCres,"degs_CONE_p60_UT_vs_p90_UT",
4288. "minpct0.25_logfcthresh0.4.csv", sep = "_"))
4290. #### DEGs in MULG Treated p60 vs p90 ####
4292. mulg_tr_tp_markers <- FindMarkers(TRseurat,
4293. logfc.threshold = 0.4,
4294. min.pct = 0.25,
4295. only.pos = F,
4296. ident.1 = "p60",
4297. group.by = "time_point",
4298. subset.ident = "MULG")
4300. write.csv(mulg_tr_tp_markers, file =
4301. paste(prefixPCres,"degs_MULG_p60_TR_vs_p90_TR",
4302. "minpct0.25_logfcthresh0.4.csv", sep = "_"))
4304. #### DEGs in MULG Untreated p60 vs p90 ####
4306. mulg_ut_tp_markers <- FindMarkers(UTseuratp60p90,
4307. logfc.threshold = 0.4,
4308. min.pct = 0.25,
4309. only.pos = F,
4310. ident.1 = "p60",
4311. group.by = "time_point",
4312. subset.ident = "MULG")
4314. write.csv(mulg_tr_tp_markers, file =
4315. paste(prefixPCres,"degs_MULG_p60_UT_vs_p90_UT",
4316. "minpct0.25_logfcthresh0.4.csv", sep = "_"))
4318. #### save subset objects ####
4319. saveRDS(p60seurat, paste0(prefixPCres,"_p60_seurat.rds"))
4321. # p60seurat <- readRDS(file = paste0(prefixPCres,"_p60_seurat.rds"))
4323. saveRDS(p90seurat, paste0(prefixPCres,"_p90_seurat.rds"))
4325. p90seurat <- readRDS(file = paste0(prefixPCres,"_p90_seurat.rds"))
4327. saveRDS(TRseurat, paste0(prefixPCres,"_cell_treated_seurat.rds"))
4329. # TRseurat <- readRDS(file = paste0(prefixPCres,"_cell_treated_seurat.rds"))
4331. saveRDS(UTseurat, paste0(prefixPCres,"_untreated_seurat.rds"))
4333. # UTseurat <- readRDS(file = paste0(prefixPCres,"_untreated_seurat.rds"))
4335. #### p60 TR vs UT DEGs within cell types ####
4336. Idents(p60seurat) <- "Cell.Type"
4338. clusterlist <- c(levels([email hidden]$Cell.Type))
4340. clusterlist
4342. for (k in 1:length(clusterlist)){
4344. nam <- paste("p60",
4345. "clusterDEG_treatment",
4346. clusterlist[k], sep = "_")
4347. assign(nam, FindMarkers(p60seurat,
4348. ident.1 = "Cell_Treated",
4349. group.by = "treatment",
4350. subset.ident = clusterlist[k]))
4352. write.csv(get(nam),
4353. file = paste(prefixPCres,
4354. nam,
4355. "list.csv",
4356. sep = "_"),
4357. row.names = TRUE)
4359. }
4361. # positive LFC means upregulation in Cell_Treated cells
4363. #### p90 TR vs UT DEGs within cell types ####
4364. Idents(p90seurat) <- "Cell.Type"
4366. clusterlist <- c(levels([email hidden]$Cell.Type))
4368. clusterlist
4370. for (k in 1:length(clusterlist)){
4372. nam <- paste("p90",
4373. "clusterDEG_treatment",
4374. clusterlist[k], sep = "_")
4375. assign(nam, FindMarkers(p90seurat,
4376. ident.1 = "Cell_Treated",
4377. group.by = "treatment",
4378. subset.ident = clusterlist[k]))
4380. write.csv(get(nam),
4381. file = paste(prefixPCres,
4382. nam,
4383. "list.csv",
4384. sep = "_"),
4385. row.names = TRUE)
4387. }
4389. # positive LFC means upregulation in Cell_Treated cells
4391. #### ACTIONet ####
4392. library(ACTIONet)
4393. library(SingleCellExperiment)
4394. # convert seurat object to a single cell experiment
4395. sce <- Seurat::as.SingleCellExperiment(seurat_sct)
4397. # convert single cell experiment to an ACTIONet experiment
4398. ace <- as.ACTIONetExperiment(sce)
4399. ace = normalize.ace(ace)
4400. ace = reduce.ace(ace)
4401. ace = runACTIONet(ace)
4403. # Annotate cell-types
4404. data("curatedMarkers_human")
4405. # markers = c(curatedMarkers_human$Brain$PFC$Mohammadi2020$marker.genes,
4406. # curatedMarkers_human$Brain$PFC$Velmeshev2019$marker.genes,
4407. # curatedMarkers_human$Brain$PFC$Schirmer2019$marker.genes,
4408. # curatedMarkers_human$Brain$PFC$MathysDavila2019$marker.genes,
4409. # curatedMarkers_human$Brain$PFC$Wang2018$marker.genes,
4410. # curatedMarkers_human$Brain$PFC$Layers$marker.genes)
4412. markers <- c(curatedMarkers_human$Retina$marker.genes)
4414. # List of human gene names
4415. markers
4417. h1 <- markers$Rods
4418. h2 <- markers$Cones
4419. h3 <- markers$RGCs
4420. h4 <- markers$BPs
4421. h5 <- markers$ACs
4422. h6 <- markers$HCs
4423. h7 <- markers$Macroglia
4424. h8 <- markers$Microglia
4425. h9 <- markers$Endo
4427. #### convert human genes to rat genes ####
4428. # Load the biomaRt library
4429. library(biomaRt)
4431. # convert human to rat gene names
4433. # listEnsemblArchives()
4434. # https://oct2022.archive.ensembl.org DIDN"T WORK
4435. # https://jul2022.archive.ensembl.org DIDN"T WORK
4436. # https://apr2022.archive.ensembl.org DIDN"T WORK
4437. # https://dec2021.archive.ensembl.org THIS ONE WORKED
4438. # https://may2021.archive.ensembl.org
4440. human <- useMart("ensembl",
4441. dataset = "hsapiens_gene_ensembl",
4442. host = "https://dec2021.archive.ensembl.org")
4444. rat <- useMart("ensembl",
4445. dataset = "rnorvegicus_gene_ensembl",
4446. host = "https://dec2021.archive.ensembl.org")
4448. r1 = getLDS(attributes = c("hgnc_symbol"),
4449. filters = "hgnc_symbol",
4450. values = h1,
4451. mart = human,
4452. attributesL = c("rgd_symbol"),
4453. martL = rat,
4454. uniqueRows=T)
4456. r2 <- getLDS(attributes = c("hgnc_symbol"),
4457. filters = "hgnc_symbol",
4458. values = h2,
4459. mart = human,
4460. attributesL = c("rgd_symbol"),
4461. martL = rat,
4462. uniqueRows=T)
4464. r3 <- getLDS(attributes = c("hgnc_symbol"),
4465. filters = "hgnc_symbol",
4466. values = h3,
4467. mart = human,
4468. attributesL = c("rgd_symbol"),
4469. martL = rat,
4470. uniqueRows=T)
4472. r4 <- getLDS(attributes = c("hgnc_symbol"),
4473. filters = "hgnc_symbol",
4474. values = h4,
4475. mart = human,
4476. attributesL = c("rgd_symbol"),
4477. martL = rat,
4478. uniqueRows=T)
4480. r5 <- getLDS(attributes = c("hgnc_symbol"),
4481. filters = "hgnc_symbol",
4482. values = h5,
4483. mart = human,
4484. attributesL = c("rgd_symbol"),
4485. martL = rat,
4486. uniqueRows=T)
4488. r6 <- getLDS(attributes = c("hgnc_symbol"),
4489. filters = "hgnc_symbol",
4490. values = h6,
4491. mart = human,
4492. attributesL = c("rgd_symbol"),
4493. martL = rat,
4494. uniqueRows=T)
4496. r7 <- getLDS(attributes = c("hgnc_symbol"),
4497. filters = "hgnc_symbol",
4498. values = h7,
4499. mart = human,
4500. attributesL = c("rgd_symbol"),
4501. martL = rat,
4502. uniqueRows=T)
4504. r8 <- getLDS(attributes = c("hgnc_symbol"),
4505. filters = "hgnc_symbol",
4506. values = h8,
4507. mart = human,
4508. attributesL = c("rgd_symbol"),
4509. martL = rat,
4510. uniqueRows=T)
4512. r9 <- getLDS(attributes = c("hgnc_symbol"),
4513. filters = "hgnc_symbol",
4514. values = h9,
4515. mart = human,
4516. attributesL = c("rgd_symbol"),
4517. martL = rat,
4518. uniqueRows=T)
4520. #### make rat reference list ####
4521. rat_markers <- list()
4523. rat_markers$Rods <- r1$RGD.symbol
4524. rat_markers$Cones <- r2$RGD.symbol
4525. rat_markers$RGCs <- r3$RGD.symbol
4526. rat_markers$BPs <- r4$RGD.symbol
4527. rat_markers$ACs <- r5$RGD.symbol
4528. rat_markers$HCs <- r6$RGD.symbol
4529. rat_markers$Macroglia <- r7$RGD.symbol
4530. rat_markers$Microglia <- r8$RGD.symbol
4531. rat_markers$Endo <- r9$RGD.symbol
4533. rat_markers
4537. rat_markers$Rods <- c("Ahi1","Rho","Gnat1","Nrl","Nr2e3",
4538. "Pde6a","Pde6b","Cngb1","Pdc","Rp1",
4539. "Guca1a","Scaper","Ppef2")
4540. rat_markers$Cones <- c("Arr3","Opn1mw","Opn1sw","Gnat2","Pde6h")
4541. rat_markers$Mueller <- c("Glul","Rlbp1","Zfp36l1","Dbi","Apoe",
4542. "Slc1a3","Sparc","Gfap","Aqp4","Gpr37",
4543. "S100b")
4544. rat_markers$BPs <- c("Scgn","Vsx1","Vsx2","Pcp2","Isl1",
4545. "Grm6","Trnp1","Tmem215","Prkca","Camk2b")
4546. rat_markers$ACs <- c("Calb1","Gad1","Gad2","C1ql1","C1ql2",
4547. "Snhg11","Tfap2b","Tfap2c","Pcsk1n","Slc6a9")
4548. rat_markers$HCs <- c("Onecut1","Onecut2","Lhx1","Slc4a3","Calb1",
4549. "Septin4","Tpm3")
4550. rat_markers$Macroglia <- c("Apoe","Slc1a3","Trpm3","Glul","Clu" )
4551. rat_markers$Microglia <- c("Aif1","Tyrobp","Cx3cr1","Tmem119","Ctsd",
4552. "Ccl4","C1qa","C1qb","C1qc","Cd163")
4553. rat_markers$RGCs <- c("Rbpms","Pou4f1","Pou4f2","Slc17a6","Nefm",
4554. "Nefl","Sncg" )
4555. # rat_markers$Endo <- c("Tm4sf1","Apold1","Adamts9","Igfbp7","Cd34",
4556. # "Rgs5","Cldn5","Cdh5")
4557. # rat_markers$RPEpi <- c("Mitf","Tjp1","Rpe65","Rlbp1","Best1")
4559. # rat_markers$Vascular <- c("Trpm1","Igfbp7")
4560. # rat_markers$Pericytes <- c("Pecam1","Acta2")
4561. rat_markers$Ret <- c("Hbb-bs")
4563. rat_markers
4565. #### annotate cells and export results to original seurat object ####
4566. annot.out <- annotate.cells.using.markers(ace, rat_markers)
4568. ace$ACEcelltypes <- annot.out$Label
4570. # Export results to the seurat object as a new metadata column
4572. seurat_sct[["ACEcelltypes"]] <- ace$ACEcelltypes
4574. # rm(h1,h2,h3,h4,h5,h6,h7,h8,h9,
4575. # r1,r2,r3,r4,r5,r6,r7,r8,r9,
4576. # markers,rat_markers,human,rat,sce,
4577. # curatedMarkers_human,annot.out,ace)
4579. #### switch to ACEcelltypes ####
4581. Idents(seurat_sct) <- "ACEcelltypes"
4583. DimPlot(seurat_sct,
4584. reduction = "umap",
4585. label = TRUE,
4586. label.size = 7,
4587. pt.size = 1)
4589. #### Module Scores for custom gene sets ####
4591. irm_1125 <- read.csv("/Users/bells/Library/CloudStorage/Box-Box/20201124_MG_scRNAseq(old)/gene_module_lists/IRM_11_25.csv", header = FALSE)
4592. irm_1125$V1 <- as.character(irm_1125$V1)
4593. irm_1125 <- list(irm_1125$V1)
4594. seurat_sct <- AddModuleScore(object = seurat_sct, features = irm_1125, search = TRUE, name = "IRM_11_25")
4596. arm_pos_1125 <- read.csv("/Users/bells/Library/CloudStorage/Box-Box/20201124_MG_scRNAseq(old)/gene_module_lists/ARM_POS_11_25.csv", header = FALSE)
4597. arm_pos_1125$V1 <- as.character(arm_pos_1125$V1)
4598. arm_pos_1125 <- list(arm_pos_1125$V1)
4599. seurat_sct <- AddModuleScore(object = seurat_sct, features = arm_pos_1125, search = TRUE, name = "ARM_POS_11_25")
4601. arm_neg_1125 <- read.csv("/Users/bells/Library/CloudStorage/Box-Box/20201124_MG_scRNAseq(old)/gene_module_lists/ARM_NEG_11_25.csv", header = FALSE)
4602. arm_neg_1125$V1 <- as.character(arm_neg_1125$V1)
4603. arm_neg_1125 <- list(arm_neg_1125$V1)
4604. seurat_sct <- AddModuleScore(object = seurat_sct, features = arm_neg_1125, search = TRUE, name = "ARM_NEG_11_25")
4606. default_FP(seurat_sct,"IRM_11_251")
4607. default_FP(seurat_sct,"ARM_POS_11_251")
4608. default_FP(seurat_sct,"ARM_NEG_11_251")
4611. #### expression values for genes of interest ####
4613. expression_data <- FetchData(object = seurat_sct,
4614. vars = c("genotype",ident_res,"Cell.Type",
4615. "ISL1","NEFL","NEFM",
4616. "NEFH","NGRN","PRPH","ELAVL3",
4617. "PFKP","RCAN1","SELPLG","STMN2","TARDBP"))
4619. write.csv(expression_data,
4620. file = paste(prefixPCres,"expression","in_all_cells","by_genotype.csv", sep = "_"))
4622. #### VENN compare DEG gene lists ####
4623. setwd(datadir)
4625. # load data and convert the output "list" to a tibble
4626. # with the rownames converted to a column with the name "symbol"
4627. p60_tr_vs_ut <- read.csv(file = "./20230327_Rat_retina_int_12PCs_res.0.2_degs_treatment_p60_minpct0.25_logfcthresh0.4.csv",
4628. header = TRUE)
4629. colnames(p60_tr_vs_ut) <- c("symbol","p_val","ave_log2FC","pct.1","pct.2","p_val_adj","cluster","gene")
4631. p90_tr_vs_ut <- read.csv(file = "./20230327_Rat_retina_int_12PCs_res.0.2_degs_treatment_p90_minpct0.25_logfcthresh0.4.csv",
4632. header = TRUE)
4633. colnames(p90_tr_vs_ut) <- c("symbol","p_val","ave_log2FC","pct.1","pct.2","p_val_adj","cluster","gene")
4635. p60_tr_vs_ut <- as_tibble(p60_tr_vs_ut)
4636. p90_tr_vs_ut <- as_tibble(p90_tr_vs_ut)
4638. # split tibbles by "cluster" and p_val_adj < 0.05
4639. p60_up_in_tr <- filter(p60_tr_vs_ut, cluster == "Cell_Treated" & p_val_adj < 0.05)
4640. p60_up_in_ut <- filter(p60_tr_vs_ut, cluster == "Untreated" & p_val_adj < 0.05)
4641. p90_up_in_tr <- filter(p90_tr_vs_ut, cluster == "Cell_Treated" & p_val_adj < 0.05)
4642. p90_up_in_ut <- filter(p90_tr_vs_ut, cluster == "Untreated" & p_val_adj < 0.05)
4644. # venn diagrams of overlapping genes
4646. ### Up in Treated condition ###
4647. treated_up <- list(p60_up_in_tr$gene,p90_up_in_tr$gene)
4648. names(treated_up) <- c("p60_Up_in_Cell_Treated","p90_Up_in_Cell_Treated")
4650. venn_treated_up <- Venn(treated_up)
4652. plot(venn_treated_up, doWeights = TRUE, type = "circles")
4654. v1 <- compute.Venn(venn_treated_up, doWeights = TRUE, type = "circles")
4656. # change default labels
4657. setlabels <- VennGetSetLabels(v1)
4658. setlabels
4659. facelabels <- VennGetFaceLabels(v1)
4660. facelabels
4661. universe <- VennGetUniverseRange(v1)
4662. universe
4663. # setlabels[setlabels$Label=="p60_Up_in_Cell_Treated","hjust"] <- "right"
4664. # setlabels[setlabels$Label=="p90_Up_in_Cell_Treated","hjust"] <- "left"
4665. # setlabels[setlabels$Label=="YoungV_vs_iMACS_13","x"] <- -1.5
4666. # setlabels
4667. # v1 <- VennSetSetLabels(v1,setlabels)
4669. grid::grid.newpage()
4670. plot(v1)
4672. pdf(paste(prefixPCres,"venn","p60","vs","p90","up","in","Cell_Treated","overlapping","DEGs.pdf", sep = "_"))
4673. print(plot(v1))
4674. dev.off()
4676. # get lists of non-overlapping and overlapping genes
4677. p60_only_tr_up <- v1@IntersectionSets$'10'
4678. p60_only_tr_up
4679. overlap_tr_up <- v1@IntersectionSets$'11'
4680. overlap_tr_up
4681. p90_only_tr_up <- v1@IntersectionSets$'01'
4682. p90_only_tr_up
4684. write.csv(p60_only_tr_up, file = paste(prefixPCres,"venn","genes","up","in","p60","cell_treated.csv",sep = "_"))
4685. write.csv(overlap_tr_up, file = paste(prefixPCres,"venn","genes","up","in","p60_and_p90","cell_treated.csv",sep = "_"))
4686. write.csv(p90_only_tr_up, file = paste(prefixPCres,"venn","genes","up","in","p90","cell_treated.csv",sep = "_"))
4688. ### Up in Untreated condition ###
4689. untreated_up <- list(p60_up_in_ut$gene,p90_up_in_ut$gene)
4690. names(untreated_up) <- c("p60_Up_in_Untreated","p90_Up_in_Untreated")
4692. venn_untreated_up <- Venn(untreated_up)
4694. plot(venn_untreated_up, doWeights = TRUE, type = "circles")
4696. v2 <- compute.Venn(venn_untreated_up, doWeights = TRUE, type = "circles")
4698. # change default labels
4699. # setlabels <- VennGetSetLabels(v2)
4700. # setlabels
4701. # facelabels <- VennGetFaceLabels(v2)
4702. # facelabels
4703. # universe <- VennGetUniverseRange(v2)
4704. # universe
4705. # setlabels[setlabels$Label=="p60_Up_in_Cell_Treated","hjust"] <- "right"
4706. # setlabels[setlabels$Label=="p90_Up_in_Cell_Treated","hjust"] <- "left"
4707. # setlabels[setlabels$Label=="YoungV_vs_iMACS_13","x"] <- -1.5
4708. # setlabels
4709. # v1 <- VennSetSetLabels(v1,setlabels)
4711. grid::grid.newpage()
4712. plot(v2)
4714. pdf(paste(prefixPCres,"venn","p60","vs","p90","up","in","Untreated","overlapping","DEGs.pdf", sep = "_"))
4715. print(plot(v2))
4716. dev.off()
4718. # get lists of non-overlapping and overlapping genes
4719. p60_only_ut_up <- v2@IntersectionSets$'10'
4720. p60_only_ut_up
4721. overlap_ut_up <- v2@IntersectionSets$'11'
4722. overlap_ut_up
4723. p90_only_ut_up <- v2@IntersectionSets$'01'
4724. p90_only_ut_up
4726. write.csv(p60_only_ut_up, file = paste(prefixPCres,"venn","genes","up","in","p60","untreated.csv",sep = "_"))
4727. write.csv(overlap_ut_up, file = paste(prefixPCres,"venn","genes","up","in","p60_and_p90","untreated.csv",sep = "_"))
4728. write.csv(p90_only_ut_up, file = paste(prefixPCres,"venn","genes","up","in","p90","untreated.csv",sep = "_"))
4729. #### Venn in ROD cluster ####
4731. celltype <- "RODs"
4733. # function to export DEG lists to csv
4734. exportListToCSV <- function(data, filename) {
4735. # Determine the maximum length among all columns
4736. max_length <- max(sapply(data, length))
4738. # Pad shorter columns with empty strings
4739. filled_data <- lapply(data, function(x) {
4740. c(x, rep("", max_length - length(x)))
4741. })
4743. # Convert the list to a data frame
4744. df <- as.data.frame(filled_data)
4746. # Set the column names
4747. colnames(df) <- names(data)
4749. # Write the data frame to a CSV file
4750. write.csv(df, filename, row.names = FALSE)
4751. }
4753. # load the data
4754. treatment_DEGs <- read.csv(file = "DEGs_p60_p90_by_treatment_RODs.csv",
4755. header = TRUE)
4757. p60_TR_UP <- treatment_DEGs$p60_TR
4758. p60_UT_UP <- treatment_DEGs$p60_UT
4759. p90_TR_UP <- treatment_DEGs$p90_TR
4760. p90_UT_UP <- treatment_DEGs$p90_UT
4762. rm(treatment_DEGs)
4764. #### p60-p90 overlap in Treated condition ####
4765. treated_up <- list(p60_TR_UP,p90_TR_UP)
4766. names(treated_up) <- c("p60_Up_in_Cell_Treated","p90_Up_in_Cell_Treated")
4768. venn_treated_up <- Venn(treated_up)
4770. plot(venn_treated_up, doWeights = TRUE, type = "circles")
4772. v1 <- compute.Venn(venn_treated_up, doWeights = TRUE, type = "circles")
4774. grid::grid.newpage()
4775. plot(v1)
4777. pdf(paste(prefixPCres,"venn","p60","p90","up","in","Cell_Treated","overlapping","DEGs",celltype,".pdf", sep = "_"))
4778. print(plot(v1))
4779. dev.off()
4781. # get lists of non-overlapping and overlapping genes
4782. p60_only_tr_up <- v1@IntersectionSets$'10'
4783. p60_only_tr_up
4784. overlap_tr_up <- v1@IntersectionSets$'11'
4785. overlap_tr_up
4786. p90_only_tr_up <- v1@IntersectionSets$'01'
4787. p90_only_tr_up
4789. p60p90_TR_up <- list(p60_only_tr_up = p60_only_tr_up,
4790. overlap_tr_up = overlap_tr_up,
4791. p90_only_tr_up = p90_only_tr_up)
4792. p60p90_TR_up
4794. exportListToCSV(p60p90_TR_up, paste(prefixPCres,"venn",
4795. "p60","p90","cell_treated",
4796. "overlapping","DEGs",celltype,
4797. ".csv",sep = "_"))
4799. #### p60 TR p90 UT overlap ####
4800. treated_up <- list(p60_TR_UP,p90_UT_UP)
4801. names(treated_up) <- c("p60_Up_in_Cell_Treated","p90_Up_in_Unreated")
4803. venn_treated_up <- Venn(treated_up)
4805. plot(venn_treated_up, doWeights = TRUE, type = "circles")
4807. v1 <- compute.Venn(venn_treated_up, doWeights = TRUE, type = "circles")
4809. grid::grid.newpage()
4810. plot(v1)
4812. pdf(paste(prefixPCres,"venn","p60","TR","p90","UT",
4813. "overlapping","DEGs",celltype,".pdf", sep = "_"))
4814. print(plot(v1))
4815. dev.off()
4817. # get lists of non-overlapping and overlapping genes
4818. p60_only_tr_up <- v1@IntersectionSets$'10'
4819. p60_only_tr_up
4820. overlap_tr_up <- v1@IntersectionSets$'11'
4821. overlap_tr_up
4822. p90_only_ut_up <- v1@IntersectionSets$'01'
4823. p90_only_ut_up
4825. p60_TR_p90_UT <- list(p60_only_tr_up = p60_only_tr_up,
4826. overlap_tr_up = overlap_tr_up,
4827. p90_only_ut_up = p90_only_ut_up)
4828. p60_TR_p90_UT
4830. exportListToCSV(p60_TR_p90_UT, paste(prefixPCres,"venn",
4831. "p60","TR","p90","UT",
4832. "overlapping","DEGs",celltype,
4833. ".csv",sep = "_"))
4835. #### p60 UT p90 TR overlap ####
4836. treated_up <- list(p60_UT_UP,p90_TR_UP)
4837. names(treated_up) <- c("p60_Up_in_Unreated","p90_Up_in_Cell_Treated")
4839. venn_treated_up <- Venn(treated_up)
4841. plot(venn_treated_up, doWeights = TRUE, type = "circles")
4843. v1 <- compute.Venn(venn_treated_up, doWeights = TRUE, type = "circles")
4845. grid::grid.newpage()
4846. plot(v1)
4848. pdf(paste(prefixPCres,"venn","p60","UT","p90","TR",
4849. "overlapping","DEGs",celltype,".pdf", sep = "_"))
4850. print(plot(v1))
4851. dev.off()
4853. # get lists of non-overlapping and overlapping genes
4854. p60_only_ut_up <- v1@IntersectionSets$'10'
4855. p60_only_ut_up
4856. overlap_tr_up <- v1@IntersectionSets$'11'
4857. overlap_tr_up
4858. p90_only_tr_up <- v1@IntersectionSets$'01'
4859. p90_only_tr_up
4861. p60_UT_p90_TR <- list(p60_only_ut_up = p60_only_ut_up,
4862. overlap_tr_up = overlap_tr_up,
4863. p90_only_tr_up = p90_only_tr_up)
4864. p60_UT_p90_TR
4866. exportListToCSV(p60_UT_p90_TR, paste(prefixPCres,"venn",
4867. "p60","UT","p90","TR",
4868. "overlapping","DEGs",celltype,
4869. ".csv",sep = "_"))
4871. #### p60-p90 overlap in Untreated condition ####
4872. untreated_up <- list(p60_UT_UP,p90_UT_UP)
4873. names(untreated_up) <- c("p60_Up_in_Untreated","p90_Up_in_Untreated")
4875. venn_untreated_up <- Venn(untreated_up)
4877. plot(venn_untreated_up, doWeights = TRUE, type = "circles")
4879. v1 <- compute.Venn(venn_untreated_up, doWeights = TRUE, type = "circles")
4881. grid::grid.newpage()
4882. plot(v1)
4884. pdf(paste(prefixPCres,"venn",
4885. "p60","p90",
4886. "up","in",
4887. "Untreated",
4888. "overlapping","DEGs",
4889. celltype,".pdf", sep = "_"))
4890. print(plot(v1))
4891. dev.off()
4893. # get lists of non-overlapping and overlapping genes
4894. p60_only_ut_up <- v1@IntersectionSets$'10'
4895. p60_only_ut_up
4896. overlap_ut_up <- v1@IntersectionSets$'11'
4897. overlap_ut_up
4898. p90_only_ut_up <- v1@IntersectionSets$'01'
4899. p90_only_ut_up
4901. p60p90_UT_up <- list(p60_only_ut_up = p60_only_ut_up,
4902. overlap_ut_up = overlap_ut_up,
4903. p90_only_ut_up = p90_only_ut_up)
4904. p60p90_UT_up
4906. exportListToCSV(p60p90_UT_up, paste(prefixPCres,"venn",
4907. "p60","p90","untreated",
4908. "overlapping","DEGs",celltype,
4909. ".csv",sep = "_"))
4911. #### Venn in MG cluster ####
4913. celltype <- "MG"
4915. # function to export DEG lists to csv
4916. exportListToCSV <- function(data, filename) {
4917. # Determine the maximum length among all columns
4918. max_length <- max(sapply(data, length))
4920. # Pad shorter columns with empty strings
4921. filled_data <- lapply(data, function(x) {
4922. c(x, rep("", max_length - length(x)))
4923. })
4925. # Convert the list to a data frame
4926. df <- as.data.frame(filled_data)
4928. # Set the column names
4929. colnames(df) <- names(data)
4931. # Write the data frame to a CSV file
4932. write.csv(df, filename, row.names = FALSE)
4933. }
4935. # load the data
4936. treatment_DEGs <- read.csv(file = "DEGs_p60_p90_by_treatment_MG.csv",
4937. header = TRUE)
4939. p60_TR_UP <- treatment_DEGs$p60_TR
4940. p60_UT_UP <- treatment_DEGs$p60_UT
4941. p90_TR_UP <- treatment_DEGs$p90_TR
4942. p90_UT_UP <- treatment_DEGs$p90_UT
4944. rm(treatment_DEGs)
4946. #### p60-p90 overlap in Treated condition ####
4947. treated_up <- list(p60_TR_UP,p90_TR_UP)
4948. names(treated_up) <- c("p60_Up_in_Cell_Treated","p90_Up_in_Cell_Treated")
4950. venn_treated_up <- Venn(treated_up)
4952. plot(venn_treated_up, doWeights = TRUE, type = "circles")
4954. v1 <- compute.Venn(venn_treated_up, doWeights = TRUE, type = "circles")
4956. grid::grid.newpage()
4957. plot(v1)
4959. pdf(paste(prefixPCres,"venn","p60","p90","up","in","Cell_Treated","overlapping","DEGs",celltype,".pdf", sep = "_"))
4960. print(plot(v1))
4961. dev.off()
4963. # get lists of non-overlapping and overlapping genes
4964. p60_only_tr_up <- v1@IntersectionSets$'10'
4965. p60_only_tr_up
4966. overlap_tr_up <- v1@IntersectionSets$'11'
4967. overlap_tr_up
4968. p90_only_tr_up <- v1@IntersectionSets$'01'
4969. p90_only_tr_up
4971. p60p90_TR_up <- list(p60_only_tr_up = p60_only_tr_up,
4972. overlap_tr_up = overlap_tr_up,
4973. p90_only_tr_up = p90_only_tr_up)
4974. p60p90_TR_up
4976. exportListToCSV(p60p90_TR_up, paste(prefixPCres,"venn",
4977. "p60","p90","cell_treated",
4978. "overlapping","DEGs",celltype,
4979. ".csv",sep = "_"))
4981. #### p60 TR p90 UT overlap ####
4982. treated_up <- list(p60_TR_UP,p90_UT_UP)
4983. names(treated_up) <- c("p60_Up_in_Cell_Treated","p90_Up_in_Untreated")
4985. venn_treated_up <- Venn(treated_up)
4987. plot(venn_treated_up, doWeights = TRUE, type = "circles")
4989. v1 <- compute.Venn(venn_treated_up, doWeights = TRUE, type = "circles")
4991. grid::grid.newpage()
4992. plot(v1)
4994. pdf(paste(prefixPCres,"venn","p60","TR","p90","UT",
4995. "overlapping","DEGs",celltype,".pdf", sep = "_"))
4996. print(plot(v1))
4997. dev.off()
4999. # get lists of non-overlapping and overlapping genes
5000. p60_only_tr_up <- v1@IntersectionSets$'10'
5001. p60_only_tr_up
5002. overlap_tr_up <- v1@IntersectionSets$'11'
5003. overlap_tr_up
5004. p90_only_ut_up <- v1@IntersectionSets$'01'
5005. p90_only_ut_up
5007. p60_TR_p90_UT <- list(p60_only_tr_up = p60_only_tr_up,
5008. overlap_tr_up = overlap_tr_up,
5009. p90_only_ut_up = p90_only_ut_up)
5010. p60_TR_p90_UT
5012. exportListToCSV(p60_TR_p90_UT, paste(prefixPCres,"venn",
5013. "p60","TR","p90","UT",
5014. "overlapping","DEGs",celltype,
5015. ".csv",sep = "_"))
5017. #### p60 UT p90 TR overlap ####
5018. treated_up <- list(p60_UT_UP,p90_TR_UP)
5019. names(treated_up) <- c("p60_Up_in_Untreated","p90_Up_in_Cell_Treated")
5021. venn_treated_up <- Venn(treated_up)
5023. plot(venn_treated_up, doWeights = TRUE, type = "circles")
5025. v1 <- compute.Venn(venn_treated_up, doWeights = TRUE, type = "circles")
5027. grid::grid.newpage()
5028. plot(v1)
5030. pdf(paste(prefixPCres,"venn","p60","UT","p90","TR",
5031. "overlapping","DEGs",celltype,".pdf", sep = "_"))
5032. print(plot(v1))
5033. dev.off()
5035. # get lists of non-overlapping and overlapping genes
5036. p60_only_ut_up <- v1@IntersectionSets$'10'
5037. p60_only_ut_up
5038. overlap_tr_up <- v1@IntersectionSets$'11'
5039. overlap_tr_up
5040. p90_only_tr_up <- v1@IntersectionSets$'01'
5041. p90_only_tr_up
5043. p60_UT_p90_TR <- list(p60_only_ut_up = p60_only_ut_up,
5044. overlap_tr_up = overlap_tr_up,
5045. p90_only_tr_up = p90_only_tr_up)
5046. p60_UT_p90_TR
5048. exportListToCSV(p60_UT_p90_TR, paste(prefixPCres,"venn",
5049. "p60","UT","p90","TR",
5050. "overlapping","DEGs",celltype,
5051. ".csv",sep = "_"))
5053. #### p60-p90 overlap in Untreated condition ####
5054. untreated_up <- list(p60_UT_UP,p90_UT_UP)
5055. names(untreated_up) <- c("p60_Up_in_Untreated","p90_Up_in_Untreated")
5057. venn_untreated_up <- Venn(untreated_up)
5059. plot(venn_untreated_up, doWeights = TRUE, type = "circles")
5061. v1 <- compute.Venn(venn_untreated_up, doWeights = TRUE, type = "circles")
5063. grid::grid.newpage()
5064. plot(v1)
5066. pdf(paste(prefixPCres,"venn",
5067. "p60","p90",
5068. "up","in",
5069. "Untreated",
5070. "overlapping","DEGs",
5071. celltype,".pdf", sep = "_"))
5072. print(plot(v1))
5073. dev.off()
5075. # get lists of non-overlapping and overlapping genes
5076. p60_only_ut_up <- v1@IntersectionSets$'10'
5077. p60_only_ut_up
5078. overlap_ut_up <- v1@IntersectionSets$'11'
5079. overlap_ut_up
5080. p90_only_ut_up <- v1@IntersectionSets$'01'
5081. p90_only_ut_up
5083. p60p90_UT_up <- list(p60_only_ut_up = p60_only_ut_up,
5084. overlap_ut_up = overlap_ut_up,
5085. p90_only_ut_up = p90_only_ut_up)
5086. p60p90_UT_up
5088. exportListToCSV(p60p90_UT_up, paste(prefixPCres,"venn",
5089. "p60","p90","untreated",
5090. "overlapping","DEGs",celltype,
5091. ".csv",sep = "_"))
5093. #### Trying ggVennDiagram and ggvenn instead of Vennerable ####
5095. # ggVennDiagram
5096. if (!require(devtools)) install.packages("devtools")
5097. devtools::install_github("gaospecial/ggVennDiagram")
5099. library("ggVennDiagram")
5101. ggVennDiagram(treated_up,label_alpha = 0) +
5102. ggplot2::scale_fill_gradient(low="blue",high = "yellow")
5104. # ggvenn
5105. if (!require(devtools)) install.packages("devtools")
5106. devtools::install_github("yanlinlin82/ggvenn")
5108. library("ggvenn")
5110. ggvenn(treated_up,
5111. # fill_color = c("green", "red"),
5112. stroke_size = 1,
5113. # set_name_size = 9,
5114. text_size = 8,
5115. show_percentage = FALSE,
5116. auto_scale = FALSE
5117. )
5119. ?ggvenn
5120. ?geom_venn
5123. #### Heatmap ####
5125. Idents(seurat_sct) <- "treatment"
5127. view([email hidden])
5129. p60_treatment_DEGs <- read_csv(file = "20230327_Rat_retina_int_12PCs_res.0.2_degs_treatment_p60_minpct0.25_logfcthresh0.4.csv")
5131. p60_treatment_DEGs %>% sort()
5133. p60_treatment_DEGs <- p60_treatment_DEGs[order(p60_treatment_DEGs$p_val_adj), ]
5135. DoHeatmap(seurat_sct, features = c(p60_treatment_DEGs$gene[1:30]),
5136. cells = WhichCells(seurat_sct,
5137. expression = time_point == "p60"))
5139. p90_treatment_DEGs <- read_csv(file = "20230327_Rat_retina_int_12PCs_res.0.2_degs_treatment_p90_minpct0.25_logfcthresh0.4.csv")
5141. DoHeatmap(seurat_sct, features = p90_treatment_DEGs$gene)
5143. #### volcano plots ####
5145. p60_treatment_DEGs <- read_csv(file = "20230327_Rat_retina_int_12PCs_res.0.2_degs_treatment_p60_minpct0.25_logfcthresh0.4.csv")
5147. TR_DEGS <- subset(p60_treatment_DEGs, subset = cluster == "Cell_Treated")
5148. UT_DEGS <- subset(p60_treatment_DEGs, subset = cluster == "Untreated")
5150. UT_DEGS$avg_log2FC <- -1 * UT_DEGS$avg_log2FC
5152. DEGS <- rbind(TR_DEGS,UT_DEGS)
5154. volc <- DEGS
5156. for (i in 1:length(volc$p_val_adj)) {
5157. if (volc$p_val_adj[i] == 0) {
5158. volc$p_val_adj[i] <- 1e-293
5159. }
5160. }
5162. EnhancedVolcano(volc,
5163. lab = volc$gene,
5164. x = 'avg_log2FC',
5165. y = 'p_val_adj',
5166. xlim = c(min(volc$avg_log2FC, na.rm = TRUE) - 0.1,
5167. max(volc$avg_log2FC, na.rm = TRUE) + 0.1),
5168. ylim = c(0, max(-log10(volc$p_val_adj), na.rm = TRUE) + 100),
5169. # title = "Aging iMPs vs Aging Veh",
5170. titleLabSize = 30,
5171. caption = NULL,
5172. subtitle = NULL,
5173. xlab = bquote(bold(~Log["2"] ~ "fold change")),
5174. ylab = bquote(bold(~"-" ~Log["10"] ~ p_val_adj)),
5175. axisLabSize = 25,
5176. # col = volc_colors,
5177. colAlpha = 1/1,
5178. legendLabels = c("NS",
5179. expression(Log[2] ~ FC),
5180. "p_val_adj",
5181. expression(p_val_adj ~ and ~ log[2] ~ FC)),
5182. legendPosition = 'bottom',
5183. legendLabSize = 25,
5184. legendIconSize = 7,
5185. drawConnectors = TRUE,
5186. min.segment.length = 1,
5187. labSize = 7,
5188. max.overlaps = Inf,
5189. pCutoff = 0.05,
5190. FCcutoff = 0.58,
5191. gridlines.major = FALSE,
5192. gridlines.minor = FALSE,
5193. arrowheads = FALSE,
5194. pointSize = 7
5195. )
5197. p90_treatment_DEGs <- read_csv(file = "20230327_Rat_retina_int_12PCs_res.0.2_degs_treatment_p90_minpct0.25_logfcthresh0.4.csv")
5199. TR_DEGS <- subset(p90_treatment_DEGs, subset = cluster == "Cell_Treated")
5200. UT_DEGS <- subset(p90_treatment_DEGs, subset = cluster == "Untreated")
5202. UT_DEGS$avg_log2FC <- -1 * UT_DEGS$avg_log2FC
5204. DEGS <- rbind(TR_DEGS,UT_DEGS)
5206. volc <- DEGS
5208. for (i in 1:length(volc$p_val_adj)) {
5209. if (volc$p_val_adj[i] == 0) {
5210. volc$p_val_adj[i] <- 1e-293
5211. }
5212. }
5214. EnhancedVolcano(volc,
5215. lab = volc$gene,
5216. x = 'avg_log2FC',
5217. y = 'p_val_adj',
5218. xlim = c(min(volc$avg_log2FC, na.rm = TRUE) - 0.1,
5219. max(volc$avg_log2FC, na.rm = TRUE) + 0.1),
5220. ylim = c(0, max(-log10(volc$p_val_adj), na.rm = TRUE) + 100),
5221. # title = "Aging iMPs vs Aging Veh",
5222. titleLabSize = 30,
5223. caption = NULL,
5224. subtitle = NULL,
5225. xlab = bquote(bold(~Log["2"] ~ "fold change")),
5226. ylab = bquote(bold(~"-" ~Log["10"] ~ p_val_adj)),
5227. axisLabSize = 25,
5228. # col = volc_colors,
5229. colAlpha = 1/1,
5230. legendLabels = c("NS",
5231. expression(Log[2] ~ FC),
5232. "p_val_adj",
5233. expression(p_val_adj ~ and ~ log[2] ~ FC)),
5234. legendPosition = 'bottom',
5235. legendLabSize = 25,
5236. legendIconSize = 7,
5237. drawConnectors = TRUE,
5238. min.segment.length = 1,
5239. labSize = 3,
5240. max.overlaps = Inf,
5241. pCutoff = 0.05,
5242. FCcutoff = 1,
5243. gridlines.major = FALSE,
5244. gridlines.minor = FALSE,
5245. arrowheads = FALSE,
5246. pointSize = 3
5247. )
5250. #### ####
5251. #### Rod subset recluster ####
5252. s_rods <- subset(seurat_sct,
5253. idents = c("ROD"))
5255. #### SCT subset ####
5256. DefaultAssay(s_rods) <- "RNA"
5258. s_rods <- SCTransform(s_rods,
5259. vars.to.regress = c("percent.mt","percent.ribo"),
5260. verbose = TRUE)
5263. #### update project variable to include subset ####
5264. project_ROD <- "Rat_retina_RODS"
5266. #### sub PCA ####
5267. # this performs PCA on the seurat object
5268. s_rods <- RunPCA(s_rods, npcs = 50, verbose = TRUE)
5270. # make PC coordinate object a data frame
5271. xx.coord <- as.data.frame(s_rods@reductions$[email hidden])
5273. # make PC feature loadings object a data frame
5274. xx.gload <- as.data.frame(s_rods@reductions$[email hidden])
5276. # calculate eigenvalues for arrays
5278. # generate squares of all sample coordinates
5279. sq.xx.coord <- as.data.frame(xx.coord^2)
5280. # create empty list for eigenvalues first
5281. eig <- c()
5282. # calculate the eigenvalue for each PC in sq.xx.coord by taking the sqrt of the sum of squares
5283. for(i in 1:ncol(sq.xx.coord))
5284. eig[i] = sqrt(sum(sq.xx.coord[,i]))
5285. # calculate the total variance by adding up all the eigenvalues
5286. sum.eig <- sum(eig)
5287. # calculate the expected contribution of all PCs if they all contribute equally to the total variance
5288. expected.contribution <- sum.eig/(length(xx.coord)-1)
5289. # return the number of principal components with an eigenvalue greater than expected by equal variance
5290. meaningful.PCs <- sum(eig > expected.contribution)
5292. # create empty list for eigenvalue percentage
5293. eig.percent <- c()
5294. # calculate the percentage of the total variance by each PC eigenvalue
5295. for(i in 1:length(eig))
5296. eig.percent[i] = 100*eig[i]/sum.eig
5297. # sum of all eig.percent should total to 100
5298. sum(eig.percent)
5299. # create empty list for scree values
5300. scree <- c()
5301. # calculate a running total of variance contribution
5302. for(i in 1:length(eig))
5303. if(i == 1) scree[i] = eig.percent[i] else scree[i] = scree[i-1] + eig.percent[i]
5305. # create data frame for eigenvalue summaries
5306. eigenvalues <- data.frame("PC" = colnames(xx.coord), "eig" = eig, "percent" = eig.percent, "scree" = scree)
5308. # write csv for eigenvalues
5309. # write.csv(eigenvalues, file = paste0("./",date,"_",project,"_PCA_eigenvalues.csv"), row.names = F)
5311. # plot scree values
5312. plot(eigenvalues$percent, ylim = c(0,100), type = "S", xlab = "PC", ylab = "Percent of variance",
5313. main = paste0(date,"_",project," scree plot all samples PCA"))
5314. points(eigenvalues$scree, ylim = c(0,100), type = "p", pch = 16)
5315. lines(eigenvalues$scree)
5316. # add red line to indicate cut-off
5317. cut.off <- 100/(length(eig)-1)
5318. abline(h = cut.off, col = "red")
5319. # add blue line to indicate which PCs are meaningful and kept
5320. abline(v = meaningful.PCs, col = "blue")
5321. text(meaningful.PCs, cut.off, label = paste("cutoff PC",meaningful.PCs),
5322. adj = c(-0.1, -0.5))
5324. dev.copy(pdf, paste0("./",date,"_",project_ROD,"_scree_plot.pdf"))
5325. dev.off()
5327. rm(eigenvalues,sq.xx.coord,xx.coord,xx.gload,cut.off,
5328. eig,eig.percent,expected.contribution,i,scree,sum.eig)
5330. # meaningful.PCs <- 9
5332. #### sub Run UMAP and look at UMAP plots ####
5334. s_rods <- RunUMAP(s_rods,
5335. reduction = "pca",
5336. dims = 1:meaningful.PCs,
5337. verbose = TRUE)
5339. # update prefixed variable
5340. ROD_prefixPC <- paste0("./",date,"_",project_ROD,"_",meaningful.PCs,"PCs")
5342. ## UMAP plot by sample name ("orig.ident")
5343. DimPlot(s_rods, reduction = "umap", label = FALSE,
5344. pt.size = .25, group.by = "orig.ident")
5345. hirestiff(paste0(ROD_prefixPC,"_UMAP_by_sample_hires.tiff"))
5346. lowrestiff(paste0(ROD_prefixPC,"_UMAP_by_sample_lowres.tiff"))
5348. DimPlot(s_rods, reduction = "umap", label = FALSE,
5349. pt.size = .25, group.by = "treatment")
5350. hirestiff(paste0(ROD_prefixPC,"_UMAP_by_treatment_hires.tiff"))
5351. lowrestiff(paste0(ROD_prefixPC,"_UMAP_by_treatment_lowres.tiff"))
5353. DimPlot(s_rods, reduction = "umap", label = FALSE,
5354. pt.size = .25, split.by = "treatment", group.by = "treatment")
5355. hirestiff(paste0(ROD_prefixPC,"_UMAP_split_by_treatment_hires.tiff"))
5356. lowrestiff(paste0(ROD_prefixPC,"_UMAP_split_by_treatment_lowres.tiff"))
5358. DimPlot(s_rods, reduction = "umap", label = FALSE,
5359. pt.size = .25, group.by = "batch")
5360. hirestiff(paste0(ROD_prefixPC,"_UMAP_by_batch_hires.tiff"))
5361. lowrestiff(paste0(ROD_prefixPC,"_UMAP_by_batch_lowres.tiff"))
5363. DimPlot(s_rods, reduction = "umap", label = FALSE,
5364. pt.size = .25, group.by = "batch", split.by = "batch")
5365. hirestiff(paste0(ROD_prefixPC,"_UMAP_split_by_batch_hires.tiff"))
5366. lowrestiff(paste0(ROD_prefixPC,"_UMAP_split_by_batch_lowres.tiff"))
5368. DimPlot(s_rods, reduction = "umap", label = FALSE,
5369. pt.size = .25, group.by = "time_point")
5370. hirestiff(paste0(ROD_prefixPC,"_UMAP_by_time_point_hires.tiff"))
5371. lowrestiff(paste0(ROD_prefixPC,"_UMAP_by_time_point_lowres.tiff"))
5373. DimPlot(s_rods, reduction = "umap", label = FALSE,
5374. pt.size = .25, split.by = "time_point", group.by = "time_point")
5375. hirestiff(paste0(ROD_prefixPC,"_UMAP_split_by_time_point_hires.tiff"))
5376. lowrestiff(paste0(ROD_prefixPC,"_UMAP_split_by_time_point_lowres.tiff"))
5378. DimPlot(s_rods, reduction = "umap", label = FALSE,
5379. pt.size = .25, group.by = "time_point_treatment")
5380. hirestiff(paste0(ROD_prefixPC,"_UMAP_by_time_point_treatment_hires.tiff"))
5381. lowrestiff(paste0(ROD_prefixPC,"_UMAP_by_time_point_treatment_lowres.tiff"))
5383. DimPlot(s_rods, reduction = "umap", label = FALSE,
5384. pt.size = 1, split.by = "time_point_treatment",
5385. group.by = "time_point_treatment")
5386. hirestiff(paste0(ROD_prefixPC,"_UMAP_split_by_time_point_treatment_hires.tiff"))
5387. lowrestiff(paste0(ROD_prefixPC,"_UMAP_split_by_time_point_treatment_lowres.tiff"))
5389. DimPlot(s_rods, reduction = "umap", label = FALSE, pt.size = .25, group.by = "time_point_batch")
5390. hirestiff(paste0(ROD_prefixPC,"_UMAP_by_time_point_batch_hires.tiff"))
5391. lowrestiff(paste0(ROD_prefixPC,"_UMAP_by_time_point_batch_lowres.tiff"))
5393. DimPlot(s_rods, reduction = "umap", label = FALSE, pt.size = .25, split.by = "time_point_batch", group.by = "time_point_batch")
5394. hirestiff(paste0(ROD_prefixPC,"_UMAP_split_by_time_point_batch_hires.tiff"))
5395. lowrestiff(paste0(ROD_prefixPC,"_UMAP_split_by_time_point_batch_lowres.tiff"))
5397. DimPlot(s_rods, reduction = "umap", label = FALSE, pt.size = .25, group.by = "treatment_batch")
5398. hirestiff(paste0(ROD_prefixPC,"_UMAP_by_treatment_batch_hires.tiff"))
5399. lowrestiff(paste0(ROD_prefixPC,"_UMAP_by_treatment_batch_lowres.tiff"))
5401. DimPlot(s_rods, reduction = "umap", label = FALSE, pt.size = .25, split.by = "treatment_batch", group.by = "treatment_batch")
5402. hirestiff(paste0(ROD_prefixPC,"_UMAP_split_by_treatment_batch_hires.tiff"))
5403. lowrestiff(paste0(ROD_prefixPC,"_UMAP_split_by_treatment_batch_lowres.tiff"))
5405. DimPlot(s_rods, reduction = "umap", label = FALSE, pt.size = .25, group.by = "time_point_treatment_batch")
5406. hirestiff(paste0(ROD_prefixPC,"_UMAP_by_time_point_treatment_batch_hires.tiff"))
5407. lowrestiff(paste0(ROD_prefixPC,"_UMAP_by_time_point_treatment_batch_lowres.tiff"))
5409. DimPlot(s_rods, reduction = "umap", label = FALSE, pt.size = .25, split.by = "time_point_treatment_batch", group.by = "time_point_treatment_batch")
5410. hirestiff(paste0(ROD_prefixPC,"_UMAP_split_by_time_point_treatment_batch_hires.tiff"))
5411. lowrestiff(paste0(ROD_prefixPC,"_UMAP_split_by_time_point_treatment_batch_lowres.tiff"))
5413. # saveRDS(s_rods, file = paste0(ROD_prefixPC,"_seurat_integrated_preclustering.rds"))
5415. # s_rods <- read_rds(file = paste0(ROD_prefixPC,"_seurat_integrated_preclustering.rds"))
5417. #### Remove outlier batch 08182 ####
5419. Idents(s_rods) <- "batch"
5421. s_rods <- subset(s_rods,
5422. idents = c("08182"),
5423. invert = TRUE)
5424. #### (rerun) SCT subset ####
5425. DefaultAssay(s_rods) <- "RNA"
5427. s_rods <- SCTransform(s_rods,
5428. vars.to.regress = c("percent.mt","percent.ribo"),
5429. verbose = TRUE)
5431. #### (rerun) sub PCA ####
5432. # this performs PCA on the seurat object
5433. s_rods <- RunPCA(s_rods, npcs = 50, verbose = TRUE)
5434. # make PC coordinate object a data frame
5435. xx.coord <- as.data.frame(s_rods@reductions$[email hidden])
5436. # make PC feature loadings object a data frame
5437. xx.gload <- as.data.frame(s_rods@reductions$[email hidden])
5438. # calculate eigenvalues for arrays
5439. # generate squares of all sample coordinates
5440. sq.xx.coord <- as.data.frame(xx.coord^2)
5441. # create empty list for eigenvalues first
5442. eig <- c()
5443. # calculate the eigenvalue for each PC in sq.xx.coord by taking the sqrt of the sum of squares
5444. for(i in 1:ncol(sq.xx.coord))
5445. eig[i] = sqrt(sum(sq.xx.coord[,i]))
5446. # calculate the total variance by adding up all the eigenvalues
5447. sum.eig <- sum(eig)
5448. # calculate the expected contribution of all PCs if they all contribute equally to the total variance
5449. expected.contribution <- sum.eig/(length(xx.coord)-1)
5450. # return the number of principal components with an eigenvalue greater than expected by equal variance
5451. ROD_meaningful.PCs <- sum(eig > expected.contribution)
5452. # create empty list for eigenvalue percentage
5453. eig.percent <- c()
5454. # calculate the percentage of the total variance by each PC eigenvalue
5455. for(i in 1:length(eig))
5456. eig.percent[i] = 100*eig[i]/sum.eig
5457. # sum of all eig.percent should total to 100
5458. sum(eig.percent)
5459. # create empty list for scree values
5460. scree <- c()
5461. # calculate a running total of variance contribution
5462. for(i in 1:length(eig))
5463. if(i == 1) scree[i] = eig.percent[i] else scree[i] = scree[i-1] + eig.percent[i]
5464. # create data frame for eigenvalue summaries
5465. eigenvalues <- data.frame("PC" = colnames(xx.coord), "eig" = eig, "percent" = eig.percent, "scree" = scree)
5467. # plot scree values
5468. plot(eigenvalues$percent, ylim = c(0,100), type = "S", xlab = "PC", ylab = "Percent of variance",
5469. main = paste0(date,"_",project," scree plot all samples PCA"))
5470. points(eigenvalues$scree, ylim = c(0,100), type = "p", pch = 16)
5471. lines(eigenvalues$scree)
5472. # add red line to indicate cut-off
5473. cut.off <- 100/(length(eig)-1)
5474. abline(h = cut.off, col = "red")
5475. # add blue line to indicate which PCs are meaningful and kept
5476. abline(v = ROD_meaningful.PCs, col = "blue")
5477. text(meaningful.PCs, cut.off, label = paste("cutoff PC",ROD_meaningful.PCs),
5478. adj = c(-0.1, -0.5))
5480. dev.copy(pdf, paste0("./",date,"_",project_ROD,"_scree_plot.pdf"))
5481. dev.off()
5483. rm(eigenvalues,sq.xx.coord,xx.coord,xx.gload,cut.off,
5484. eig,eig.percent,expected.contribution,i,scree,sum.eig)
5486. # ROD_meaningful.PCs <- 11
5488. #### (rerun) sub Run UMAP and look at UMAP plots ####
5490. s_rods <- RunUMAP(s_rods,
5491. reduction = "pca",
5492. dims = 1:ROD_meaningful.PCs,
5493. verbose = TRUE)
5495. # update prefixed variable
5496. ROD_prefixPC <- paste0("./",date,"_",project_ROD,"_",ROD_meaningful.PCs,"PCs")
5498. ## UMAP plot by sample name ("orig.ident")
5499. DimPlot(s_rods, reduction = "umap", label = FALSE,
5500. pt.size = .25, group.by = "orig.ident")
5501. hirestiff(paste0(ROD_prefixPC,"_UMAP_by_sample_hires.tiff"))
5502. lowrestiff(paste0(ROD_prefixPC,"_UMAP_by_sample_lowres.tiff"))
5504. DimPlot(s_rods, reduction = "umap", label = FALSE,
5505. pt.size = .25, group.by = "treatment")
5506. hirestiff(paste0(ROD_prefixPC,"_UMAP_by_treatment_hires.tiff"))
5507. lowrestiff(paste0(ROD_prefixPC,"_UMAP_by_treatment_lowres.tiff"))
5509. DimPlot(s_rods, reduction = "umap", label = FALSE,
5510. pt.size = .25, split.by = "treatment", group.by = "treatment")
5511. hirestiff(paste0(ROD_prefixPC,"_UMAP_split_by_treatment_hires.tiff"))
5512. lowrestiff(paste0(ROD_prefixPC,"_UMAP_split_by_treatment_lowres.tiff"))
5514. DimPlot(s_rods, reduction = "umap", label = FALSE,
5515. pt.size = .25, group.by = "batch")
5516. hirestiff(paste0(ROD_prefixPC,"_UMAP_by_batch_hires.tiff"))
5517. lowrestiff(paste0(ROD_prefixPC,"_UMAP_by_batch_lowres.tiff"))
5519. DimPlot(s_rods, reduction = "umap", label = FALSE,
5520. pt.size = .25, group.by = "batch", split.by = "batch")
5521. hirestiff(paste0(ROD_prefixPC,"_UMAP_split_by_batch_hires.tiff"))
5522. lowrestiff(paste0(ROD_prefixPC,"_UMAP_split_by_batch_lowres.tiff"))
5524. DimPlot(s_rods, reduction = "umap", label = FALSE,
5525. pt.size = .25, group.by = "time_point")
5526. hirestiff(paste0(ROD_prefixPC,"_UMAP_by_time_point_hires.tiff"))
5527. lowrestiff(paste0(ROD_prefixPC,"_UMAP_by_time_point_lowres.tiff"))
5529. DimPlot(s_rods, reduction = "umap", label = FALSE,
5530. pt.size = .25, split.by = "time_point", group.by = "time_point")
5531. hirestiff(paste0(ROD_prefixPC,"_UMAP_split_by_time_point_hires.tiff"))
5532. lowrestiff(paste0(ROD_prefixPC,"_UMAP_split_by_time_point_lowres.tiff"))
5534. DimPlot(s_rods, reduction = "umap", label = FALSE,
5535. pt.size = .25, group.by = "time_point_treatment")
5536. hirestiff(paste0(ROD_prefixPC,"_UMAP_by_time_point_treatment_hires.tiff"))
5537. lowrestiff(paste0(ROD_prefixPC,"_UMAP_by_time_point_treatment_lowres.tiff"))
5539. DimPlot(s_rods, reduction = "umap", label = FALSE,
5540. pt.size = 1, split.by = "time_point_treatment",
5541. group.by = "time_point_treatment")
5542. hirestiff(paste0(ROD_prefixPC,"_UMAP_split_by_time_point_treatment_hires.tiff"))
5543. lowrestiff(paste0(ROD_prefixPC,"_UMAP_split_by_time_point_treatment_lowres.tiff"))
5545. DimPlot(s_rods, reduction = "umap", label = FALSE, pt.size = .25, group.by = "time_point_batch")
5546. hirestiff(paste0(ROD_prefixPC,"_UMAP_by_time_point_batch_hires.tiff"))
5547. lowrestiff(paste0(ROD_prefixPC,"_UMAP_by_time_point_batch_lowres.tiff"))
5549. DimPlot(s_rods, reduction = "umap", label = FALSE, pt.size = .25, split.by = "time_point_batch", group.by = "time_point_batch")
5550. hirestiff(paste0(ROD_prefixPC,"_UMAP_split_by_time_point_batch_hires.tiff"))
5551. lowrestiff(paste0(ROD_prefixPC,"_UMAP_split_by_time_point_batch_lowres.tiff"))
5553. DimPlot(s_rods, reduction = "umap", label = FALSE, pt.size = .25, group.by = "treatment_batch")
5554. hirestiff(paste0(ROD_prefixPC,"_UMAP_by_treatment_batch_hires.tiff"))
5555. lowrestiff(paste0(ROD_prefixPC,"_UMAP_by_treatment_batch_lowres.tiff"))
5557. DimPlot(s_rods, reduction = "umap", label = FALSE, pt.size = .25, split.by = "treatment_batch", group.by = "treatment_batch")
5558. hirestiff(paste0(ROD_prefixPC,"_UMAP_split_by_treatment_batch_hires.tiff"))
5559. lowrestiff(paste0(ROD_prefixPC,"_UMAP_split_by_treatment_batch_lowres.tiff"))
5561. DimPlot(s_rods, reduction = "umap", label = FALSE, pt.size = .25, group.by = "time_point_treatment_batch")
5562. hirestiff(paste0(ROD_prefixPC,"_UMAP_by_time_point_treatment_batch_hires.tiff"))
5563. lowrestiff(paste0(ROD_prefixPC,"_UMAP_by_time_point_treatment_batch_lowres.tiff"))
5565. DimPlot(s_rods, reduction = "umap", label = FALSE, pt.size = .25, split.by = "time_point_treatment_batch", group.by = "time_point_treatment_batch")
5566. hirestiff(paste0(ROD_prefixPC,"_UMAP_split_by_time_point_treatment_batch_hires.tiff"))
5567. lowrestiff(paste0(ROD_prefixPC,"_UMAP_split_by_time_point_treatment_batch_lowres.tiff"))
5569. # saveRDS(s_rods, file = paste0(ROD_prefixPC,"_seurat_integrated_preclustering.rds"))
5571. # s_rods <- read_rds(file = paste0(ROD_prefixPC,"_seurat_integrated_preclustering.rds"))
5573. #### sub Clustering and Resolution ####
5574. # DefaultAssay(s_rods) <- "integrated"
5576. # Determine the K-nearest neighbor graph
5577. s_rods <- FindNeighbors(object = s_rods, reduction = "pca", dims = 1:ROD_meaningful.PCs)
5579. # Determine the clusters
5580. s_rods <- FindClusters(object = s_rods,
5581. resolution = c(0.1,0.2,0.3,0.4,0.5))
5583. ROD_res <- "_res.0.1"
5584. ROD_ident_res <- paste0("SCT_snn",ROD_res)
5586. Idents(s_rods) <- ROD_ident_res
5587. DimPlot(s_rods, reduction = "umap", label = TRUE, label.size = 5, pt.size = 0.8) +
5588. NoLegend() +
5589. ggtitle(paste0(ROD_ident_res))
5591. #update prefix
5592. ROD_prefixPCres <- paste0(ROD_prefixPC,ROD_res)
5594. # Plot the UMAP
5595. DimPlot(s_rods, reduction = "umap", label = TRUE, label.size = 6, pt.size = 1) +
5596. NoLegend() +
5597. ggtitle(paste0(ROD_ident_res))
5598. hirestiff(paste0(ROD_prefixPCres,"_UMAP","_by_","cluster","_hires.tiff"))
5599. lowrestiff(paste0(ROD_prefixPCres,"_UMAP","_by_","cluster","_lowres.tiff"))
5601. # UMAP of cells in each cluster by treatment without cluster labels
5602. DimPlot(s_rods, reduction = "umap", label = FALSE, split.by = "treatment", pt.size = 1) + NoLegend()
5603. hirestiff(paste0(ROD_prefixPCres,"_UMAP","_by_","treatment","_no_labels","_hires.tiff"))
5604. lowrestiff(paste0(ROD_prefixPCres,"_UMAP","_by_","treatment","_no_labels","_lowres.tiff"))
5606. # UMAP of cells in each cluster by treatment with cluster labels
5607. DimPlot(s_rods, reduction = "umap", label = TRUE,
5608. split.by = "treatment", pt.size = 1, label.size = 6) + NoLegend()
5609. hirestiff(paste0(ROD_prefixPCres,"_UMAP","_by_","treatment","_with_clusters_labels","_hires.tiff"))
5610. lowrestiff(paste0(ROD_prefixPCres,"_UMAP","_by_","treatment","_with_clusters_labels","_lowres.tiff"))
5612. # UMAP of cells in each cluster by batch without cluster labels
5613. DimPlot(s_rods, reduction = "umap", label = FALSE, split.by = "batch", pt.size = 1) +
5614. NoLegend() +
5615. ggtitle("Batch")
5616. hirestiff(paste0(ROD_prefixPCres,"_UMAP","_by_","batch","_no_labels","_hires.tiff"))
5617. lowrestiff(paste0(ROD_prefixPCres,"_UMAP","_by_","batch","_no_labels","_lowres.tiff"))
5619. # UMAP of cells in each cluster by batch with cluster labels
5620. DimPlot(s_rods, reduction = "umap", label = TRUE,
5621. split.by = "batch", pt.size = 1, label.size = 6) +
5622. NoLegend() +
5623. ggtitle("Batch")
5624. hirestiff(paste0(ROD_prefixPCres,"_UMAP","_by_","batch","_with_clusters_labels","_hires.tiff"))
5625. lowrestiff(paste0(ROD_prefixPCres,"_UMAP","_by_","batch","_with_clusters_labels","_lowres.tiff"))
5627. # UMAP of cells in each cluster by time_point without cluster labels
5628. DimPlot(s_rods, reduction = "umap", label = FALSE, split.by = "time_point", pt.size = 1) + NoLegend()
5629. hirestiff(paste0(ROD_prefixPCres,"_UMAP","_by_","time_point","_no_labels","_hires.tiff"))
5630. lowrestiff(paste0(ROD_prefixPCres,"_UMAP","_by_","time_point","_no_labels","_lowres.tiff"))
5632. # UMAP of cells in each cluster by time_point with cluster labels
5633. DimPlot(s_rods, reduction = "umap", label = TRUE,
5634. split.by = "time_point", pt.size = 1, label.size = 6) + NoLegend()
5635. hirestiff(paste0(ROD_prefixPCres,"_UMAP","_by_","time_point","_with_clusters_labels","_hires.tiff"))
5636. lowrestiff(paste0(ROD_prefixPCres,"_UMAP","_by_","time_point","_with_clusters_labels","_lowres.tiff"))
5638. # UMAP of cells in each cluster by time_point without cluster labels
5639. DimPlot(s_rods, reduction = "umap", label = FALSE, split.by = "time_point_treatment", pt.size = 1) + NoLegend()
5640. hirestiff(paste0(ROD_prefixPCres,"_UMAP","_by_","time_point_treatment","_no_labels","_hires.tiff"))
5641. lowrestiff(paste0(ROD_prefixPCres,"_UMAP","_by_","time_point_treatment","_no_labels","_lowres.tiff"))
5643. # UMAP of cells in each cluster by time_point with cluster labels
5644. DimPlot(s_rods, reduction = "umap", label = TRUE,
5645. split.by = "time_point_treatment", pt.size = 1, label.size = 6) + NoLegend()
5646. hirestiff(paste0(ROD_prefixPCres,"_UMAP","_by_","time_point_treatment","_with_clusters_labels","_hires.tiff"))
5647. lowrestiff(paste0(ROD_prefixPCres,"_UMAP","_by_","time_point_treatment","_with_clusters_labels","_lowres.tiff"))
5649. #### sub Extract number of cells per cluster per orig.ident ####
5650. n_cells <- FetchData(s_rods, vars = c("ident", "orig.ident")) %>%
5651. dplyr::count(ident, orig.ident) %>%
5652. tidyr::spread(ident, n)
5653. n_cells
5654. write.csv(n_cells, file = paste0(ROD_prefixPCres,"_cells_per_cluster.csv"))
5656. rm(n_cells)
5658. #### sub save RDS containing reduction and cluster idents ####
5659. # saveRDS(s_rods, paste0(ROD_prefixPCres,"_seurat_after_clustering.rds"))
5661. # s_rods <- readRDS(file = paste0(ROD_prefixPCres,"_seurat_after_clustering.rds"))
5663. #### sub normalize rna slot ####
5665. # Select the RNA counts slot to be the default assay for visualization purposes
5666. DefaultAssay(s_rods) <- "RNA"
5668. # Normalize, find variable features, scale data
5669. s_rods <- NormalizeData(s_rods)
5670. s_rods <- FindVariableFeatures(s_rods)
5671. all.genes <- rownames(s_rods)
5672. s_rods <- ScaleData(s_rods, features = all.genes)
5674. # Export normalized counts
5675. # norm_counts <- GetAssayData(s_rods, slot = "data")
5676. # save(norm_counts, file = paste0(ROD_prefixPCres,"_normalized_counts_RNA_sparse_matrix.rdata"))
5678. # write.csv(norm_counts, file = paste0(ROD_prefixPCres,"_normalized_counts_RNA.csv"))
5680. # save RDS containing RNA normalized data
5681. # saveRDS(s_rods, paste0(ROD_prefixPCres,"_seurat_after_RNAnorm.rds"))
5682. # s_rods <- readRDS(file = paste0(ROD_prefixPCres,"_seurat_after_RNAnorm.rds"))
5684. #### Plots ####
5686. # PLOTS
5688. VlnPlot(s_rods,
5689. features = Rods,
5690. stack = TRUE,
5691. flip = TRUE) +
5692. NoLegend() +
5693. ggtitle("Rod Photoreceptors")
5695. hirestiff(paste(ROD_prefixPCres,"VLN","Rod","genes","hires.tiff", sep = "_"))
5696. lowrestiff(paste(ROD_prefixPCres,"VLN","Rod","genes","lowres.tiff", sep = "_"))
5698. FeaturePlot(s_rods,
5699. reduction = "umap",
5700. features = c(Rods1),
5701. order = TRUE,
5702. min.cutoff = 'q10',
5703. label = TRUE,
5704. label.size = 5,
5705. pt.size = 0.8)
5707. hirestiff(paste(ROD_prefixPCres,"FTR","Rod","genes","1","hires.tiff", sep = "_"))
5708. lowrestiff(paste(ROD_prefixPCres,"FTR","Rod","genes","1","lowres.tiff", sep = "_"))
5710. FeaturePlot(s_rods,
5711. reduction = "umap",
5712. features = c(Rods2),
5713. order = TRUE,
5714. min.cutoff = 'q10',
5715. label = TRUE,
5716. label.size = 5,
5717. pt.size = 0.8)
5719. hirestiff(paste(ROD_prefixPCres,"FTR","Rod","genes","2","hires.tiff", sep = "_"))
5720. lowrestiff(paste(ROD_prefixPCres,"FTR","Rod","genes","2","lowres.tiff", sep = "_"))
5722. #### DEGs ####
5724. #### Find Cluster Markers for ####
5726. # find all markers
5727. cluster_markers <- FindAllMarkers(s_rods, only.pos = TRUE, min.pct = 0.50, logfc.threshold = 0.4)
5729. write.csv(cluster_markers,
5730. file = paste(ROD_prefixPCres,"clustermarkers_pos_only_minpct0.50_logfcthresh0.4.csv", sep = "_"))
5732. #### find DEGs between treatments ####
5734. Idents(object = s_rods) <- "treatment"
5736. cluster_markers <- FindAllMarkers(s_rods, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.4)
5738. write.csv(cluster_markers, file =
5739. paste(ROD_prefixPCres,"degs_treatment_all_cells_minpct0.25_logfcthresh0.4.csv", sep = "_"))
5741. Idents(object = s_rods) <- ident_res
5743. #### find DEGs between time points ####
5745. Idents(object = s_rods) <- "time_point"
5747. cluster_markers <- FindAllMarkers(s_rods, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.4)
5749. write.csv(cluster_markers, file =
5750. paste(ROD_prefixPCres,"degs_timepoint_all_cells_minpct0.25_logfcthresh0.4.csv", sep = "_"))
5752. Idents(object = s_rods) <- ident_res
5754. #### find DEGs between time points and treatment ####
5756. Idents(object = s_rods) <- "time_point_treatment"
5758. cluster_markers <- FindAllMarkers(s_rods, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.4)
5760. write.csv(cluster_markers, file =
5761. paste(ROD_prefixPCres,"degs_timepoint_treatment_all_cells_minpct0.25_logfcthresh0.4.csv", sep = "_"))
5763. Idents(object = s_rods) <- ident_res
5765. #### Cone subset recluster ####
5766. s_cones <- subset(seurat_sct,
5767. idents = c("CONE9"))
5769. #### SCT subset ####
5770. DefaultAssay(s_cones) <- "RNA"
5772. s_cones <- SCTransform(s_cones,
5773. vars.to.regress = c("percent.mt","percent.ribo"),
5774. verbose = TRUE)
5777. #### update project variable to include subset ####
5778. project <- "Rat_retina_CONES"
5780. #### sub PCA ####
5781. # this performs PCA on the seurat object
5782. s_cones <- RunPCA(s_cones, npcs = 50, verbose = TRUE)
5784. # make PC coordinate object a data frame
5785. xx.coord <- as.data.frame(s_cones@reductions$[email hidden])
5787. # make PC feature loadings object a data frame
5788. xx.gload <- as.data.frame(s_cones@reductions$[email hidden])
5790. # calculate eigenvalues for arrays
5792. # generate squares of all sample coordinates
5793. sq.xx.coord <- as.data.frame(xx.coord^2)
5794. # create empty list for eigenvalues first
5795. eig <- c()
5796. # calculate the eigenvalue for each PC in sq.xx.coord by taking the sqrt of the sum of squares
5797. for(i in 1:ncol(sq.xx.coord))
5798. eig[i] = sqrt(sum(sq.xx.coord[,i]))
5799. # calculate the total variance by adding up all the eigenvalues
5800. sum.eig <- sum(eig)
5801. # calculate the expected contribution of all PCs if they all contribute equally to the total variance
5802. expected.contribution <- sum.eig/(length(xx.coord)-1)
5803. # return the number of principal components with an eigenvalue greater than expected by equal variance
5804. meaningful.PCs <- sum(eig > expected.contribution)
5806. # create empty list for eigenvalue percentage
5807. eig.percent <- c()
5808. # calculate the percentage of the total variance by each PC eigenvalue
5809. for(i in 1:length(eig))
5810. eig.percent[i] = 100*eig[i]/sum.eig
5811. # sum of all eig.percent should total to 100
5812. sum(eig.percent)
5813. # create empty list for scree values
5814. scree <- c()
5815. # calculate a running total of variance contribution
5816. for(i in 1:length(eig))
5817. if(i == 1) scree[i] = eig.percent[i] else scree[i] = scree[i-1] + eig.percent[i]
5819. # create data frame for eigenvalue summaries
5820. eigenvalues <- data.frame("PC" = colnames(xx.coord), "eig" = eig, "percent" = eig.percent, "scree" = scree)
5822. # write csv for eigenvalues
5823. # write.csv(eigenvalues, file = paste0("./",date,"_",project,"_PCA_eigenvalues.csv"), row.names = F)
5825. # plot scree values
5826. plot(eigenvalues$percent, ylim = c(0,100), type = "S", xlab = "PC", ylab = "Percent of variance",
5827. main = paste0(date,"_",project," scree plot all samples PCA"))
5828. points(eigenvalues$scree, ylim = c(0,100), type = "p", pch = 16)
5829. lines(eigenvalues$scree)
5830. # add red line to indicate cut-off
5831. cut.off <- 100/(length(eig)-1)
5832. abline(h = cut.off, col = "red")
5833. # add blue line to indicate which PCs are meaningful and kept
5834. abline(v = meaningful.PCs, col = "blue")
5835. text(meaningful.PCs, cut.off, label = paste("cutoff PC",meaningful.PCs),
5836. adj = c(-0.1, -0.5))
5838. dev.copy(pdf, paste0("./",date,"_",project,"_scree_plot.pdf"))
5839. dev.off()
5841. rm(eigenvalues,sq.xx.coord,xx.coord,xx.gload,cut.off,
5842. eig,eig.percent,expected.contribution,i,scree,sum.eig)
5844. # meaningful.PCs <- 10
5846. #### sub Run UMAP and look at UMAP plots ####
5848. s_cones <- RunUMAP(s_cones,
5849. reduction = "pca",
5850. dims = 1:meaningful.PCs,
5851. verbose = TRUE)
5853. # update prefixed variable
5854. prefixPC <- paste0("./",date,"_",project,"_",meaningful.PCs,"PCs")
5856. ## UMAP plot by sample name ("orig.ident")
5857. DimPlot(s_cones, reduction = "umap", label = FALSE,
5858. pt.size = .25, group.by = "orig.ident")
5859. hirestiff(paste0(prefixPC,"_UMAP_by_sample_hires.tiff"))
5860. lowrestiff(paste0(prefixPC,"_UMAP_by_sample_lowres.tiff"))
5862. DimPlot(s_cones, reduction = "umap", label = FALSE,
5863. pt.size = .25, group.by = "treatment")
5864. hirestiff(paste0(prefixPC,"_UMAP_by_treatment_hires.tiff"))
5865. lowrestiff(paste0(prefixPC,"_UMAP_by_treatment_lowres.tiff"))
5867. DimPlot(s_cones, reduction = "umap", label = FALSE,
5868. pt.size = .25, split.by = "treatment", group.by = "treatment")
5869. hirestiff(paste0(prefixPC,"_UMAP_split_by_treatment_hires.tiff"))
5870. lowrestiff(paste0(prefixPC,"_UMAP_split_by_treatment_lowres.tiff"))
5872. DimPlot(s_cones, reduction = "umap", label = FALSE,
5873. pt.size = .25, group.by = "batch")
5874. hirestiff(paste0(prefixPC,"_UMAP_by_batch_hires.tiff"))
5875. lowrestiff(paste0(prefixPC,"_UMAP_by_batch_lowres.tiff"))
5877. DimPlot(s_cones, reduction = "umap", label = FALSE,
5878. pt.size = .25, group.by = "batch", split.by = "batch")
5879. hirestiff(paste0(prefixPC,"_UMAP_split_by_batch_hires.tiff"))
5880. lowrestiff(paste0(prefixPC,"_UMAP_split_by_batch_lowres.tiff"))
5882. DimPlot(s_cones, reduction = "umap", label = FALSE,
5883. pt.size = .25, group.by = "time_point")
5884. hirestiff(paste0(prefixPC,"_UMAP_by_time_point_hires.tiff"))
5885. lowrestiff(paste0(prefixPC,"_UMAP_by_time_point_lowres.tiff"))
5887. DimPlot(s_cones, reduction = "umap", label = FALSE,
5888. pt.size = .25, split.by = "time_point", group.by = "time_point")
5889. hirestiff(paste0(prefixPC,"_UMAP_split_by_time_point_hires.tiff"))
5890. lowrestiff(paste0(prefixPC,"_UMAP_split_by_time_point_lowres.tiff"))
5892. DimPlot(s_cones, reduction = "umap", label = FALSE,
5893. pt.size = .25, group.by = "time_point_treatment")
5894. hirestiff(paste0(prefixPC,"_UMAP_by_time_point_treatment_hires.tiff"))
5895. lowrestiff(paste0(prefixPC,"_UMAP_by_time_point_treatment_lowres.tiff"))
5897. DimPlot(s_cones, reduction = "umap", label = FALSE,
5898. pt.size = 1, split.by = "time_point_treatment",
5899. group.by = "time_point_treatment")
5900. hirestiff(paste0(prefixPC,"_UMAP_split_by_time_point_treatment_hires.tiff"))
5901. lowrestiff(paste0(prefixPC,"_UMAP_split_by_time_point_treatment_lowres.tiff"))
5903. DimPlot(s_cones, reduction = "umap", label = FALSE, pt.size = .25, group.by = "time_point_batch")
5904. hirestiff(paste0(prefixPC,"_UMAP_by_time_point_batch_hires.tiff"))
5905. lowrestiff(paste0(prefixPC,"_UMAP_by_time_point_batch_lowres.tiff"))
5907. DimPlot(s_cones, reduction = "umap", label = FALSE, pt.size = .25, split.by = "time_point_batch", group.by = "time_point_batch")
5908. hirestiff(paste0(prefixPC,"_UMAP_split_by_time_point_batch_hires.tiff"))
5909. lowrestiff(paste0(prefixPC,"_UMAP_split_by_time_point_batch_lowres.tiff"))
5911. DimPlot(s_cones, reduction = "umap", label = FALSE, pt.size = .25, group.by = "treatment_batch")
5912. hirestiff(paste0(prefixPC,"_UMAP_by_treatment_batch_hires.tiff"))
5913. lowrestiff(paste0(prefixPC,"_UMAP_by_treatment_batch_lowres.tiff"))
5915. DimPlot(s_cones, reduction = "umap", label = FALSE, pt.size = .25, split.by = "treatment_batch", group.by = "treatment_batch")
5916. hirestiff(paste0(prefixPC,"_UMAP_split_by_treatment_batch_hires.tiff"))
5917. lowrestiff(paste0(prefixPC,"_UMAP_split_by_treatment_batch_lowres.tiff"))
5919. DimPlot(s_cones, reduction = "umap", label = FALSE, pt.size = .25, group.by = "time_point_treatment_batch")
5920. hirestiff(paste0(prefixPC,"_UMAP_by_time_point_treatment_batch_hires.tiff"))
5921. lowrestiff(paste0(prefixPC,"_UMAP_by_time_point_treatment_batch_lowres.tiff"))
5923. DimPlot(s_cones, reduction = "umap", label = FALSE, pt.size = .25, split.by = "time_point_treatment_batch", group.by = "time_point_treatment_batch")
5924. hirestiff(paste0(prefixPC,"_UMAP_split_by_time_point_treatment_batch_hires.tiff"))
5925. lowrestiff(paste0(prefixPC,"_UMAP_split_by_time_point_treatment_batch_lowres.tiff"))
5927. # saveRDS(s_cones, file = paste0(prefixPC,"_seurat_integrated_preclustering.rds"))
5929. # s_cones <- read_rds(file = paste0(prefixPC,"_seurat_integrated_preclustering.rds"))
5931. #### sub Clustering and Resolution ####
5932. # DefaultAssay(s_cones) <- "integrated"
5934. # Determine the K-nearest neighbor graph
5935. s_cones <- FindNeighbors(object = s_cones, reduction = "pca", dims = 1:meaningful.PCs)
5937. # Determine the clusters
5938. s_cones <- FindClusters(object = s_cones,
5939. resolution = c(0.1,0.2,0.3,0.4,0.5))
5940. view([email hidden])
5941. res <- "_res.0.1"
5942. ident_res <- paste0("SCT_snn",res)
5944. Idents(s_cones) <- ident_res
5945. DimPlot(s_cones, reduction = "umap", label = TRUE, label.size = 5, pt.size = 0.8) +
5946. NoLegend() +
5947. ggtitle(paste0(ident_res))
5949. #update prefix
5950. prefixPCres <- paste0(prefixPC,res)
5952. # Plot the UMAP
5953. DimPlot(s_cones, reduction = "umap", label = TRUE, label.size = 6, pt.size = 0.8) +
5954. NoLegend() +
5955. ggtitle(paste0(ident_res))
5956. hirestiff(paste0(prefixPCres,"_UMAP","_by_","cluster","_hires.tiff"))
5957. lowrestiff(paste0(prefixPCres,"_UMAP","_by_","cluster","_lowres.tiff"))
5959. # UMAP of cells in each cluster by treatment without cluster labels
5960. DimPlot(s_cones, reduction = "umap", label = FALSE,

20230922_rat_retina_with_GEO.R at commit 0ca5d25, under Apache-2.0 · at the source