SingleCellAnalysisTutorial/scRNAseqAnalysis.R at main · EngreitzLab/SingleCellAnalysisTutorial · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
## Single Cell Analysis Tutorial
## 220713

## Data used in this tutorial: Multiome hESC-EC data generated by Dulguun and Chad
## Also see: https://satijalab.org/seurat/articles/pbmc3k_tutorial.html

## load packages
suppressPackageStartupMessages({
  library(optparse)
  library(dplyr)
  library(tidyr)
  library(reshape2)
  library(ggplot2)
  library(cowplot)
  library(ggpubr) ## ggarrange
  library(gplots) ## heatmap.2
  library(scales) ## geom_tile gradient rescale
  library(ggrepel)
  library(stringr)
  library(Seurat)
  library(SeuratObject)
  # library(ggseqlogo)
  # library(universalmotif)
})


##########################################################################################
## parse in command line options
option.list <- list(
  make_option("--datadir",type="character", default="/Volumes/engreitz/Users/kangh/process_sequencing_data/210109_scATAC.GEX/", help="Data directory"),
  make_option("--sampleName",type="character",default="multiome_fresh", help="Sample name"),
  make_option("--geneSet", type="character",default="/Volumes/engreitz/Users/kangh/multiome_hESC.EC_fresh.frozen/data/marker.gene.list.xlsx", help="Gene set file location"),
  make_option("--maxMt", type="numeric", default=50, help="filter out cells with percent mitochondrial gene higher than this threhsold"),
  make_option("--maxCount", type="numeric", default=25000, help="filter out cells with UMI count more than this threshold"),
  make_option("--minUniqueGenes", type="numeric", default=0, help="filter out cells with unique gene detected less than this threshold"),
  make_option("--project",type="character",default="/Volumes/engreitz/Users/kangh/tutorials/SingleCellAnalysisTutorial/scRNAseqAnalysis/",help="Project Directory")
)

opt <- parse_args(OptionParser(option_list=option.list))


##########################################################################################
## directories and constants
PROJECT=opt$project
DATADIR=opt$datadir
OUTDIR=paste0(PROJECT, "/outputs/")
FIGDIR= paste0(PROJECT, "/figures/")
SAMPLE=opt$sampleName
# create dir if not already
check.dir <- c(OUTDIR, FIGDIR)
invisible(lapply(check.dir, function(x) { if(!dir.exists(x)) dir.create(x, recursive=T) }))


##########################################################################################
## load data

## Read 10X Matrix
my.data <- Read10X(data.dir = paste0(DATADIR, SAMPLE, "/outs/filtered_feature_bc_matrix/")) # two matrices: "Gene Expression", "Peaks" ## takes about a minute
gex.mtx <- my.data[["Gene Expression"]]
atac.mtx <- my.data[["Peaks"]]


##########################################################################################
## Prepare the data
## Create Seurat Object
s <- CreateSeuratObject(
    counts = gex.mtx,
    project = SAMPLE
)

## Choose filters, use cells with Good_singlet, and filter MT/RP
s[["percent.mt"]] <- PercentageFeatureSet(s, pattern = "^MT-")
s[["percent.ribo"]] <- PercentageFeatureSet(s, pattern = "^RPS|^RPL")

## UMAP on gex and atac separately
s <- SCTransform(s) ## takes about two minutes (~7,000 cells)
s <- RunPCA(s, verbose = FALSE)
## need to determine the dimensionality of the data (recommend >10)
ElbowPlot(s)
selectedDim <- 10
selectedResolution <- 0.06
s <- FindNeighbors(s, dims = 1:selectedDim)
s <- FindClusters(s, resolution = selectedResolution) ## The 'resolution' parameter controls how many UMAP clusters there are, but it does not change the shape of the cluster
s <- RunUMAP(s, dims = 1:selectedDim)

## ## save Seurat Object
## saveRDS(s, file=paste0(OUTDIR, SAMPLE,"_SeuratObjectWithUMAP.RDS"))


##########################################################################################
## Plots

## Visualize UMAP
pdf(paste0(FIGDIR, SAMPLE, "_PCADim", selectedDim, "_Resolution", selectedResolution, "_UMAP.pdf"))
DimPlot(s, reduction = "umap", label=TRUE) %>% print()
dev.off()

## QC metrics
FeaturePlot(s, reduction = "umap", features='nCount_RNA') %>% print() ## UMI / cell
FeaturePlot(s, reduction = "umap", features='nFeature_RNA') %>% print() ## number of genes detected / cell
FeaturePlot(s, reduction = "umap", features='percent.mt') %>% print() ## % mitochondrial gene among all genes detected
FeaturePlot(s, reduction = "umap", features='percent.ribo') %>% print() ## % ribosomal gene among all genes detected
QC.metrics <- c("nCount_RNA", "nFeature_RNA", "percent.mt", "percent.ribo")
FeaturePlot(s, reduction = "umap", features = QC.metrics, ncol = 2) %>%

## other ways of viewing QC metrics
## Violin plot
VlnPlot(s, features = 'nCount_RNA', pt.size=0) %>% print()

## Visualize Gene Expression
gene <- "NOS3"
FeaturePlot(s, reduction = "umap", features=gene) %>% print()
gene.set <- c("NOS3", "PECAM1", "SOX2", "NANOG")
FeaturePlot(s, reduction = "umap", features=gene.set, ncol = 2) %>% print()
RidgePlot(s, features = gene.set, ncol = 2) %>% print()


##########################################################################################
## Further process the data
maxMt <- opt$maxMt ## % mitochondrial gene filter
maxCount <- opt$maxCount ## maximum UMI / cell filter
minUniqueGenes <- opt$minUniqueGenes ## minimum number of unique gene / cell filter

## create another Seurat Object that has the filtered cells
sf <- subset(s, subset = percent.mt < maxMt & nCount_RNA < maxCount & nFeature_RNA > minUniqueGenes)

## save the filtered data
saveRDS(sf, file=paste0(OUTDIR, SAMPLE,"_filteredSeuratObjectWithUMAP.RDS"))


## Differential Expression Analysis
DEAnalysis.df <- FindMarkers(sf, ident.1=c(0,3), ident.2=c(1,2)) %>% mutate(Gene=rownames(.))