Single cell differential composition analysis

using bootstrapping and generalized linear mixed model

1 Data preprocessing

1.1 Load BgrB Seurat object

Code

# load in the genome annotation
# replace "_" with "-" as Seurat cannot recognize "_"
annotations_mz <- readGFFAsGRanges('/Users/haowenhe/Desktop/GCF_000238955.4_M_zebra_UMD2a_genomic.gtf')
annotations_mz$gene_id <- str_replace_all(annotations_mz$gene_id, "_", "-")
annotations_mz$gene_name <- annotations_mz$gene_id
annotations_mz$gene <- str_replace_all(annotations_mz$gene, "_", "-")
annotations_mz$transcript_id <- str_replace_all(annotations_mz$transcript_id, "_", "-")
annotations_mz$tx_id <- annotations_mz$transcript_id
# export(annotations_mz, 'annotations_mz.gtf')

# load in the gene annotation
mz_genome_file <- "/users/haowenhe/Desktop/GCF_000238955.4_M_zebra_UMD2a_genomic.fna"
mz_genome <- read_genome(format = 'fasta', file = mz_genome_file)

# remove the long text description after the feature names in annotation
my_mz_genome <- mz_genome
names(my_mz_genome) <- reshape2::colsplit(names(my_mz_genome), " ", c("1", '2'))[,1]

# load tel&ob data
tel_final <- readRDS("/users/haowenhe/Downloads/tel.ob_integrated_atac.gex_FINAL_04.17.24_KL.rds")

1.2 Update file path for Fragment objects

Code

DefaultAssay(tel_final) <- "peaks"

# get list of fragment objects
frags <- Fragments(tel_final)  
# remove fragment information from assay
Fragments(tel_final) <- NULL  

GetFragmentData(object = frags[[1]], slot = "path")

[1] "/Users/kathrynleatherbury/Dropbox (GaTech)/BioSci-Streelman/Katie/BgrB_atac_gex_multiome/CellRanger/CellRanger_MT_final/tob9_MTcontig/atac_fragments.tsv.gz"

Code

newpath_9 <- "/users/haowenhe/Desktop/BgrB_fragments/Pool9/atac_fragments.tsv.gz"
frags[[1]] <- UpdatePath(frags[[1]], new.path = newpath_9)  

Computing hash

Code

GetFragmentData(object = frags[[2]], slot = "path")

[1] "/Users/kathrynleatherbury/Dropbox (GaTech)/BioSci-Streelman/Katie/BgrB_atac_gex_multiome/CellRanger/CellRanger_MT_final/tob13_MTcontig/atac_fragments.tsv.gz"

Code

newpath_13 <- "/users/haowenhe/Desktop/BgrB_fragments/Pool13/atac_fragments.tsv.gz"
frags[[2]] <- UpdatePath(frags[[2]], new.path = newpath_13)

Computing hash

Code

GetFragmentData(object = frags[[3]], slot = "path")

[1] "/Users/kathrynleatherbury/Dropbox (GaTech)/BioSci-Streelman/Katie/BgrB_atac_gex_multiome/CellRanger/CellRanger_MT_final/tob10_MTcontig/atac_fragments.tsv.gz"

Code

newpath_10 <- "/users/haowenhe/Desktop/BgrB_fragments/Pool10/atac_fragments.tsv.gz"
frags[[3]] <- UpdatePath(frags[[3]], new.path = newpath_10)

Computing hash

Code

GetFragmentData(object = frags[[4]], slot = "path")

[1] "/Users/kathrynleatherbury/Dropbox (GaTech)/BioSci-Streelman/Katie/BgrB_atac_gex_multiome/CellRanger/CellRanger_MT_final/tob14_MTcontig/atac_fragments.tsv.gz"

Code

newpath_14 <- "/users/haowenhe/Desktop/BgrB_fragments/Pool14/atac_fragments.tsv.gz"
frags[[4]] <- UpdatePath(frags[[4]], new.path = newpath_14)

Computing hash

Code

GetFragmentData(object = frags[[5]], slot = "path")

[1] "/Users/kathrynleatherbury/Dropbox (GaTech)/BioSci-Streelman/Katie/BgrB_atac_gex_multiome/CellRanger/CellRanger_MT_final/tob11_MTcontig/atac_fragments.tsv.gz"

Code

newpath_11 <- "/users/haowenhe/Desktop/BgrB_fragments/Pool11/atac_fragments.tsv.gz"
frags[[5]] <- UpdatePath(frags[[5]], new.path = newpath_11)

Computing hash

Code

GetFragmentData(object = frags[[6]], slot = "path")

[1] "/Users/kathrynleatherbury/Dropbox (GaTech)/BioSci-Streelman/Katie/BgrB_atac_gex_multiome/CellRanger/CellRanger_MT_final/tob12_MTcontig/atac_fragments.tsv.gz"

Code

newpath_12 <- "/users/haowenhe/Desktop/BgrB_fragments/Pool12/atac_fragments.tsv.gz"
frags[[6]] <- UpdatePath(frags[[6]], new.path = newpath_12)

Computing hash

Code

GetFragmentData(object = frags[[7]], slot = "path")

[1] "/Users/kathrynleatherbury/Dropbox (GaTech)/BioSci-Streelman/Katie/BgrB_atac_gex_multiome/CellRanger/CellRanger_MT_final/tob15_MTcontig/atac_fragments.tsv.gz"

Code

newpath_15 <- "/users/haowenhe/Desktop/BgrB_fragments/Pool15/atac_fragments.tsv.gz"
frags[[7]] <- UpdatePath(frags[[7]], new.path = newpath_15)

Computing hash

Code

GetFragmentData(object = frags[[8]], slot = "path")

[1] "/Users/kathrynleatherbury/Dropbox (GaTech)/BioSci-Streelman/Katie/BgrB_atac_gex_multiome/CellRanger/CellRanger_MT_final/tob16_MTcontig/atac_fragments.tsv.gz"

Code

newpath_16 <- "/users/haowenhe/Desktop/BgrB_fragments/Pool16/atac_fragments.tsv.gz"
frags[[8]] <- UpdatePath(frags[[8]], new.path = newpath_16)

Computing hash

Code

# assign update list of fragment objects back to the assay
Fragments(tel_final) <- frags  

1.3 Prepare sex- and condition-specifc data

Code

# remove peaks assay to reduce file size
DefaultAssay(tel_final) <- "RNA"
tel_final[['peaks']] <- NULL

# subset sex-specific data
Idents(tel_final) <- tel_final$SEX
male <- subset(x = tel_final, idents = "MALE")
female <- subset(x = tel_final, idents = "FEMALE")

### prepare data for female-specific analysis ###
exprsMat_female <- female[["RNA"]]$counts
subject_female <- female$Individual
cellTypes_female <- female$primary_celltypes
cond_female <- female$COND

# check input data for female-specific proportion analysis
knitr::kable(table(subject_female, cellTypes_female)) %>%
  kable_material(c("striped", "hover"))

	1_Astro/MG	2_OPC/Oligo	3_Peri	4_In	5_In	6_In	7_In	8_Ex	9_Ex	10_Ex	11_Ex	12_Ex	13_Ex	14_Ex	15_In/Ex
DF1B	77	46	2	147	26	35	0	418	183	88	91	27	18	11	59
DF1C	212	69	15	409	96	73	6	720	364	204	132	34	20	6	195
DF2B	76	29	2	161	1	26	0	332	117	73	52	22	8	2	74
DF2C	72	37	5	261	24	35	0	449	185	108	88	55	15	6	105
DF3B	373	73	19	319	103	55	3	613	280	132	117	58	14	13	212
DF3C	235	69	14	406	102	68	6	773	388	226	151	72	18	15	209
DF5B	307	89	17	360	104	71	1	684	362	182	158	47	17	11	128
DF5C	273	64	14	355	2	60	0	688	298	188	126	54	17	11	160
DF6B	111	32	15	168	50	28	4	257	131	70	57	22	10	2	91
DF6C	89	61	6	339	83	68	4	613	254	173	126	58	8	3	182
DF7B	72	15	1	103	24	19	4	209	103	61	41	12	2	4	78
DF7C	62	37	5	262	7	43	1	441	192	136	102	29	11	3	93
DF8B	181	69	12	312	4	71	0	697	240	162	122	54	19	16	178
DF8C	442	88	10	355	101	53	2	562	281	171	114	37	18	13	182
DF9B	84	14	3	138	34	24	2	204	125	54	48	17	2	2	60
DF9C	74	33	6	216	1	36	0	437	193	107	68	40	13	2	104

Code

knitr::kable(table(cellTypes_female, cond_female)) %>%
  kable_material(c("striped", "hover"))

	BHVE	CTRL
1_Astro/MG	1281	1459
2_OPC/Oligo	367	458
3_Peri	71	75
4_In	1708	2603
5_In	346	416
6_In	329	436
7_In	14	19
8_Ex	3414	4683
9_Ex	1541	2155
10_Ex	822	1313
11_Ex	686	907
12_Ex	259	379
13_Ex	90	120
14_Ex	61	59
15_In/Ex	880	1230

Code

clusters <- female$primary_celltypes
Idents(female) <- female$Individual
sample <- Idents(female)
prop.list <- getTransformedProps(clusters, sample)
Proportions <- as.vector(t(prop.list$Proportions))
Samples <- rep(colnames(prop.list$Proportions), nrow(prop.list$Proportions))
Clusters <- rep(rownames(prop.list$Proportions),each=ncol(prop.list$Proportions))
plotdf <- data.frame(Cluster.Name=Samples, Clusters=Clusters,Proportions=Proportions)
plotdf$Cond <- ifelse(substr(plotdf$Cluster.Name, nchar(plotdf$Cluster.Name), nchar(plotdf$Cluster.Name)) == "B", "BHVE",
                        ifelse(substr(plotdf$Cluster.Name, nchar(plotdf$Cluster.Name), nchar(plotdf$Cluster.Name)) == "C", "CTRL", NA))
plotdf$Clusters <- factor(plotdf$Clusters, levels = c("1_Astro/MG", "2_OPC/Oligo", "3_Peri", "4_In",
                                                      "5_In", "6_In", "7_In",
                                                      "8_Ex", "9_Ex", "10_Ex", "11_Ex", "12_Ex", 
                                                      "13_Ex", "14_Ex", "15_In/Ex"))

ggplot(plotdf, aes(x=as.factor(Cluster.Name),y=Proportions,fill=Clusters)) +
    geom_bar(stat="identity") +
    theme(axis.text.x = element_text(size=12),
          axis.text.y = element_text(size=12),
          axis.title = element_text(size=14),
          legend.text = element_text(size=12),
          legend.title = element_text(size=14)) +
  theme_bw() +
  xlab("Individual") +
  theme(axis.text.x = ggtext::element_markdown(size = 9))+ 
          ggtitle("Proportions of primary cell types by individual") + 
  theme(plot.title = element_text(size = 18, hjust = 0)) + 
  theme(panel.grid.major = element_blank(),
                        panel.grid.minor = element_blank(),
          legend.position = "bottom") +
    guides(fill = guide_legend(ncol = 10)) +
  facet_grid(. ~ Cond, scales = "free", space = 'free') +
  theme(strip.background =element_rect(fill='white'))+
  theme(strip.text = element_text(colour = 'black'))

Code

res_scDC_noClust <- scDC_noClustering(cellTypes_female, subject_female, calCI = TRUE, 
                                     calCI_method = c("percentile", "BCa", "multinom"),
                                     nboot = 100)

[1] "Calculating sample proportion..."
[1] "Calculating bootstrap proportion..."
[1] "Calculating percentile ..." "Calculating BCa ..."       
[3] "Calculating multinom ..."  
[1] "Calculating z0 ..."
[1] "Calculating acc ..."

Code

barplotCI(res_scDC_noClust, c("BHVE","CTRL","BHVE","CTRL",
                              "BHVE","CTRL","BHVE","CTRL",
                              "BHVE","CTRL","BHVE","CTRL",
                              "BHVE","CTRL","BHVE","CTRL"))

Code

densityCI(res_scDC_noClust, c("BHVE","CTRL","BHVE","CTRL",
                              "BHVE","CTRL","BHVE","CTRL",
                              "BHVE","CTRL","BHVE","CTRL",
                              "BHVE","CTRL","BHVE","CTRL"))

Code

# res_GLM <- fitGLM(res_scDC_noClust, c("BHVE","CTRL","BHVE","CTRL",
#                               "BHVE","CTRL","BHVE","CTRL",
#                               "BHVE","CTRL","BHVE","CTRL",
#                               "BHVE","CTRL","BHVE","CTRL"), 
#                   pairwise = TRUE)
# 
# summary(res_GLM$pool_res_fixed)
# summary(res_GLM$pool_res_random)

Code

Idents(female) <- female$primary_celltypes
female_8Glut <- subset(x = female, idents = "8_Ex")

clusters <- female_8Glut$secondary_celltypes
Idents(female_8Glut) <- female_8Glut$Individual
sample <- Idents(female_8Glut)
prop.list <- getTransformedProps(clusters, sample)

Performing logit transformation of proportions

Code

Proportions <- as.vector(t(prop.list$Proportions))
Samples <- rep(colnames(prop.list$Proportions), nrow(prop.list$Proportions))
Clusters <- rep(rownames(prop.list$Proportions),each=ncol(prop.list$Proportions))
plotdf <- data.frame(Cluster.Name=Samples, Clusters=Clusters,Proportions=Proportions)
plotdf$Cond <- ifelse(substr(plotdf$Cluster.Name, nchar(plotdf$Cluster.Name), nchar(plotdf$Cluster.Name)) == "B", "BHVE",
                        ifelse(substr(plotdf$Cluster.Name, nchar(plotdf$Cluster.Name), nchar(plotdf$Cluster.Name)) == "C", "CTRL", NA))
plotdf <- plotdf[plotdf$Proportions != 0,]

# install.packages("ggrain")
library(ggrain)

Registered S3 methods overwritten by 'ggpp':
  method                  from   
  heightDetails.titleGrob ggplot2
  widthDetails.titleGrob  ggplot2

Code

# remotes::install_github('jorvlan/raincloudplots')
library(raincloudplots)

#Define colors
mentored_col = '#0B8AAD'
control_col = "#21731B"
#Mentored Group Colors
m_pre = "#52D1F4"
m_during = '#0B8AAD'
m_post = '#0B8AAD'
#Control Group Colors
c_pre = "darkcyan"
c_during = "#21731B"
c_post = "gold"


BHVE <- plotdf %>% dplyr::filter(Clusters == "8.1_Glut" & Cond == "BHVE")
CTRL <- plotdf %>% dplyr::filter(Clusters == "8.1_Glut" & Cond == "CTRL")

BHVE$Cluster.Name <- factor(BHVE$Cluster.Name, levels = c("DF1B", "DF2B", "DF3B", "DF5B", "DF6B", "DF7B", "DF8B", "DF9B"))
CTRL$Cluster.Name <- factor(CTRL$Cluster.Name, levels = c("DF1C", "DF2C", "DF3C", "DF5C", "DF6C", "DF7C", "DF8C", "DF9C"))

BHVE <- BHVE[order(BHVE$Cluster.Name),]
CTRL <- CTRL[order(CTRL$Cluster.Name),]
# Create Mentored Pre-Post 1x1
Mentored_pre_post <- data_1x1( 
  array_1 = CTRL$Proportions, #first set of values
  array_2 = BHVE$Proportions, #second set of values
  jit_distance = .09,
  jit_seed = 321) 

# Mentored Pre-Post 1x1 raincloud plot - ties
M_pre_post_ties <- raincloud_1x1_repmes(
  data = Mentored_pre_post,
  colors = (c(c_pre, c_post)), 
  fills = (c(c_pre, c_post)), 
  line_color = 'gray',
  line_alpha = .3,
  size = 1.5,
  alpha = .5,
  align_clouds = FALSE) +
  scale_x_continuous(breaks=c(1,2), labels=c("CTRL", "BHVE"), limits=c(0, 3)) +
  xlab("") + 
  ylab("8.1_Glut proportion") +
  labs(title = "8.1_Glut proportion within 8_Glut (Female)") +
  theme_classic()
M_pre_post_ties

Code

Idents(male) <- male$primary_celltypes
male_8Glut <- subset(x = male, idents = "8_Ex")

clusters <- male_8Glut$secondary_celltypes
Idents(male_8Glut) <- male_8Glut$Individual
sample <- Idents(male_8Glut)
prop.list <- getTransformedProps(clusters, sample)

Performing logit transformation of proportions

Code

Proportions <- as.vector(t(prop.list$Proportions))
Samples <- rep(colnames(prop.list$Proportions), nrow(prop.list$Proportions))
Clusters <- rep(rownames(prop.list$Proportions),each=ncol(prop.list$Proportions))
plotdf <- data.frame(Cluster.Name=Samples, Clusters=Clusters,Proportions=Proportions)
plotdf$Cond <- ifelse(substr(plotdf$Cluster.Name, nchar(plotdf$Cluster.Name), nchar(plotdf$Cluster.Name)) == "B", "BHVE",
                        ifelse(substr(plotdf$Cluster.Name, nchar(plotdf$Cluster.Name), nchar(plotdf$Cluster.Name)) == "C", "CTRL", NA))
plotdf <- plotdf[plotdf$Proportions != 0,]

BHVEm <- plotdf %>% dplyr::filter(Clusters == "8.1_Glut" & Cond == "BHVE")
CTRLm <- plotdf %>% dplyr::filter(Clusters == "8.1_Glut" & Cond == "CTRL")

BHVEm$Cluster.Name <- factor(BHVEm$Cluster.Name, levels = c("DM1B", "DM2B", "DM3B", "DM5B", "DM6B", "DM7B", "DM8B", "DM9B"))
CTRLm$Cluster.Name <- factor(CTRLm$Cluster.Name, levels = c("DM1C", "DM2C", "DM3C", "DM5C", "DM6C", "DM7C", "DM8C", "DM9C"))

BHVEm <- BHVEm[order(BHVEm$Cluster.Name),]
CTRLm <- CTRLm[order(CTRLm$Cluster.Name),]
# Create Mentored Pre-Post 1x1
Mentored_pre_post <- data_1x1( 
  array_1 = CTRLm$Proportions, #first set of values
  array_2 = BHVEm$Proportions, #second set of values
  jit_distance = .09,
  jit_seed = 321) 

# Mentored Pre-Post 1x1 raincloud plot - ties
M_pre_post_ties <- raincloud_1x1_repmes(
  data = Mentored_pre_post,
  colors = (c(c_pre, c_post)), 
  fills = (c(c_pre, c_post)), 
  line_color = 'gray',
  line_alpha = .3,
  size = 1.5,
  alpha = .5,
  align_clouds = FALSE) +
  scale_x_continuous(breaks=c(1,2), labels=c("CTRL", "BHVE"), limits=c(0, 3)) +
  xlab("") + 
  ylab("8.1_Glut proportion") +
  labs(title = "8.1_Glut proportion within 8_Glut (Male)") +
  theme_classic()
M_pre_post_ties

Code

m_pre = "#D7C6CA"
m_post = '#713948'
c_pre = "#C9CBD6"
c_post = "#495373"

#Create 2x2 together dataframe
Both_pre_post_together <- data_2x2(
  array_1 = CTRLm$Proportions, #first set of values - group 1
  array_2 = BHVEm$Proportions, #first set of values - group 2
  array_3 = CTRL$Proportions, #second set of values - group 1
  array_4 = BHVE$Proportions, #second set of values - group 2
  labels = (c('CTRL','BHVE')),
  jit_distance = .09,
  jit_seed = 321,
  spread_x_ticks = FALSE) #We'll be using this data in our not-spread chart

# 2x2 raincloud plot - together
pre_post_together <- raincloud_2x2_repmes(
  data = Both_pre_post_together,
  colors = c(m_pre, m_post, c_pre, c_post),
  fills = c(m_pre, m_post, c_pre, c_post),
  size = 1.5,
  alpha = .5,
  spread_x_ticks = FALSE) +
  scale_x_continuous(breaks=c(1,2), labels=c("Male", "Female"), limits=c(0, 3)) +
  xlab("") + 
  ylab("8.1_Glut proportion") +
  labs(title = "8.1_Glut proportion within 8_Glut (Male & Female)") +
  theme_classic() 
pre_post_together