-
Notifications
You must be signed in to change notification settings - Fork 0
/
3_DApeaksInTBC.Rmd
119 lines (102 loc) · 5.15 KB
/
3_DApeaksInTBC.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
---
title: "scATAC-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 scATAC-seq data, generated from rat metrial gland tissues on gestational day (GD) 15.5 and 19.5. <br>
## GD 15.5:
```{r}
library(Signac)
library(Seurat)
library(ggplot2)
library(GenomicRanges)
library(patchwork)
set.seed(1234)
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/rnor6.rda")
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/rnor6.ranges.rda")
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd15.5-ATACwithRNAlables.rda")
Idents(rats) <- [email protected]$predicted.id # change from old identities to those predicted by scRNA-seq
DefaultAssay(rats) <- 'MACSpeaks'
Annotation(rats) <- rnor6.ranges
```
```{r,eval=F}
rats1.markers <- FindAllMarkers(rats, test.use = "LR", latent.vars = 'peak_region_fragments')
#write.table(rats1.markers, "/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd15.5-unfilteredMarkers.txt", quote = F, sep = '\t')
```
```{r}
rats1.markers <- read.table("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd15.5-unfilteredMarkers.txt", header = T, sep = "\t")
tbc15.5 <- rats1.markers[rats1.markers$cluster == "Invasive trophoblast cells" & rats1.markers$avg_log2FC >= log2(1.5) & rats1.markers$p_val_adj <= 0.05,]
#write.table(tbc15.5, "/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd15.5-tbc.txt", quote = F, sep = '\t')
dim(tbc15.5)
closest_genes <- ClosestFeature(rats, regions = rownames(tbc15.5))
length(unique(closest_genes$gene_name))
#write.table(closest_genes, "/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd15.5-tbc-toGenes.txt", quote = F, sep = '\t')
```
In BASH, run the following command to obtain a bed file of the TBC peaks:
```
cut -f 8 gd15.5-tbc.txt | tail -n +2 | sed 's/-/\t/g' | awk '{print $1"\t"$2"\t"$3"\t"$1"-"$2"-"$3}' > gd15.5-tbc.bed
```
We can look at the distribution of the TBC-enriched peaks:
```{r}
library(ChIPseeker)
library(GenomicFeatures)
library(RMariaDB)
txdb <- makeTxDbFromEnsembl(organism= "Rattus norvegicus", release = 98)
file <- "/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd15.5-tbc.bed"
peakAnno <- annotatePeak(file, tssRegion=c(-3000, 3000), TxDb=txdb)
#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS3-15.5-tbc.#pdf", height = 5, width = 6)
plotDistToTSS(peakAnno,
title="Distribution of transcription factor-binding loci relative to TSS")
#dev.off()
```
## GD 19.5:
```{r}
set.seed(1234)
load("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd19.5-ATACwithRNAlables.rda")
Idents(rats) <- [email protected]$predicted.id # change from old identities to those predicted by scRNA-seq
DefaultAssay(rats) <- 'MACSpeaks'
Annotation(rats) <- rnor6.ranges
```
```{r, eval=F}
rats1.markers <- FindAllMarkers(rats, test.use = "LR", latent.vars = 'peak_region_fragments')
#write.table(rats1.markers, "/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd19.5-unfilteredMarkers.txt", quote = F, sep = '\t')
```
```{r}
rats1.markers <- read.table("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd19.5-unfilteredMarkers.txt", header = T, sep = "\t")
tbc19.5 <- rats1.markers[rats1.markers$cluster == "Invasive trophoblast cells" & rats1.markers$avg_log2FC >= log2(1.5) & rats1.markers$p_val_adj <= 0.05,]
dim(tbc19.5)
#write.table(tbc19.5, "/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd19.5-tbc.txt", quote = F, sep = '\t')
closest_genes <- ClosestFeature(rats, regions = rownames(tbc19.5))
length(unique(closest_genes$gene_name))
#write.table(closest_genes, "/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd19.5-tbc-toGenes.txt", quote = F, sep = '\t')
```
In BASH:
```
cut -f 1 gd19.5-tbc.txt | tail -n +2 | sed 's/-/\t/g' | awk '{print $1"\t"$2"\t"$3"\t"$1"-"$2"-"$3}' > gd19.5-tbc.bed
```
```{r}
file <- "/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/4_RNAintegration/2_MACS2Peaks/annotation_scRNAseq/gd19.5-tbc.bed"
peakAnno <- annotatePeak(file, tssRegion=c(-3000, 3000), TxDb=txdb)
#pdf("/work/LAS/geetu-lab/hhvu/project3_scATAC/scATAC-seq-analysis/FIGURES/20220829_draftV3/figureS3-19.5-tbc.#pdf", height = 5, width = 6)
plotDistToTSS(peakAnno,
title="Distribution of transcription factor-binding loci relative to TSS")
#dev.off()
```
```{r}
sessionInfo()
```