Example analysis

Author

Roujin An

Published

November 11, 2024

1 Question to be answered

  1. Which genes are most differentially expressed between specific cell clusters, and how do these genes inform the functional differences between these cell types?
  2. What biological pathways are enriched in differentially expressed genes, and how might these pathways be relevant to disease mechanisms?
  3. How do specific clusters (e.g., vascular smooth muscle cells) differ in their expression of key marker genes compared to other clusters?

2 Intended Audiences

The intended audience for this analysis includes academic researchers and clinicians focusing on vascular biology and human cardiovascular health. This analysis is particularly relevant for those studying cellular heterogeneity and gene expression patterns in normal human aorta tissue, providing a baseline for comparison with diseased states such as atherosclerosis. Additionally, bioinformaticians and data scientists interested in single-cell RNA-seq methodologies can benefit from the data wrangling and visualization techniques applied here. Finally, biotech and pharmaceutical industry professionals may find value in understanding baseline gene expression profiles to guide early-stage research in cardiovascular therapeutics and biomarker discovery.

3 Original data

GSE216860 is a publicly available dataset from the Gene Expression Omnibus (GEO) that provides single-cell RNA sequencing (scRNA-seq) data of normal human ascending aortic tissues from donors of varying ages.

“The study aimed to investigate age-dependent changes in cellular composition, phenotypes, and cell-cell communication within the aorta. The dataset includes samples from six organ donors, categorized into a young group (ages 3 months, 1 year, 3 years, and 12 years) and an old group (ages 57 and 58 years). Analysis of 65,470 cells identified ten cell types, including endothelial cells, smooth muscle cells, fibroblasts, mesenchymal cells, macrophages, and various immune cells.”

link to the dataset: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE216860(Zhang et al. 2022)

Here we work with the pre-processed DEG dataset after data QC, normalization, dimension reduction and DEG analysis (FindAllMarkers) using Seurat pipeline(Satija et al. 2015). Given the original Seurat object is very large to load, we are directly loading the DEG marker genes .csv file.

4 Data Dictionary

Column Description
gene Gene name
avg_log2FC Average log2 fold change between clusters
pct.1 Percentage of cells expressing the gene in Cluster 1/percentage of cells expressing the gene in the testing cluster (markers)
pct.2 Percentage of cells expressing the gene in Cluster 7/percentage of cells expressing the gene in all the other clusters (markers)
p_val_adj Adjusted p-value for differential expression
cell_type Type of cell where the gene is expressed
p_val P-value for differential expression
cluster Cluster in which the gene is analyzed

5 Data loading

Code 
# Loading packages
library(dplyr)

Attaching package: 'dplyr'
The following objects are masked from 'package:stats':

    filter, lag
The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union
Code 
library(tidyverse)
── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
✔ forcats   1.0.0     ✔ readr     2.1.5
✔ ggplot2   3.5.1     ✔ stringr   1.5.1
✔ lubridate 1.9.3     ✔ tibble    3.2.1
✔ purrr     1.0.2     ✔ tidyr     1.3.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
Code 
library(patchwork)
library(openxlsx)
library(ggplot2)
library(stringr)
library(clusterProfiler)

clusterProfiler v4.12.6 Learn more at https://yulab-smu.top/contribution-knowledge-mining/

Please cite:

Guangchuang Yu, Li-Gen Wang, Yanyan Han and Qing-Yu He.
clusterProfiler: an R package for comparing biological themes among
gene clusters. OMICS: A Journal of Integrative Biology. 2012,
16(5):284-287

Attaching package: 'clusterProfiler'

The following object is masked from 'package:purrr':

    simplify

The following object is masked from 'package:stats':

    filter
Code 
library(org.Hs.eg.db)
Loading required package: AnnotationDbi
Loading required package: stats4
Loading required package: BiocGenerics

Attaching package: 'BiocGenerics'

The following objects are masked from 'package:lubridate':

    intersect, setdiff, union

The following objects are masked from 'package:dplyr':

    combine, intersect, setdiff, union

The following objects are masked from 'package:stats':

    IQR, mad, sd, var, xtabs

The following objects are masked from 'package:base':

    anyDuplicated, aperm, append, as.data.frame, basename, cbind,
    colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
    get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
    match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
    Position, rank, rbind, Reduce, rownames, sapply, setdiff, table,
    tapply, union, unique, unsplit, which.max, which.min

Loading required package: Biobase
Welcome to Bioconductor

    Vignettes contain introductory material; view with
    'browseVignettes()'. To cite Bioconductor, see
    'citation("Biobase")', and for packages 'citation("pkgname")'.

Loading required package: IRanges
Loading required package: S4Vectors

Attaching package: 'S4Vectors'

The following object is masked from 'package:clusterProfiler':

    rename

The following objects are masked from 'package:lubridate':

    second, second<-

The following object is masked from 'package:tidyr':

    expand

The following objects are masked from 'package:dplyr':

    first, rename

The following object is masked from 'package:utils':

    findMatches

The following objects are masked from 'package:base':

    expand.grid, I, unname


Attaching package: 'IRanges'

The following object is masked from 'package:clusterProfiler':

    slice

The following object is masked from 'package:lubridate':

    %within%

The following object is masked from 'package:purrr':

    reduce

The following objects are masked from 'package:dplyr':

    collapse, desc, slice


Attaching package: 'AnnotationDbi'

The following object is masked from 'package:clusterProfiler':

    select

The following object is masked from 'package:dplyr':

    select
Code 
library(enrichplot)
library(DOSE)
DOSE v3.30.5  For help: https://yulab-smu.top/biomedical-knowledge-mining-book/

If you use DOSE in published research, please cite:
Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an R/Bioconductor package for Disease Ontology Semantic and Enrichment analysis. Bioinformatics 2015, 31(4):608-609

The UMAP of the Seurat object after scTransform:

UMAP of the normal human aorta scRNA-seq data
Code 
# loading data
markers <- read.csv('/Users/roujinan/Desktop/StatsComp_project1/DataAnalysis/markers_data.csv')
head(markers)
  p_val avg_log2FC pct.1 pct.2 p_val_adj cluster   gene
1     0   2.711262 0.871 0.178         0       0  SFRP2
2     0   2.898204 0.838 0.169         0       0 SCARA5
3     0   2.488215 0.903 0.285         0       0   MMP2
4     0   2.379625 0.865 0.247         0       0  ABCA8
5     0   2.546620 0.838 0.225         0       0   PODN
6     0   2.728942 0.834 0.222         0       0 ABI3BP

We have another marker genes and cell type mapping dictionary, we will load this data for mapping cluster identity. This is an empirical mapping of markers genes of major aortic cell types. However this .csv is based on mouse gene, we will need to convert all lowercase letter in gene names into uppercase, and change one gene H2-aa to HLA-DRA to completely change the list of mouse markers into human markers.

Code 
# Loading cell type marker genes dictionary
Mapper <- read.csv('/Users/roujinan/Desktop/StatsComp_project1/DataAnalysis/atherosclerotic_plaque_marker_genes.csv')
head(Mapper)
    gene           cell_type
1 Pecam1   Endothelial Cells
2   Cdh5   Endothelial Cells
3    Vwf   Endothelial Cells
4   Nos3   Endothelial Cells
5  Acta2 Smooth Muscle Cells
6  Myh11 Smooth Muscle Cells
Code 
# Use mutate to replace 'H2-aa' with 'HLA-DRA' in the 'gene' column
Mapper<- Mapper %>%
  mutate(gene = ifelse(gene == "H2-aa", "HLA-DRA", gene))

# Convert 'gene' column to uppercase, 1st usage of dplyr function mutate
Mapper <- Mapper %>%
  mutate(gene = toupper(gene))

print(Mapper)
      gene           cell_type
1   PECAM1   Endothelial Cells
2     CDH5   Endothelial Cells
3      VWF   Endothelial Cells
4     NOS3   Endothelial Cells
5    ACTA2 Smooth Muscle Cells
6    MYH11 Smooth Muscle Cells
7     CNN1 Smooth Muscle Cells
8    TAGLN Smooth Muscle Cells
9     CD68         Macrophages
10    CD14         Macrophages
11    MSR1         Macrophages
12   ITGAM         Macrophages
13    APOE          Foam Cells
14   ACAT1          Foam Cells
15    CD36          Foam Cells
16   PLIN2          Foam Cells
17    CD3E             T Cells
18     CD4             T Cells
19    CD8A             T Cells
20    CD28             T Cells
21    CD19             B Cells
22   MS4A1             B Cells
23    CD22             B Cells
24    IGHM             B Cells
25   ITGAX     Dendritic Cells
26    CD83     Dendritic Cells
27   CD209     Dendritic Cells
28 HLA-DRA     Dendritic Cells
29     VIM         Fibroblasts
30     FAP         Fibroblasts
31     FN1         Fibroblasts
32  PDGFRA         Fibroblasts
33  ITGA2B           Platelets
34   ITGB3           Platelets
35   GP1BA           Platelets
36    SELP           Platelets
37 CEACAM6         Neutrophils
38     MPO         Neutrophils
39    FUT4         Neutrophils
40   ELANE         Neutrophils

Data preparation is done. Next we want to select the genes that are present in both dataset, and to map the identity of each cluster.

Code 
# Innerjoin two dataset and arrange the joint output. The use of 3 more dolyr and tidyr function: inner_join, select and arrange.

IJ <- inner_join(markers, Mapper,by= 'gene')

ClusterIdentity<- IJ %>% dplyr::select(c('gene','cell_type','p_val_adj','cluster')) %>% arrange(cluster)

ClusterIdentity
       gene           cell_type     p_val_adj cluster
1    PDGFRA         Fibroblasts  0.000000e+00       0
2       FN1         Fibroblasts  0.000000e+00       0
3     MYH11 Smooth Muscle Cells  0.000000e+00       1
4      CNN1 Smooth Muscle Cells  0.000000e+00       1
5     ACTA2 Smooth Muscle Cells  0.000000e+00       1
6     TAGLN Smooth Muscle Cells  0.000000e+00       1
7       VIM         Fibroblasts  0.000000e+00       1
8     ACAT1          Foam Cells  4.128197e-38       1
9    PDGFRA         Fibroblasts  0.000000e+00       2
10     APOE          Foam Cells  0.000000e+00       2
11     APOE          Foam Cells  0.000000e+00       3
12    ACTA2 Smooth Muscle Cells  0.000000e+00       3
13    TAGLN Smooth Muscle Cells  0.000000e+00       3
14      VIM         Fibroblasts 5.494357e-233       3
15     CD36          Foam Cells 1.152690e-197       3
16     CD14         Macrophages  0.000000e+00       4
17     MSR1         Macrophages  0.000000e+00       4
18  HLA-DRA     Dendritic Cells  0.000000e+00       4
19      CD4             T Cells  0.000000e+00       4
20     CD83     Dendritic Cells  0.000000e+00       4
21    ITGAX     Dendritic Cells  0.000000e+00       4
22     CD36          Foam Cells  0.000000e+00       4
23    ITGAM         Macrophages  0.000000e+00       4
24    PLIN2          Foam Cells 3.164482e-198       4
25      VWF   Endothelial Cells  0.000000e+00       5
26   PECAM1   Endothelial Cells  0.000000e+00       5
27  HLA-DRA     Dendritic Cells  0.000000e+00       5
28     SELP           Platelets  0.000000e+00       5
29     CDH5   Endothelial Cells  0.000000e+00       5
30     CD36          Foam Cells  0.000000e+00       5
31     CD3E             T Cells  0.000000e+00       6
32     CD8A             T Cells  0.000000e+00       6
33     CD83     Dendritic Cells 5.592862e-163       6
34    MYH11 Smooth Muscle Cells  0.000000e+00       7
35      FN1         Fibroblasts  0.000000e+00       7
36     CNN1 Smooth Muscle Cells  0.000000e+00       7
37    ACTA2 Smooth Muscle Cells  0.000000e+00       7
38    TAGLN Smooth Muscle Cells  0.000000e+00       7
39      VIM         Fibroblasts  0.000000e+00       7
40     CD14         Macrophages  0.000000e+00       8
41     MSR1         Macrophages  0.000000e+00       8
42     CD83     Dendritic Cells  0.000000e+00       8
43      CD4             T Cells  0.000000e+00       8
44  HLA-DRA     Dendritic Cells  0.000000e+00       8
45    PLIN2          Foam Cells  0.000000e+00       8
46    CD209     Dendritic Cells  0.000000e+00       8
47    ITGAM         Macrophages  0.000000e+00       8
48     CD68         Macrophages  0.000000e+00       8
49    ITGAX     Dendritic Cells  0.000000e+00       8
50     APOE          Foam Cells 2.820754e-276       8
51     CD3E             T Cells  0.000000e+00       9
52     CD8A             T Cells  0.000000e+00       9
53    MYH11 Smooth Muscle Cells  0.000000e+00      10
54     CNN1 Smooth Muscle Cells  0.000000e+00      10
55    ACTA2 Smooth Muscle Cells  0.000000e+00      10
56    TAGLN Smooth Muscle Cells  0.000000e+00      10
57      VIM         Fibroblasts  0.000000e+00      10
58    ACAT1          Foam Cells  5.187792e-15      10
59   PECAM1   Endothelial Cells  0.000000e+00      11
60     CDH5   Endothelial Cells  0.000000e+00      11
61      VWF   Endothelial Cells  0.000000e+00      11
62     CD36          Foam Cells 1.475665e-143      11
63      FN1         Fibroblasts  3.366928e-78      11
64      VIM         Fibroblasts  3.851649e-62      11
65   PDGFRA         Fibroblasts  0.000000e+00      12
66      VWF   Endothelial Cells  0.000000e+00      13
67   PECAM1   Endothelial Cells  0.000000e+00      13
68  HLA-DRA     Dendritic Cells  0.000000e+00      13
69     SELP           Platelets  0.000000e+00      13
70     CDH5   Endothelial Cells  0.000000e+00      13
71     CD36          Foam Cells  6.298348e-79      13
72     CD14         Macrophages  0.000000e+00      14
73     MSR1         Macrophages  0.000000e+00      14
74     CD83     Dendritic Cells  0.000000e+00      14
75      CD4             T Cells  0.000000e+00      14
76  HLA-DRA     Dendritic Cells  0.000000e+00      14
77    ITGAM         Macrophages  0.000000e+00      14
78    PLIN2          Foam Cells  0.000000e+00      14
79    CD209     Dendritic Cells  0.000000e+00      14
80    ITGAX     Dendritic Cells  0.000000e+00      14
81    MYH11 Smooth Muscle Cells 8.202856e-170      15
82      FN1         Fibroblasts 1.531678e-125      15
83   PECAM1   Endothelial Cells  0.000000e+00      16
84     CDH5   Endothelial Cells  0.000000e+00      16
85      VIM         Fibroblasts 2.193968e-262      16
86      VWF   Endothelial Cells 1.595629e-223      16
87    ACAT1          Foam Cells  1.813367e-32      16
88     CD36          Foam Cells  0.000000e+00      17
89      VWF   Endothelial Cells  0.000000e+00      17
90   PECAM1   Endothelial Cells  0.000000e+00      17
91     CDH5   Endothelial Cells  0.000000e+00      17
92    ITGAM         Macrophages  0.000000e+00      18
93     CD14         Macrophages  0.000000e+00      18
94     CD83     Dendritic Cells  0.000000e+00      18
95    ITGAX     Dendritic Cells  0.000000e+00      18
96  HLA-DRA     Dendritic Cells 4.406305e-278      18
97    PLIN2          Foam Cells 1.762404e-204      18
98      CD4             T Cells 8.497471e-175      18
99    PLIN2          Foam Cells  0.000000e+00      19
100  PDGFRA         Fibroblasts 9.634330e-300      19
101     FN1         Fibroblasts 2.269406e-190      19
102    APOE          Foam Cells  1.641942e-91      19
103     FN1         Fibroblasts 6.222499e-129      20
104     VIM         Fibroblasts  1.991899e-26      20
105     VIM         Fibroblasts 7.041632e-222      21
106     VWF   Endothelial Cells  0.000000e+00      22
107  PECAM1   Endothelial Cells  0.000000e+00      22
108    SELP           Platelets  0.000000e+00      22
109    CDH5   Endothelial Cells  0.000000e+00      22
110    NOS3   Endothelial Cells  0.000000e+00      22
111 HLA-DRA     Dendritic Cells 1.520953e-288      22
112    CD36          Foam Cells  1.175719e-88      22
113  PDGFRA         Fibroblasts 1.120461e-277      23
114   PLIN2          Foam Cells 5.042019e-273      23
115    APOE          Foam Cells 1.420083e-161      23
116     VIM         Fibroblasts  2.443916e-61      23
117     VWF   Endothelial Cells  0.000000e+00      24
118  PECAM1   Endothelial Cells  0.000000e+00      24
119 HLA-DRA     Dendritic Cells  0.000000e+00      24
120    SELP           Platelets  0.000000e+00      24
121    CDH5   Endothelial Cells 6.599280e-203      24
122    CD36          Foam Cells  1.413439e-21      24
123   PLIN2          Foam Cells  1.215111e-16      24
124    CNN1 Smooth Muscle Cells  0.000000e+00      25
125   MYH11 Smooth Muscle Cells  0.000000e+00      25
126   ACTA2 Smooth Muscle Cells  0.000000e+00      25
127   TAGLN Smooth Muscle Cells 9.055596e-217      25
128   ACAT1          Foam Cells  1.190212e-44      25
129   PLIN2          Foam Cells  7.636017e-32      25
130    CNN1 Smooth Muscle Cells  0.000000e+00      26
131   MYH11 Smooth Muscle Cells  0.000000e+00      26
132     FN1         Fibroblasts  0.000000e+00      26
133   ACTA2 Smooth Muscle Cells  0.000000e+00      26
134   TAGLN Smooth Muscle Cells 9.374771e-295      26
135     VIM         Fibroblasts 7.151773e-206      26
136   ACAT1          Foam Cells  9.007617e-09      26
137     FN1         Fibroblasts 2.924736e-268      27
138     VIM         Fibroblasts 3.288365e-204      27
139   TAGLN Smooth Muscle Cells  1.181996e-96      27
140  PDGFRA         Fibroblasts  2.584394e-34      27
141    CNN1 Smooth Muscle Cells  2.683319e-31      27
142   ACTA2 Smooth Muscle Cells  2.454618e-26      27
143    APOE          Foam Cells  1.929503e-10      27
144   MYH11 Smooth Muscle Cells  4.001013e-08      27
145   ACAT1          Foam Cells  1.517133e-01      27
146   MYH11 Smooth Muscle Cells 9.290329e-252      28
147    CNN1 Smooth Muscle Cells 2.979034e-205      28
148   ACTA2 Smooth Muscle Cells 4.232966e-200      28
149   TAGLN Smooth Muscle Cells 2.452872e-162      28
150  PDGFRA         Fibroblasts  3.483233e-91      28
151     VIM         Fibroblasts  2.443648e-62      28
152   ACAT1          Foam Cells  2.927823e-06      28
153     FN1         Fibroblasts 7.932934e-150      30
154    APOE          Foam Cells 8.316590e-149      30
155    IGHM             B Cells 3.126366e-191      31
156 CEACAM6         Neutrophils  0.000000e+00      32
157    CD36          Foam Cells  2.484222e-56      32
158   ACAT1          Foam Cells  3.477542e-10      32
159    CD3E             T Cells  5.035155e-66      33
160     CD4             T Cells  6.617962e-25      33
161    MSR1         Macrophages  8.401097e-24      33
162    CD14         Macrophages  6.839104e-21      33
163    CD83     Dendritic Cells  1.614827e-19      33
164 HLA-DRA     Dendritic Cells  6.177279e-07      33
165   PLIN2          Foam Cells  5.388189e-04      33
166   ACAT1          Foam Cells  1.000000e+00      33
167   MS4A1             B Cells  0.000000e+00      34
168    CD22             B Cells  0.000000e+00      34
169    CD19             B Cells  0.000000e+00      34
170    IGHM             B Cells 7.497608e-239      34
171 HLA-DRA     Dendritic Cells 7.394811e-110      34
172    CD83     Dendritic Cells  5.014749e-91      34
173  PECAM1   Endothelial Cells  4.975202e-04      35
174     VWF   Endothelial Cells  3.472254e-96      36
175  PECAM1   Endothelial Cells  8.053721e-58      36
176    SELP           Platelets  4.135159e-55      36
177    CNN1 Smooth Muscle Cells  1.545159e-50      36
178   MYH11 Smooth Muscle Cells  1.059195e-48      36
179   ACTA2 Smooth Muscle Cells  2.045518e-40      36
180   TAGLN Smooth Muscle Cells  7.744618e-36      36
181    CDH5   Endothelial Cells  1.363577e-25      36
182 HLA-DRA     Dendritic Cells  2.637106e-23      36
183     VIM         Fibroblasts  2.608712e-13      36
184    CD36          Foam Cells  7.394107e-02      36
185    CD83     Dendritic Cells  5.555551e-03      37
186   PLIN2          Foam Cells  1.000000e+00      37
187    CD36          Foam Cells  1.735400e-29      38
188   PLIN2          Foam Cells  7.951951e-13      38
189   ACAT1          Foam Cells  4.338791e-10      38
190   ACTA2 Smooth Muscle Cells  7.462093e-36      39
191   TAGLN Smooth Muscle Cells  9.980134e-31      39
192    CNN1 Smooth Muscle Cells  2.027748e-30      39
193   MYH11 Smooth Muscle Cells  7.070574e-25      39
194    APOE          Foam Cells  2.252329e-08      39
195    CD36          Foam Cells  7.756982e-01      39
196     VIM         Fibroblasts  1.000000e+00      39
197    MSR1         Macrophages  1.016540e-77      40
198    CD83     Dendritic Cells  6.783886e-60      40
199    CD14         Macrophages  3.000753e-50      40
200    CD3E             T Cells  6.089443e-37      40
201    APOE          Foam Cells  1.953506e-34      40
202   PLIN2          Foam Cells  6.374785e-30      40
203 HLA-DRA     Dendritic Cells  2.184687e-26      40
204   ITGAX     Dendritic Cells  2.655588e-18      40
205     CD4             T Cells  1.950063e-09      40
206     VIM         Fibroblasts  2.423523e-06      40

Based on the chart, it will be easy to map the cluster identity with the markers. Combining with the original markers data, question 1 can be answered. Since with this results, it was apparent that some clusters, for example cluster 18 contains multiple immune cell type markers such as Macrophages, Dendritic cells, Foam cells and T cells, it may be wise to conduct further subclustering to the original data. We didn’t load the data here, thus we don’t assert cluster identity for these “mixed clusters”.

However, cluster 0 shows clear Fibroblast markers and cluster1 shows clear vascular smooth muscle cells markers expression. At the same time, cluster 7 also seems to be a smooth muscle cell cluster, yet it was quite separated from cluster 1 on the UMAP. Thus we will continue here with these three clusters.

Next, we want to know the differentially expressed genes between the two smooth muscle cell cluster 1 and 7. We ran DEG between the two clusters in seurat on the original data, which is not shown here. The data output is VSMCmarkers.csv, which will be loaded here for further analysis.

Code 
# Loading VSMCmarkers
VSMCmarkers <- read.csv('/Users/roujinan/Desktop/StatsComp_project1/DataAnalysis/VSMCmarkers_data.csv')
head(VSMCmarkers)
  p_val avg_log2FC pct.1 pct.2 p_val_adj       gene
1     0  -7.010129 0.041 0.867         0       VCAN
2     0   7.635419 0.835 0.011         0 AL450405.1
3     0  -5.018332 0.131 0.934         0        OGN
4     0   5.068551 0.887 0.092         0      RERGL
5     0  -5.048445 0.106 0.895         0      MYH10
6     0  -5.720245 0.055 0.842         0       CCN3
Code 
# Filter significant markers with adjusted p-value < 0.05, the application of the 5th dplyr and tidyr function: filter
significant_VSMCmarkers <- VSMCmarkers %>%
  dplyr::filter(p_val_adj < 0.05)

dim(VSMCmarkers)
[1] 12591     6
Code 
dim(significant_VSMCmarkers)
[1] 6335    6
Note

Note: The analysis identified key differentially expressed genes between Cluster 1 and Cluster 7. Both clusters are VSMC clusters. We are interested in exploring the difference between these two clusteres, which is very likely to me batch/sample difference due to different sample age.

6 Visualization of DEG between VSMC clusters with volcano plots

Plots below were created using ggplot2, a powerful data visualization package(Wickham 2016).

Code 
# Replace zero p-values with a very small number to avoid issues in the plot
significant_VSMCmarkers <- significant_VSMCmarkers %>%
  mutate(p_val_adj = ifelse(p_val_adj == 0, 1e-350, p_val_adj))

# Volcano plot
ggplot(significant_VSMCmarkers, aes(x = avg_log2FC, y = -log10(p_val_adj))) +
  geom_point(aes(color = p_val_adj < 0.001 & abs(avg_log2FC) > 3), alpha = 0.5) +
  labs(
    title = "Volcano Plot of Cluster 1 vs Cluster 7",
    subtitle = "Highlighting significant markers (p_val_adj < 0.001 and avg_log2FC > 3)",
    x = "Log2 Fold Change",
    y = "-Log10 Adjusted P-Value",
    caption = "Data source: GSE216860"
  ) +
  theme_minimal() +
  scale_color_manual(values = c("grey", "red"))

Code 
# Grey for non-significant, red for significant

7 Barplot

Code 
# Select top 10 markers by absolute log fold change
top_markers <- markers %>%
  arrange(desc(abs(avg_log2FC))) %>%
  head(10)

# Bar plot
ggplot(top_markers, aes(x = reorder(gene, abs(avg_log2FC)), y = abs(avg_log2FC), fill = gene)) +
  geom_col(show.legend = FALSE) +
  labs(
    title = "Top 10 Markers by Absolute Log2 Fold Change",
    subtitle = "Cluster 1 vs Cluster 7",
    x = "Gene",
    y = "Absolute Log2 Fold Change",
    caption = "Top genes based on |log2FC|"
  ) +
  coord_flip() +
  theme_minimal()

8 Dotplots demonstrating DEG expression between clusters

This plot will visualize the relationship between pct.1 and pct.2 for multiple genes, faceted by each gene.

Code 
# Select a few significant genes for visualization
selected_genes <- significant_VSMCmarkers %>%
  filter(gene %in% c("ACTB",'MYH10', 'IGFBP2'))  # Example gene names

# Faceted scatter plot
ggplot(selected_genes, aes(x = pct.1, y = pct.2)) +
  geom_point(aes(color = abs(avg_log2FC)), size = 4, alpha = 0.7) +  # Color points by log2FC
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey") +  # Equality line
  facet_wrap(~ gene) +  # Facet by gene
  labs(
    title = "Faceted Scatter Plot of % Cells Expressing in Cluster 1 vs Cluster 7",
    subtitle = "Colored by Absolute Log2 Fold Change",
    x = "% Cells Expressing (Cluster 1)",
    y = "% Cells Expressing (Cluster 7)",
    color = "Abs(Log2FC)",
    caption = "Dashed line indicates equal expression percentages"
  ) +
  scale_color_gradient(low = "blue", high = "red") +  # Color gradient
  theme_minimal()

9 Heatmap (geom_tile()) for Expression percentage of genes with top50 fold_change

Code 
# Reshape data for heatmap
all_genes_long <- significant_VSMCmarkers %>%
  arrange(desc(abs(avg_log2FC))) %>%
  head(50) %>%
  pivot_longer(cols = c(pct.1, pct.2), names_to = "cluster", values_to = "expression_pct")

# Heatmap for all genes
ggplot(all_genes_long, aes(x = cluster, y = gene, fill = expression_pct)) +
  geom_tile(color = "white") +  # Heatmap tiles
  scale_fill_gradient(low = "blue", high = "red") +  # Gradient color scale
  labs(
    title = "Heatmap of Expression Percentages for top50 Genes(fold change)",
    subtitle = "Comparison of Cluster 1(pct.1) and Cluster 7(pct.2)",
    x = "Cluster",
    y = "Gene",
    fill = "% Cells Expressing",
    caption = 'Heatmap shows percentage of cells expressing each gene per cluster'
  ) +
  theme_minimal() +
  theme(
    axis.text.y = element_text(size = 5),  # Adjust for many genes
    axis.text.x = element_text(angle = 45, hjust = 1)
  )

Warning

Warning: The percentages shown here for pct.1 and pct.2 are based on seurat FindMarkers function, pct.1 and pct.2 are separately set in the function argument. Do not assume the identities of two clusters here. The pct.1 and pct.2 from the markers data was generated with FindAllMarkers function in seurat, where pct.1 is the cluster of interests, pct.2 is all the other clusters.

10 Enrichment analysis

Code 
deg_genes <- significant_VSMCmarkers$gene
# Convert gene symbols to Entrez IDs
entrez_ids <- bitr(deg_genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb = org.Hs.eg.db)
'select()' returned 1:1 mapping between keys and columns
Warning in bitr(deg_genes, fromType = "SYMBOL", toType = "ENTREZID", OrgDb =
org.Hs.eg.db): 7.75% of input gene IDs are fail to map...
Code 
go_enrichment <- enrichGO(
  gene = entrez_ids$ENTREZID,
  OrgDb = org.Hs.eg.db,
  keyType = "ENTREZID",
  ont = "BP",  # Biological Process
  pAdjustMethod = "BH",
  pvalueCutoff = 0.001
)
Code 
dotplot(go_enrichment, showCategory = 20, size = "GeneRatio") + 
  ggtitle("GO Enrichment Analysis") +
  scale_size(range = c(1, 4)) +  # Adjusts the minimum and maximum dot sizes
  theme_minimal() +
  theme(
    plot.title = element_text(size = 12),
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 4),
    legend.text = element_text(size = 8),
    legend.title = element_text(size = 10)
  )
Scale for size is already present.
Adding another scale for size, which will replace the existing scale.

Code 
kegg_enrichment <- enrichKEGG(
  gene = entrez_ids$ENTREZID,
  organism = 'hsa',  # Homo sapiens
  pvalueCutoff = 0.001
)
Reading KEGG annotation online: "https://rest.kegg.jp/link/hsa/pathway"...
Reading KEGG annotation online: "https://rest.kegg.jp/list/pathway/hsa"...
Code 
dotplot(kegg_enrichment, showCategory = 20, size = "GeneRatio") +  # You can use another numeric value column like "Count"
  ggtitle("GO Enrichment Analysis") +
  scale_size(range = c(1, 4)) +  # Adjusts the minimum and maximum dot sizes
  theme_minimal() +
  theme(
    plot.title = element_text(size = 12),
    axis.text.x = element_text(size = 8),
    axis.text.y = element_text(size = 5),
    legend.text = element_text(size = 8),
    legend.title = element_text(size = 10)
  )
Scale for size is already present.
Adding another scale for size, which will replace the existing scale.

Code 
# Convert gene symbols to Entrez IDs
entrez_ids2 <- bitr(
  significant_VSMCmarkers$gene,
  fromType = "SYMBOL",    # Input type (gene symbols)
  toType = "ENTREZID",    # Output type (Entrez IDs)
  OrgDb = org.Hs.eg.db    # Database for human genes
)
'select()' returned 1:1 mapping between keys and columns
Warning in bitr(significant_VSMCmarkers$gene, fromType = "SYMBOL", toType =
"ENTREZID", : 7.75% of input gene IDs are fail to map...
Code 
# Merge Entrez IDs back into your original data
deg_data <- merge(significant_VSMCmarkers, entrez_ids, by.x = "gene", by.y = "SYMBOL")
# View the data with Entrez IDs
ranked_genes <- deg_data %>%
  arrange(desc(avg_log2FC)) %>%
  pull(avg_log2FC)  # Pull log2FoldChange values as vector

# Assign Entrez IDs as names
names(ranked_genes) <- deg_data %>%
  arrange(desc(avg_log2FC)) %>%
  pull(ENTREZID)  # Pull corresponding Entrez IDs

# Verify structure
head(ranked_genes)
    64102     57582       776      3051      4004     84467 
10.663336  8.564719  8.421073  7.670534  7.630421  7.612656
Code 
gsea_go <- gseGO(
  geneList = ranked_genes,
  OrgDb = org.Hs.eg.db,
  keyType = "ENTREZID",
  ont = "BP",
  pvalueCutoff = 0.05
)
using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
preparing geneSet collections...
GSEA analysis...
leading edge analysis...
done...
Code 
gsea_kegg <- gseKEGG(
  geneList = ranked_genes,
  organism = 'hsa',
  pvalueCutoff = 0.05
)
using 'fgsea' for GSEA analysis, please cite Korotkevich et al (2019).
preparing geneSet collections...
GSEA analysis...
leading edge analysis...
done...
Code 
# Dotplot for GO
dotplot(gsea_go, showCategory = 10) +
  ggtitle("GSEA for GO Biological Process") +
  theme_minimal()

Code 
# Dotplot for KEGG
dotplot(gsea_kegg, showCategory = 10) +
  ggtitle("GSEA for KEGG Pathways") +
  theme_minimal()

Key Insight: Single-cell RNA sequencing enables detailed exploration of cellular heterogeneity within complex tissues like the aorta. This analysis is based on Seurat output and ggplot2 to identify and visualize distinct marker genes and defferentially expressed genes across clusters.

11 Summary

This analysis explored the cellular heterogeneity of human aortic tissues using single-cell RNA sequencing data from GSE216860. Differential expression analysis among all UMAP clusters identified cell type identity of each cluster, and DEG analysis between two vascular smooth muscle cell clusters Cluster 1 and Cluster 7 identified key marker genes associated with extracellular matrix and actin filament organization, and cytoskeleton function. Visualization techniques such as volcano plots, bar plots, and heatmaps highlighted significant trends and provided insights into the relative expression and ratio across clusters. These findings demonstrate the utility of computational methods in uncovering potential targets for understanding vascular diseases like atherosclerosis. This study emphasizes the power of integrating data science with biological research to gain actionable insights.

12 Functions Used

Function Purpose Package
mutate() Adds new columns or modifies existing ones. dplyr
filter() Filters rows based on specific conditions. dplyr
arrange() Arranges rows in ascending or descending order based on column(s). dplyr
inner_join() Joins two datasets based on matching keys in both tables. dplyr
select() Selects specific columns from a dataset. dplyr
pivot_longer() Converts wide-format data to long-format (e.g., column to rows). tidyr
geom_point() Creates scatter plots or adds points to plots. ggplot2
geom_bar() Creates bar plots to visualize categorical data. ggplot2
geom_tile() Creates heatmaps by filling tiles with values. ggplot2
facet_wrap() Splits a plot into multiple panels by a factor variable. ggplot2

13 References

Satija, Rahul, Jeffrey A. Farrell, David Gennert, Alexander F. Schier, and Aviv Regev. 2015. “Spatial Reconstruction of Single-Cell Gene Expression Data.” Nature Biotechnology 33 (5): 495–502. https://doi.org/10.1038/nbt.3192.
Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org/.
Zhang, B., K. Zeng, R. Guan, and Y. Yang. 2022. “Single-Cell RNA Sequencing of the Normal Human Ascending Aortic Tissues from Donors of Different Ages Reveals Significant Age-Dependent Changes in the Cellular Composition, Cellular Phenotypes, and Cell-Cell Communication.” https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE216860.