Mercurial > repos > gaelcge > r_signac_galaxy

comparison Signac.R @ 0:6e0b320d8b6a draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"planemo upload commit dc808171975d0012e25bd7b32adc7a5a5c56a145-dirty"
author gaelcge
date Tue, 02 Aug 2022 19:11:27 +0000
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:6e0b320d8b6a
1 rm(list = ls())
2
3 library(Signac)
4 library(Seurat)
5 library(GenomeInfoDb)
6 library(EnsDb.Hsapiens.v75)
7 library(ggplot2)
8 library(patchwork)
9
10 library(future)
11 plan("multicore", workers = (availableCores()-1))
12 options(future.globals.maxSize = 3000000 * 1024^2)
13
14 setwd("/home/gaelcge/projects/def-jsjoyal/gaelcge/scATACseq/10XData/atac_v1_pbmc_10k_output/")
15
16 counts <- Read10X_h5(filename = "./atac_v1_pbmc_10k_filtered_peak_bc_matrix.h5")
17
18
19 metadata <- read.csv(
20 file = "./atac_v1_pbmc_10k_singlecell.csv",
21 header = TRUE,
22 row.names = 1
23 )
24
25 metadata <- metadata[colnames(counts),]
26
27 chrom_assay <- CreateChromatinAssay(
28 counts = counts,
29 sep = c(":", "-"),
30 genome = 'hg19',
31 fragments = './atac_v1_pbmc_10k_fragments.tsv.gz',
32 min.cells = 10,
33 min.features = 200
34 )
35
36 pbmc <- CreateSeuratObject(
37 counts = chrom_assay,
38 assay = "peaks",
39 meta.data = metadata
40 )
41
42 setwd("/home/gaelcge/projects/def-jsjoyal/gaelcge/scATACseq/Signac_analysis/atac_pbmc_500_nextgem")
43
44 pbmc[['peaks']]
45
46 granges(pbmc)
47
48 # extract gene annotations from EnsDb
49 annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v75)
50
51 # change to UCSC style since the data was mapped to GRCh38
52 seqlevelsStyle(annotations) <- 'UCSC'
53 genome(annotations) <- "GRCh38"
54
55 # add the gene information to the object
56 Annotation(pbmc) <- annotations
57
58
59 # compute nucleosome signal score per cell
60 pbmc <- NucleosomeSignal(object = pbmc)
61
62 # compute TSS enrichment score per cell
63 pbmc <- TSSEnrichment(object = pbmc, fast = FALSE)
64
65 # add blacklist ratio and fraction of reads in peaks
66 pbmc$pct_reads_in_peaks <- pbmc$peak_region_fragments / pbmc$passed_filters * 100
67 pbmc$blacklist_ratio <- pbmc$blacklist_region_fragments / pbmc$peak_region_fragments
68
69
70
71 pbmc$high.tss <- ifelse(pbmc$TSS.enrichment > 2, 'High', 'Low')
72
73 png("Tssplot.png")
74 TSSPlot(pbmc, group.by = 'high.tss') + NoLegend()
75 dev.off()
76
77
78
79 pbmc$nucleosome_group <- ifelse(pbmc$nucleosome_signal > 4, 'NS > 4', 'NS < 4')
80 png("FragmentHistogram.png")
81 FragmentHistogram(object = pbmc, group.by = 'nucleosome_group')
82 dev.off()
83
84 png("VlnPlot_QC.png", width=1000)
85 VlnPlot(
86 object = pbmc,
87 features = c('pct_reads_in_peaks', 'peak_region_fragments',
88 'TSS.enrichment', 'blacklist_ratio', 'nucleosome_signal'),
89 pt.size = 0.1,
90 ncol = 5
91 )
92 dev.off()
93
94
95 pbmc <- subset(
96 x = pbmc,
97 subset = peak_region_fragments > 3000 &
98 peak_region_fragments < 20000 &
99 pct_reads_in_peaks > 15 &
100 blacklist_ratio < 0.05 &
101 nucleosome_signal < 2 &
102 TSS.enrichment > 1
103 )
104 pbmc
105
106
107 pbmc <- RunTFIDF(pbmc)
108 pbmc <- FindTopFeatures(pbmc, min.cutoff = 'q0')
109 pbmc <- RunSVD(pbmc)
110
111 png("DepthCor.png")
112 DepthCor(pbmc)
113 dev.off()
114
115 pbmc <- RunUMAP(object = pbmc, reduction = 'lsi', dims = 2:30)
116 pbmc <- FindNeighbors(object = pbmc, reduction = 'lsi', dims = 2:30)
117 pbmc <- FindClusters(object = pbmc, verbose = FALSE, algorithm = 3)
118
119
120 png("UMAP.png")
121 DimPlot(object = pbmc, label = TRUE) + NoLegend()
122 dev.off()
123
124
125 gene.activities <- GeneActivity(pbmc)
126
127
128 # add the gene activity matrix to the Seurat object as a new assay and normalize it
129 pbmc[['RNA']] <- CreateAssayObject(counts = gene.activities)
130 pbmc <- NormalizeData(
131 object = pbmc,
132 assay = 'RNA',
133 normalization.method = 'LogNormalize',
134 scale.factor = median(pbmc$nCount_RNA)
135 )
136
137
138 DefaultAssay(pbmc) <- 'RNA'
139
140 png("FeaturePlot_knownMarkers.png", width=1000)
141 FeaturePlot(
142 object = pbmc,
143 features = c('MS4A1', 'CD3D', 'LEF1', 'NKG7', 'TREM1', 'LYZ'),
144 pt.size = 0.1,
145 max.cutoff = 'q95',
146 ncol = 3
147 )
148 dev.off()
149
150
151 saveRDS(pbmc, paste("Seurat_object.rds", sep="."))
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170 q("no")