-
Notifications
You must be signed in to change notification settings - Fork 0
/
4_DEG_GD15.5vsGD19.5.Rmd
103 lines (80 loc) · 4.74 KB
/
4_DEG_GD15.5vsGD19.5.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
---
title: "scRNA-seq Rat Metrial Glands"
author: "Ha T. H. Vu"
output: html_document
---
```{r setup, include=FALSE}
options(max.print = "75")
knitr::opts_chunk$set(
echo = TRUE,
collapse = TRUE,
comment = "#>",
fig.path = "Files/",
fig.width = 15,
prompt = FALSE,
tidy = FALSE,
message = FALSE,
warning = TRUE
)
knitr::opts_knit$set(width = 75)
```
This is a documentation for analyses of scRNA-seq data, generated from rat metrial gland tissues on gestational day (GD) 15.5 and 19.5. <br>
```{r}
library(dplyr)
library(Seurat)
library(patchwork)
library(ggplot2)
set.seed(1234)
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scRNA-seq-analysis/3_cluster_DEG/rmChrMTgenes/GD15.5/gd15.5_res0.8.rda")
gd15.5 <- data
gd15.5$sample <- substr(names(gd15.5$orig.ident), 1, 6)
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scRNA-seq-analysis/3_cluster_DEG/rmChrMTgenes/GD19.5/gd19.5-4-5-6-7_res0.8.rda")
gd19.5 <- data2
gd19.5$sample <- ifelse(names(gd19.5$orig.ident) == "19.5_7", names(gd19.5$orig.ident),
substr(names(gd19.5$orig.ident), 1, 6))
DefaultAssay(gd15.5) <- "RNA"
DefaultAssay(gd19.5) <- "RNA"
gd15.5 <- RenameIdents(gd15.5, `0` = "gd15.5_0", `1` = "gd15.5_1", `2` = "gd15.5_2", `3` = "gd15.5_3", `4` = "gd15.5_4", `5` = "gd15.5_5", `6` = "gd15.5_6", `7` = "gd15.5_7", `8` = "gd15.5_8", `9` = "gd15.5_9", `10` = "gd15.5_10", `11` = "gd15.5_11", `12` = "gd15.5_12", `13` = "gd15.5_13", `14` = "gd15.5_14", `15` = "gd15.5_15", `16` = "gd15.5_16", `17` = "gd15.5_17", `18` = "gd15.5_18", `19` = "gd15.5_19", `20` = "gd15.5_20", `21` = "gd15.5_21", `22` = "gd15.5_22", `23` = "gd15.5_23-TBC", `24` = "gd15.5_24", `25` = "gd15.5_25", `26` = "gd15.5_26", `27` = "gd15.5_27", `28` = "gd15.5_28")
gd19.5 <- RenameIdents(gd19.5, `0` = "gd19.5_0", `1` = "gd19.5_1", `2` = "gd19.5_2", `3` = "gd19.5_3", `4` = "gd19.5_4", `5` = "gd19.5_5", `6` = "gd19.5_6-TBC", `7` = "gd19.5_7", `8` = "gd19.5_8", `9` = "gd19.5_9", `10` = "gd19.5_10", `11` = "gd19.5_11", `12` = "gd19.5_12", `13` = "gd19.5_13", `14` = "gd19.5_14", `15` = "gd19.5_15", `16` = "gd19.5_16", `17` = "gd19.5_17", `18` = "gd19.5_18", `19` = "gd19.5_19", `20` = "gd19.5_20", `21` = "gd19.5_21", `22` = "gd19.5_22", `23` = "gd19.5_23", `24` = "gd19.5_24", `25` = "gd19.5_25", `26` = "gd19.5_26")
twotimepoint <- merge(gd15.5, gd19.5)
DefaultAssay(twotimepoint) <- "RNA"
#save(twotimepoint, file = "/work/LAS/geetu-lab/hhvu/project3_scATAC/scRNA-seq-analysis/4_twoTimepoints/rmChrMTgenes/twotimepoint.rda")
trophoblast.markers <- FindMarkers(twotimepoint, ident.1 = "gd19.5_6-TBC", ident.2 = "gd15.5_23-TBC", verbose = FALSE, logfc.threshold = 0)
#save(trophoblast.markers, file = "/work/LAS/geetu-lab/hhvu/project3_scATAC/scRNA-seq-analysis/4_twoTimepoints/rmChrMTgenes/trophoblast.markers.rda")
gd15.5tbc <- subset(trophoblast.markers, trophoblast.markers$p_val_adj <= 0.05 & trophoblast.markers$avg_log2FC <= -log2(1.5))
dim(gd15.5tbc)
head(gd15.5tbc)
#write.table(gd15.5tbc, "/work/LAS/geetu-lab/hhvu/project3_scATAC/scRNA-seq-analysis/4_twoTimepoints/rmChrMTgenes/DEG/gd15.5TBC.txt", sep = "\t", row.names = T, quote = F)
gd19.5tbc <- subset(trophoblast.markers, trophoblast.markers$p_val_adj <= 0.05 & trophoblast.markers$avg_log2FC >= log2(1.5))
dim(gd19.5tbc)
head(gd19.5tbc)
#write.table(gd19.5tbc, "/work/LAS/geetu-lab/hhvu/project3_scATAC/scRNA-seq-analysis/4_twoTimepoints/rmChrMTgenes/DEG/gd19.5TBC.txt", sep = "\t", row.names = T, quote = F)
```
We can assess the functional relevance of the differentially expressed genes at each timepoint using gene ontology analysis.
```{r}
library("org.Rn.eg.db")
library("clusterProfiler")
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/rnor6.rda")
go <- function(genes) {
GO <- enrichGO(gene = genes, OrgDb=org.Rn.eg.db, ont = "BP", pvalueCutoff = 0.05, qvalueCutoff = 0.05, maxGSSize=500, readable = T, keyType = "ENSEMBL")
goSimplify <- simplify(GO)
res <- as.data.frame(goSimplify)
res$Rank <- seq(1, nrow(res))
res$GeneRatio <- sapply(strsplit(res$GeneRatio, "/"), function(x) as.numeric(x[1])/as.numeric(x[2]))
res$BgRatio <- sapply(strsplit(res$BgRatio, "/"), function(x) as.numeric(x[1])/as.numeric(x[2]))
res$Fold <- as.numeric(res$GeneRatio)/as.numeric(res$BgRatio)
res <- res[res$qvalue <= 0.05 & res$Fold >= 2 & res$Count >= 5,]
return(res)
}
g <- rnor6[rnor6$gene_name %in% rownames(gd15.5tbc),]
r <- go(g$ensembl_gene_id)
head(r)
#save(r, file="/work/LAS/geetu-lab/hhvu/project3_scATAC/scRNA-seq-analysis/4_twoTimepoints/rmChrMTgenes/DEG/gd15.5TBC_goRats-max500-simplified.rda")
g <- rnor6[rnor6$gene_name %in% rownames(gd19.5tbc),]
r <- go(g$ensembl_gene_id)
head(r)
#save(r, file="/work/LAS/geetu-lab/hhvu/project3_scATAC/scRNA-seq-analysis/4_twoTimepoints/rmChrMTgenes/DEG/gd19.5TBC_goRats-max500-simplified.rda")
```
```{r}
sessionInfo()
```