Alelle specific gene expression analysis for single cell RNA sequencing

This vignette uses gene and single nucleotide variant (SNV) expression files generated from the Annotate Various Features for Alignment.

0. Perform cell clustering with Seurat

We begin by performing cell clustering based on gene expression, which is identical to the first step outlined in the Alternative splicing analysis for scRNA-seq. If you’re continuing from the previous vignette, you are free to skip this section.

require(Yano)

Loading required package: Yano

── Attaching packages ────────────────────────────────────────────── Yano 1.0 ──
✔ dplyr   1.1.4     ✔ Seurat  5.2.1
✔ ggplot2 3.5.1

exp <- ReadPISA("./exp/")

obj <- CreateSeuratObject(exp, min.features = 1000, min.cells = 10)
obj[["percent.mt"]] <- PercentageFeatureSet(obj, pattern = "^MT-")
obj <- subset(obj, nFeature_RNA < 9000 & percent.mt < 20)

# Downsampling to 2000 cells for fast testing
obj <- obj[, sample(colnames(obj),2000)]

# To improve visualization, we reduced the resolution when identifying cell clusters.
obj <- NormalizeData(obj) %>% FindVariableFeatures() %>% ScaleData() %>%  RunPCA(verbose=FALSE) %>% FindNeighbors(dims = 1:10, verbose=FALSE) %>% FindClusters(resolution = 0.1, verbose=FALSE) %>% RunUMAP(dims=1:10, verbose=FALSE)

Normalizing layer: counts
Finding variable features for layer counts
Centering and scaling data matrix

DimPlot(obj, label=TRUE, label.size = 5, label.box = TRUE)

1. Perform allele-specific gene expression analysis by using SNV anchors

In this vignette, when a genetic variant is detected, all possible alleles at that location, including both alternative alleles and the reference allele, are annotated and counted. Therefore, genetic variants refer to both alternative alleles and the reference allele in this context.

var <- ReadPISA("./var")
obj[['var']] <- CreateAssayObject(var[,colnames(obj)], min.cells = 10)
# Switch working assay to 'var'
DefaultAssay(obj) <- "var"

obj <- NormalizeData(obj)

# you can skip to read GTF again if you already load it
gtf <- gtf2db("./gencode.v44.annotation.gtf.gz")

[2025-02-20 18:17:15] GTF loading..
[2025-02-20 18:17:27] Load 62700 genes.

# At this moment the meta infomation for genetic variant assay is still empty
obj[['var']][[]] %>% head

data frame with 0 columns and 6 rows

obj <- annoVAR(obj, gtf = gtf)

Parse names ..

# After annotation, the locations are parsed from the feature name, and overlapped gene are also annotated.
obj[['var']][[]] %>% head %>% knitr::kable()

	chr	start	ref	alt	strand	locus	gene_name	type
chr1:631862G>A/+	chr1	631862	G	A	+	chr1:631862/+	MTCO1P12	utr5
chr1:633192G=/+	chr1	633192	G	G	+	chr1:633192/+	MTCO2P12	utr5
chr1:634337T>C/+	chr1	634337	T	C	+	chr1:634337/+	MTATP6P1	utr5
chr1:826352C>T/-	chr1	826352	C	T	-	chr1:826352/-	LINC00115	utr5
chr1:853876T>C/+	chr1	853876	T	C	+	chr1:853876/+	LINC01128	utr5
chr1:944296G>A/-	chr1	944296	G	A	-	chr1:944296/-	NOC2L	utr3

# Select autocorrlated genetic variants 
obj <- RunAutoCorr(obj)

Working on assay : var

Run autocorrelation test for 39228 features.

Runtime : 2.759464 secs

obj <- SetAutoCorrFeatures(obj)

13280 autocorrelated features.

In the alternative splicing analysis, we use mode 1 (the default mode). However, for allele-specific gene expression, we switch to mode 2 to compare the expression of an SNV with other SNVs at the same location. Since no binding assay is specified, the records at the same locus are first aggregated, and then the test feature is subtracted from the aggregated data (mode 2).

obj <- RunBlockCorr(object = obj, bind.name = "locus", mode = 2)

Working on assay var.

Processing 13280 features.

Processing 11274 blocks.

No bind.assay specified, update min.features.per.block to 2.

Use "counts" layer for test features.

Aggregate counts for binding features.

Using 7 threads.

Using method "D" with mode 2.

Use predefined weight matrix "pca_wm".

Runtime : 26.9589 secs.

sel.chrs <- c(1:21, "X")
FbtPlot(obj, val = "locus.padj", remove.chr = TRUE, sel.chrs = sel.chrs)

# We could also zoom in specific region
# Interesting, we can found the peak at chromosome 6 is located at 3' of gene HLA-J
FbtPlot(obj, val = "locus.padj", chr="chr6", start = 30000000, end = 30020000, gtf = gtf)

Zoom in chr6 : 29999000-30021000

# Let's random pick a pair of features at the same locus
obj[['var']][[]] %>% filter(locus.padj < 1e-6) %>% knitr::kable()

	chr	start	ref	alt	strand	locus	gene_name	type	moransi.pval	moransi	autocorr.variable	locus.D	locus.r	locus.mean	locus.var
chr1:87328836C>G/+	chr1	87328836	C	G	+	chr1:87328836/+	LMO4	upstream	0.0e+00	0.3093687	TRUE	0.2969442	0.1340188	0.1327753	0.0206212
chr10:3136689C>T/+	chr10	3136689	C	T	+	chr10:3136689/+	PFKP	utr5	0.0e+00	0.0980815	TRUE	0.3112325	0.1149604	0.1365351	0.0211855
chr10:3136723C>T/+	chr10	3136723	C	T	+	chr10:3136723/+	PFKP	utr5	0.0e+00	0.0834913	TRUE	0.2760167	0.0442716	0.1339272	0.0174969
chr10:12250345A>G/+	chr10	12250345	A	G	+	chr10:12250345/+	CDC123	utr5	0.0e+00	0.0434644	TRUE	0.3402460	-0.0284169	0.1339235	0.0251529
chr10:12250345A=/+	chr10	12250345	A	A	+	chr10:12250345/+	CDC123	utr5	1.6e-06	0.0189449	TRUE	0.3046426	-0.0234226	0.1303622	0.0240605
chr10:17237326C>T/+	chr10	17237326	C	T	+	chr10:17237326/+	VIM	utr5	0.0e+00	0.6865959	TRUE	0.9230209	-0.1926175	0.1331239	0.0246671
chr10:17237326C=/+	chr10	17237326	C	C	+	chr10:17237326/+	VIM	utr5	0.0e+00	0.5822921	TRUE	0.8484525	-0.0404702	0.1364212	0.0263815
chr10:62069780C>T/+	chr10	62069780	C	T	+	chr10:62069780/+	ARID5B	exon	0.0e+00	0.1100910	TRUE	0.4042234	0.0574161	0.1340167	0.0226507
chr10:62069780C=/+	chr10	62069780	C	C	+	chr10:62069780/+	ARID5B	exon	0.0e+00	0.1000460	TRUE	0.4048788	-0.0518440	0.1309692	0.0223058
chr10:71816550C=/-	chr10	71816550	C	C	-	chr10:71816550/-	PSAP	utr3	0.0e+00	0.2918100	TRUE	0.4736578	0.1078816	0.1291249	0.0205535
chr19:5233834T>C/+	chr19	5233834	T	C	+	chr19:5233834/+	PTPRS	antisense_intron	0.0e+00	0.0867366	TRUE	0.3007367	0.1882152	0.1327901	0.0203429
chr2:74969604A>G/+	chr2	74969604	A	G	+	chr2:74969604/+	POLE4	utr5	0.0e+00	0.1260407	TRUE	0.5337651	-0.0247983	0.1319195	0.0192812
chr20:63521953A>C/+	chr20	63521953	A	C	+	chr20:63521953/+	PPDPF	utr5	0.0e+00	0.1082245	TRUE	0.2782363	0.1246119	0.1372135	0.0193050
chr6:30009771C>G/+	chr6	30009771	C	G	+	chr6:30009771/+	ENSG00000290574	utr5	0.0e+00	0.0603756	TRUE	0.3051727	0.0230214	0.1314449	0.0189238
chr6:30009735T>C/+	chr6	30009735	T	C	+	chr6:30009735/+	ENSG00000290574	utr5	0.0e+00	0.0497517	TRUE	0.2930629	-0.0469905	0.1322509	0.0187106
chr6:30009746C>T/+	chr6	30009746	C	T	+	chr6:30009746/+	ENSG00000290574	utr5	0.0e+00	0.0473356	TRUE	0.2759231	-0.0380521	0.1334138	0.0174720
chr6:31268945C>T/-	chr6	31268945	C	T	-	chr6:31268945/-	HLA-C	utr3	0.0e+00	0.2866787	TRUE	0.3188955	-0.1811709	0.1328688	0.0211667
chrM:8860A>G/+	chrM	8860	A	G	+	chrM:8860/+	MT-ATP6	exon	0.0e+00	0.1573206	TRUE	0.3285938	-0.0306096	0.1330598	0.0181256

FeaturePlot(obj, features = c("chr10:17237326C>T/+", "chr10:17237326C=/+", "VIM" ), order=TRUE, pt.size = 1, ncol=3)

Warning: Could not find VIM in the default search locations, found in 'RNA'
assay instead

It appears that the T allele of the gene VIM is expressed in most cell groups, while the C allele seems to be primarily expressed in groups 1 and 4.

Using our spatial dissimilarity analysis method, we have identified many cases of allele-specific gene expression. This phenomenon can be influenced by several known and unknown mechanisms, including gene imprinting, genetic recombination, and somatic variants during cell or disease development. While we won’t explore these mechanisms in detail here, you can refer to our manuscript, which includes several case studies dedicated to this topic.

# Save the Seurat object for further use
saveRDS(obj, file = "glioblastoma5k.rds")

Command(obj)

 [1] "NormalizeData.RNA"        "FindVariableFeatures.RNA"
 [3] "ScaleData.RNA"            "RunPCA.RNA"              
 [5] "FindNeighbors.RNA.pca"    "FindClusters"            
 [7] "RunUMAP.RNA.pca"          "NormalizeData.var"       
 [9] "ParseVarName.var"         "annoVAR.var"             
[11] "RunAutoCorr.var.pca"      "SetAutoCorrFeatures.var" 
[13] "RunBlockCorr.var.pca"