Identifying cell-states

workflowr

• Summary
• Checks
• Past versions

Last updated: 2022-12-16

Checks: 7 0

Knit directory: Multimodal-Plasmacell_manuscript/

This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

---

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

Seed: set.seed(20211005)

The command set.seed(20211005) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: 67dc2d2

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.

The results in this page were generated with repository version 67dc2d2. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.

Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:

Ignored files:
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    output/MOFA_analysis_Donorgroup.hdf5
    Ignored:    output/MOFA_analysis_Donorgroup.rds
    Ignored:    output/MOFA_analysis_Donorgroup_clustered.rds
    Ignored:    output/MOFA_analysis_Donorgroup_noIg.hdf5
    Ignored:    output/MOFA_analysis_Donorgroup_noIg2.hdf5
    Ignored:    output/extra plots.docx
    Ignored:    output/seu.PROT_fix.rds
    Ignored:    output/seu.PROT_live.rds
    Ignored:    output/seu.RNA.rds
    Ignored:    output/seu.fix_norm.rds
    Ignored:    output/seu.fix_norm_cellstate.rds
    Ignored:    output/seu.fix_norm_plasmacells.rds
    Ignored:    output/seu.live_norm.rds
    Ignored:    output/seu.live_norm_cellstate.rds
    Ignored:    output/seu.live_norm_plasmacells.rds
    Ignored:    output/seu.live_norm_plasmacells_RNA.rds
    Ignored:    output/top-PROT-loadings_IgA.tsv
    Ignored:    output/top-PROT-loadings_IgG.tsv
    Ignored:    output/top-PROT-loadings_IgM.tsv
    Ignored:    output/top-gene-loadings_IgA.tsv
    Ignored:    output/top-gene-loadings_IgG.tsv
    Ignored:    output/top-gene-loadings_IgM.csv
    Ignored:    output/top-gene-loadings_IgM.tsv

Untracked files:
    Untracked:  analysis/cellstate_sidetest.Rmd
    Untracked:  analysis/hallmarks.Rmd
    Untracked:  analysis/supplements.Rmd
    Untracked:  code/Import_and_create_seuratobject.R
    Untracked:  data/Seq2Science/
    Untracked:  data/azimuth_PBMCs/
    Untracked:  data/azimuth_bonemarrow/
    Untracked:  data/citeseqcount_htseqcount.zip
    Untracked:  data/genelist.plots.diffmarkers.txt
    Untracked:  data/genelist.plots.diffmarkers2.txt
    Untracked:  data/raw/
    Untracked:  data/supplementary/
    Untracked:  output/paper_figures/

Unstaged changes:
    Modified:   .gitignore
    Modified:   code/load_packages.R

Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.

---

These are the previous versions of the repository in which changes were made to the R Markdown (analysis/cellstate.Rmd) and HTML (docs/cellstate.html) files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view the files as they were in that past version.

File	Version	Author	Date	Message
Rmd	67dc2d2	Jessie van Buggenum	2022-12-16	Fig 1 and 2 and supplementary figures documentation
html	b6ab01b	jessievb	2021-12-06	Build site.
Rmd	791178a	jessievb	2021-12-06	Plasmacell and cell-cycle gating analysis

---

seu.fix <- readRDS( file = "output/seu.fix_norm.rds")
seu.live <- readRDS(file = "output/seu.live_norm.rds")

Summary

• Cell-Cycle scoring
  
• Plasma cell markers (CD27high & IgDlow)
  
• Ig-subtype classification

Cell-cycle state

Cell-cycle scoring was performed scoring algorithm UCell. (For citation see (this preprint)

s.genes <- cc.genes.updated.2019$s.genes
g2m.genes <- cc.genes.updated.2019$g2m.genes

###          --- Live cells ------

# Cell-cycle scoring was performed novel scoring algorithm https://carmonalab.github.io/UCell/UCell_Seurat_vignette.html
# Genesets used also in Seurat Cell-cycle scoring
seu.live <- AddModuleScore_UCell(seu.live, assay = "SCT",slot = "data",features = cc.genes.updated.2019)

## Fetch data for plotting
data.CC.live <- FetchData(seu.live, vars = c('G2M.Score', 'S.Score', "Phase", "donor", "s.genes_UCell", "g2m.genes_UCell"))
data.CC.live$MKI67 <- FetchData(seu.live[["SCT"]], "MKI67", slot = "scale.data")$MKI67

## Plot CC score

p.CCscore.MKI67.live <- ggplot(data = data.CC.live, aes(x = g2m.genes_UCell, y =s.genes_UCell , color = MKI67)) +
  geom_point( position = 'jitter', alpha = 0.8) +
  scale_color_gradientn(colors = c('lightgrey', 'dodgerblue2')) +
  geom_rect(mapping=aes(xmin=-0.01, xmax=0.1, ymin=-0.01, ymax=0.1), color="black", alpha=0)+
  cowplot::theme_cowplot() + 
  labs(title = "RNA-based G2M and S-score mark \ncycling and non-cycling cells", x = "G2M-score", y = "S-score")+
  add.textsize+
  #facet_wrap(.~donor)+
  annotate(geom = "text", label = "Non-\ncycling \n~90%", x=0.01, y = 0.09,hjust = 0,vjust = 1, size = 2) +
  theme(legend.position = c(0.85,0.85), legend.key.size = unit(2, 'mm')) 

## Catagorize cells based on CCscore and MKI67 as cycling or non-cycling. 

data.CC.live <- mutate(data.CC.live, Proliferation_state = ifelse(g2m.genes_UCell >= 0.1 | s.genes_UCell >=0.1 | MKI67 >=1 , "cycling","non-cycling"))

seu.live <- AddMetaData(seu.live, metadata = data.CC.live$Proliferation_state, col.name = "Cellcycle_state")

Percentage_differentiated.live <- data.CC.live %>% dplyr::count(Proliferation_state = factor(Proliferation_state),group = factor(donor)) %>% 
  group_by(group)%>%
    mutate(pct = prop.table(n)) 

## Plot percentage cycling per donor
p.percentage.cyclinglive <- ggplot(Percentage_differentiated.live, aes(x = Proliferation_state , y = pct, fill = group)) + #label = scales::percent(pct,accuracy = 0.1)
    geom_col(position = 'dodge') + 
    scale_y_continuous(labels = scales::percent) +
  cowplot::theme_cowplot() + 
  labs(title="~90% non-cycling cells \nacross three donors",x="", y = "Percentage (%)")+
scale_color_manual("Donor", values=c("#999999", "#E69F00", "#56B4E9"))+
scale_fill_manual("Donor",values=c("#999999", "#E69F00", "#56B4E9"))+
  add.textsize+
  theme(legend.position = c(0.05,0.85), legend.key.size = unit(4, 'mm')) 


###          --- fix cells ------

# Cell-cycle scoring was performed novel scoring algorithm https://carmonalab.github.io/UCell/UCell_Seurat_vignette.html
# Genesets used also in Seurat Cell-cycle scoring
seu.fix <- AddModuleScore_UCell(seu.fix, assay = "SCT",slot = "data",features = cc.genes.updated.2019)

## Fetch data for plotting
data.CC.fix <- FetchData(seu.fix, vars = c('G2M.Score', 'S.Score', "Phase", "donor", "s.genes_UCell", "g2m.genes_UCell"))
data.CC.fix$MKI67 <- FetchData(seu.fix[["SCT"]], "MKI67", slot = "scale.data")$MKI67
data.CC.fix$MKI67.PROT <- FetchData(seu.fix[["PROT"]], "Ki67", slot = "scale.data")$Ki67
data.CC.fix$`p-Rb` <- FetchData(seu.fix[["PROT"]], "p-Rb", slot = "scale.data")$`p-Rb`
data.CC.fix$CyclinA <- FetchData(seu.fix[["PROT"]], "Cyclin A", slot = "scale.data")$`Cyclin`
## Plot CC score

p.CCscore.MKI67.fix <- ggplot(data = data.CC.fix, aes(x = g2m.genes_UCell, y =s.genes_UCell , color = MKI67)) +
  geom_point( position = 'jitter', alpha = 0.8) +
  scale_color_gradientn(colors = c('lightgrey', 'dodgerblue2')) +
  cowplot::theme_cowplot() + 
  labs(title = "RNA-based G2M and S-score mark \ncycling and non-cycling cells", x = "G2M-score", y = "S-score")+
  add.textsize+
  #facet_wrap(.~donor)+
  geom_vline(xintercept = 0.1) +
  geom_hline(yintercept = 0.07) +
  annotate(geom = "text", label = "Non-\ncycling", x=0.01, y = 0.025, size = 2) +
  theme(legend.position = c(0.85,0.85), legend.key.size = unit(2, 'mm')) 

##The advantage of dataset is intracellular phospho-protein markers detected in protein dataset, which can be used to filter filter non-deviding cells. p-Rb and Cyclin A are used:

p.CCscore.pRb <- ggplot(data = data.CC.fix, aes(x = `p-Rb`, y =CyclinA , color = MKI67.PROT)) +
  geom_point( position = 'jitter', alpha = 0.8) +
  scale_color_gradientn(colors = c('lightgrey', 'deeppink3')) +
  cowplot::theme_cowplot() + 
  labs(title = "Intracellular protein markers classify fixed  \ndataset in cycling and non-cycling cells")+
  add.textsize+
  geom_rect(mapping=aes(xmin=-1.8, xmax=1.6, ymin=-2.8, ymax=2.3), color="black", alpha=0)+
  annotate(geom = "text", label = "Non-\ncycling \n>95%", x=-1.7, y = 0.6,hjust = 0,vjust = 1, size = 2)+
  theme(legend.position = c(0.05,0.85), legend.key.size = unit(2, 'mm')) 


## Catagorize cells based on CCscore and MKI67 as cycling or non-cycling. 

data.CC.fix <- mutate(data.CC.fix, Proliferation_state = ifelse( `p-Rb` >= 1.6 | CyclinA >= 2.3, "cycling","non-cycling"))


seu.fix <- AddMetaData(seu.fix, metadata = data.CC.fix$Proliferation_state, col.name = "Cellcycle_state")

Percentage_differentiated.fix <- data.CC.fix %>% dplyr::count(Proliferation_state = factor(Proliferation_state),group = factor(donor)) %>% 
  group_by(group)%>%
    mutate(pct = prop.table(n)) 

## Plot percentage cycling per donor
p.percentage.cycling.fix <- ggplot(Percentage_differentiated.fix, aes(x = Proliferation_state , y = pct, fill = group, label = scales::percent(pct,accuracy = 0.1))) + 
    geom_col(position = 'dodge') + 
    geom_text(position = position_dodge(width = .9),    # move to center of bars
              vjust = -0.5,    # nudge above top of bar
              size = 1.8) + 
    scale_y_continuous(labels = scales::percent) +
  cowplot::theme_cowplot() + 
  labs(title="Fixed dataset classification per donor \n>95% non-cycling cells",x="", y = "Percentage (%)")+
scale_color_manual("Donor", values=c("#E69F00", "#56B4E9"))+
scale_fill_manual("Donor",values=c("#E69F00", "#56B4E9"))+
  add.textsize+
  theme(legend.position = c(0.05,0.85), legend.key.size = unit(2, 'mm')) 

Differentiated Plasma-Blast/Cells

data.CC.fix$CD27 <- FetchData(seu.fix[["PROT"]], "CD27", slot = "scale.data")$CD27
data.CC.fix$IgD <- FetchData(seu.fix[["PROT"]], "IgD", slot = "scale.data")$IgD

data.CC.live$CD27 <- FetchData(seu.live[["PROT"]], "CD27", slot = "scale.data")$CD27
data.CC.live$IgD <- FetchData(seu.live[["PROT"]], "IgD", slot = "scale.data")$IgD

## ------ live
p.scatter.CD27.IgD.live <- ggplot(data = data.CC.live,aes(x = IgD, y = CD27, color = donor)) +
  geom_point( position = 'jitter', alpha = 0.8) +
  geom_rect(mapping=aes(xmin=-2.5, xmax=3, ymin=-2, ymax=2.8), color="black", alpha=0)+
  cowplot::theme_cowplot() + 
  labs(title = "Protein-based gate of CD27+IgD- cells", x = "IgD", y = "CD27")+
  theme(legend.position = "none") +
scale_color_manual("Donor", values=c("#999999", "#E69F00", "#56B4E9"))+
  annotate(geom = "text", label = "CD27+IgD- \n>95%", x=-2.3, y = 2.5,hjust = 0,vjust = 1, size = 2)+
  add.textsize  

data.CC.live <- mutate(data.CC.live, class_switch_plasma = ifelse(IgD <= 3 & CD27 >= -2.5 , "CD27+IgD-","other"))
seu.live <- AddMetaData(seu.live, metadata = data.CC.live$class_switch_plasma, col.name = "CD27_IgD_state")

Percentage_class_switch_plasma.live <- data.CC.live %>% dplyr::count(class_switch_plasma = factor(class_switch_plasma),group = factor(donor)) %>% 
  group_by(group)%>%
    mutate(pct = prop.table(n)) 

p.percentage.class_switch_plasma.live <- ggplot(Percentage_class_switch_plasma.live, aes(x = class_switch_plasma , y = pct, fill = group, label = scales::percent(pct, accuracy = 0.1))) + 
    geom_col(position = 'dodge') + 
    geom_text(position = position_dodge(width = .9),    # move to center of bars
              vjust = -0.5,    # nudge above top of bar
              size = 1.8) + 
    scale_y_continuous(labels = scales::percent) +
  cowplot::theme_cowplot() + 
  labs(title="~95% cells in all three donors \nare Plasma Blast/Cells",x="", y = "Percentage (%)")+
scale_fill_manual(name = "Donor",values=c("#999999", "#E69F00", "#56B4E9")) +
  theme(legend.position = c(0.7,0.85), legend.key.size = unit(2, "mm"))+
  add.textsize

## --- fix
## fix
p.scatter.CD27.IgD.fix <- ggplot(data = data.CC.fix,aes(x = IgD, y = CD27, color = donor)) +
  geom_point( position = 'jitter', alpha = 0.8) +
  geom_rect(mapping=aes(xmin=-4, xmax=4.2, ymin=-1.3, ymax=2), color="black", alpha=0)+
  cowplot::theme_cowplot() + 
  labs(title = "Protein-based gate of CD27+IgD- cells", x = "IgD", y = "CD27")+
theme(legend.position = c(0.85,0.85)) +
scale_color_manual("Donor", values=c( "#E69F00", "#56B4E9"))+
  add.textsize  +
  xlim(c(NA,7.5))

data.CC.fix <- mutate(data.CC.fix, class_switch_plasma = ifelse(IgD <= 4.2 & CD27 >= -1.3 , "CD27+IgD-","other"))
seu.fix <- AddMetaData(seu.fix, metadata = data.CC.fix$class_switch_plasma, col.name = "CD27_IgD_state")

Percentage_class_switch_plasma.fix <- data.CC.fix %>% dplyr::count(class_switch_plasma = factor(class_switch_plasma),group = factor(donor)) %>% 
  group_by(group)%>%
    mutate(pct = prop.table(n)) 

p.percentage.class_switch_plasma.fix <- ggplot(Percentage_class_switch_plasma.fix, aes(x = class_switch_plasma , y = pct, fill = group, label = scales::percent(pct, accuracy = 0.1))) + 
    geom_col(position = 'dodge') + 
    geom_text(position = position_dodge(width = .9),    # move to center of bars
              vjust = -0.5,    # nudge above top of bar
              size = 1.8) + 
    scale_y_continuous(labels = scales::percent) +
  cowplot::theme_cowplot() + 
  labs(title="Fixed dataset % cells CD27+IgD-",x="", y = "Percentage (%)")+
scale_fill_manual(name = "Donor",values=c( "#E69F00", "#56B4E9")) +
  theme(legend.position = c(0.7,0.9), legend.key.size = unit(2, 'mm'))+
  add.textsize

CD27+IgDlow

To determine the cultures ‘differentiation state’ CD27-high & IgD-low cells are gated, representing differentiated plasmablast/cells.

seu.live <- SetIdent(seu.live,value = "CD27_IgD_state")

## RNA differential expression
markers.CD27IgD.RNA <- FindMarkers(seu.live, ident.1 = "CD27+IgD-", ident.2 = "other", assay = "SCT", logfc.threshold = 0.01, test.use = "wilcox", only.pos = T)
markers.CD27IgD.RNA <- filter(markers.CD27IgD.RNA, p_val <= 0.05)

p.CD27IgD.dotplot.genesign <-DotPlot(seu.live,assay = "SCT", features = rev(rownames(markers.CD27IgD.RNA)[1:20]), split.by = "donor", group.by =  "CD27_IgD_state", cols = "RdBu",dot.scale = 3, scale.max = 100, scale.min = 0, scale = T, scale.by = "radius", col.min = -0.5, col.max = 0.5) +
  labs(x = "Differentiation markers - RNA", y = "") +
  cowplot::theme_cowplot() + 
  labs(title=paste0("Top 20 upregulated genes \n(", nrow(markers.CD27IgD.RNA)," total)"),x="Differential expressed genes \n(p-val < 0.05, logfc >= 0.01)", y = "")+
  coord_flip()+
    add.textsize +
  theme(legend.position = "right", legend.key.size = unit(2, 'mm')) +
  scale_y_discrete(labels = c("D1 \n","D2 \nCD27+IgD-", "D3","D1 \n","D2 \nother", "D3")) 
#+  guides(size = "none",color = "none")

p.CD27IgD.Vln.gene.markers <- VlnPlot(seu.live, assay="SCT", features = c("CD27","IRF4", "MZB1", "PTPRC", "MS4A1"), split.by = "donor", group.by =  "CD27_IgD_state", ncol =3, cols = c(colors.donors,colors.donors), combine = T, stack = T, flip = T) +
    add.textsize +
  theme(legend.position = "", legend.key.size = unit(2, 'mm')) + 
  labs(title="CD27+IgD express plasmablast/cell \nRNA markers",x="", y = "Expression")  +
  theme(axis.text.x = element_text(angle = 0, hjust = 0.5))
 
## RNA differential expression
markers.CD27IgD.RNA.neg <- FindMarkers(seu.live, ident.1 = "CD27+IgD-", ident.2 = "other", assay = "SCT", logfc.threshold = -0.01, test.use = "wilcox", only.pos = F)
markers.CD27IgD.RNA.neg <- filter(markers.CD27IgD.RNA.neg, p_val <= 0.05, avg_log2FC < -0.01)

p.CD27IgD.dotplot.genesign.neg <-DotPlot(seu.live,assay = "SCT", features = rev(rownames(markers.CD27IgD.RNA.neg)[1:20]), split.by = "donor", group.by =  "CD27_IgD_state", cols = "RdBu",dot.scale = 3, scale.max = 100, scale.min = 0, scale = T, scale.by = "size", col.min = -0.5, col.max = 0.5) +
  labs(x = "Differentiation markers - RNA", y = "") +
  cowplot::theme_cowplot() + 
  labs(title=paste0("Top 20 downregulated genes \n(", nrow(markers.CD27IgD.RNA.neg)," total)"),x="Differential expressed genes \n(p-val < 0.05, logfc <= -0.01)", y = "")+
  coord_flip()+
    add.textsize +
  theme(legend.position = "right", legend.key.size = unit(2, 'mm')) +
  scale_y_discrete(labels = c("D1 \n","D2 \nCD27+IgD-", "D3","D1 \n","D2 \nother", "D3")) +
    guides(size =  guide_legend(title = "Percent \nexpressed"),color = guide_colorbar(title = "Average scaled\nexpression"))
#+  guides(size = "none",color = "none")

 
## Surface proteins
markers.CD27IgD.PROT <- FindMarkers(seu.live, ident.1 = "CD27+IgD-", ident.2 = "other", assay = "PROT", logfc.threshold = 0.01, test.use = "wilcox", only.pos = T)
markers.CD27IgD.PROT <- filter(markers.CD27IgD.PROT, p_val <= 0.05)

p.CD27IgD.dotplot.PROTsign <-DotPlot(seu.live,assay = "PROT", features = rev(rownames(markers.CD27IgD.PROT)), split.by = "donor", group.by =  "CD27_IgD_state", cols = "RdBu",dot.scale = 3, scale.max = 100, scale.min = 0, scale = T, scale.by = "radius", col.min = -0.5, col.max = 0.5) +
  cowplot::theme_cowplot() + 
  labs(title="Differential upregulated \nSurface protein markers",x="Surface Proteins \n(p-val < 0.005, logfc >= 0.01)", y = "")+
  coord_flip()+
    add.textsize +
  theme(legend.position = "right", legend.key.size = unit(2, 'mm')) +
  scale_y_discrete(labels = c("D1 \n","D2 \nCD27+IgD-", "D3","D1 \n","D2 \nother", "D3"))+
  guides(size = "none",color = "none")

markers.CD27IgD.PROT.neg <- FindMarkers(seu.live, ident.1 = "CD27+IgD-", ident.2 = "other", assay = "PROT", logfc.threshold = 0.01, test.use = "wilcox", only.pos = F)
markers.CD27IgD.PROT.neg <- filter(markers.CD27IgD.PROT.neg, p_val <= 0.005, avg_log2FC < -0.01)

p.CD27IgD.dotplot.PROTsign.neg <-DotPlot(seu.live,assay = "PROT", features = rev(rownames(markers.CD27IgD.PROT.neg)), split.by = "donor", group.by =  "CD27_IgD_state", cols = "RdBu",dot.scale = 3, scale.max = 100, scale.min = 0, scale = T, scale.by = "radius", col.min = -0.5, col.max = 0.5) +
  cowplot::theme_cowplot() + 
  labs(title="Differential downregulated \nSurface protein markers",x="Surface Proteins \n(p-val < 0.005, logfc <= -0.01)", y = "")+
  coord_flip()+
    add.textsize +
  theme(legend.position = "right", legend.key.size = unit(2, 'mm')) +
  scale_y_discrete(labels = c("D1 \n","D2 \nCD27+IgD-", "D3","D1 \n","D2 \nother", "D3"))+
    guides(size =  "none",color = guide_colorbar(title = "Average scaled\nexpression"))


## Intracellular proteins
seu.fix <- SetIdent(seu.fix,value = "CD27_IgD_state")

markers.CD27IgD.PROT.intra <- FindMarkers(seu.fix, ident.1 = "CD27+IgD-", ident.2 = "other", assay = "PROT", logfc.threshold = 0.05, test.use = "wilcox", only.pos = T)
markers.CD27IgD.PROT.intra <- filter(markers.CD27IgD.PROT.intra, p_val <= 0.05)


p.CD27IgD.dotplot.PROT.intra.sign <- DotPlot(seu.fix,assay = "PROT", features = rev(rownames(markers.CD27IgD.PROT.intra)), split.by = "donor", group.by =  "CD27_IgD_state", cols = "RdBu",dot.scale = 3, scale.max = 100, scale.min = 0, scale = T, scale.by = "radius", col.min = -0.5, col.max = 0.5) +
  cowplot::theme_cowplot() + 
  labs(title="Differential upregulated \nIntracellular protein markers",x="Intracellular proteins \n(p-val < 0.005, logfc >= 0.01)", y = "")+
  coord_flip()+
    add.textsize +
  theme(legend.position = "none", legend.key.size = unit(2, 'mm')) +
  scale_y_discrete(labels = c("D2 \nCD27+IgD-", "D3","D2 \nother", "D3"))

markers.CD27IgD.PROT.intra.neg <- FindMarkers(seu.fix, ident.1 = "CD27+IgD-", ident.2 = "other", assay = "PROT", logfc.threshold = 0.01, test.use = "wilcox", only.pos = F)
markers.CD27IgD.PROT.intra.neg <- filter(markers.CD27IgD.PROT.intra.neg, p_val <= 0.005, avg_log2FC < -0.01)

p.CD27IgD.dotplot.PROTsign.intra.neg <-DotPlot(seu.fix,assay = "PROT", features = rev(rownames(markers.CD27IgD.PROT.intra.neg)), split.by = "donor", group.by =  "CD27_IgD_state", cols = "RdBu",dot.scale = 3, scale.max = 100, scale.min = 0, scale = T, scale.by = "radius", col.min = -0.5, col.max = 0.5) +
  cowplot::theme_cowplot() + 
  labs(title="Differential downregulated \nSurface protein markers",x="Intracellular Proteins \n(p-val < 0.005, logfc <= -0.01)", y = "")+
  coord_flip()+
    add.textsize +
  theme(legend.position = "right", legend.key.size = unit(2, 'mm')) +
  scale_y_discrete(labels = c("D1 \n","D2 \nCD27+IgD-", "D3","D1 \n","D2 \nother", "D3"))+
    guides(size =  "none",color = guide_colorbar(title = "Average scaled\nexpression"))



markers.forviolin <- c("CD138","BLIMP1","IRF4", "IRF8", "XBP1", "CD27","CD24")[c("CD138","BLIMP1","IRF4", "IRF8", "XBP1", "CD27","CD24") %in% rev(rownames(markers.CD27IgD.PROT.intra))]


p.CD27IgD.Vln.PROT.markers <- VlnPlot(seu.fix, assay="PROT", features = rev(c(markers.forviolin)), split.by = "donor", group.by =  "CD27_IgD_state", ncol =3, cols = c(colors.donors[2:3],colors.donors[2:3]), combine = T, stack = T, flip = T) +
    add.textsize +
  theme(legend.position = "right", legend.key.size = unit(2, 'mm')) + 
  labs(title="CD27+IgD express plasmablast/cell \nIntracellular protein markers",x="", y = "Expression")  +
  theme(axis.text.x = element_text(angle = 0, hjust = 0.5))

p.CD27IgD.dotplot.PROT.intra.markers <- DotPlot(seu.fix,assay = "PROT", features = rev(c(markers.forviolin)), split.by = "donor", group.by =  "CD27_IgD_state", cols = "RdBu",dot.scale = 3, scale.max = 100, scale.min = 0, scale = T, scale.by = "radius", col.min = -1, col.max = 1) +
  labs(x = "Differentiation markers - Surface Proteins", y = "") +
  cowplot::theme_cowplot() + 
  labs(title="CD27+IgD express plasmablast/cell \nIntracellular protein markers",x="Differentiation proteins", y = "", colour = "Legend Title\n")+
  coord_flip()+
    add.textsize +
  theme(legend.position = "right", legend.key.size = unit(2, 'mm')) +
  scale_y_discrete(labels = c("D2 \nCD27+IgD-", "D3","D2 \nother", "D3")) +
  guides(size = "none",color = guide_colorbar(title = "Average \nExpression"))

Filter & Save Seurat dataobjects

Fixed-cell dataset: (original and CD27+IgD- object )

message("Fixed dataset:")
seu.fix

An object of class Seurat 
10114 features across 1038 samples within 3 assays 
Active assay: PROT (76 features, 76 variable features)
 2 other assays present: RNA, SCT

seu.fix.filtered <- subset(seu.fix, Cellcycle_state == "non-cycling" & CD27_IgD_state == "CD27+IgD-")
message("CD27+IgD- object:")
seu.fix.filtered

An object of class Seurat 
10114 features across 940 samples within 3 assays 
Active assay: PROT (76 features, 76 variable features)
 2 other assays present: RNA, SCT

Live-cell dataset: (original and CD27+IgD- object )

seu.live

An object of class Seurat 
20366 features across 1433 samples within 3 assays 
Active assay: PROT (50 features, 50 variable features)
 2 other assays present: RNA, SCT

seu.live.filtered <- subset(seu.live, Cellcycle_state == "non-cycling" & CD27_IgD_state == "CD27+IgD-")

seu.live.filtered

An object of class Seurat 
20366 features across 1231 samples within 3 assays 
Active assay: PROT (50 features, 50 variable features)
 2 other assays present: RNA, SCT

Save seurat objects with all cells & information on CC , CD27, and Ig-subclasses.

#saveRDS(seu.fix, "output/seu.fix_norm_cellstate.rds")
#saveRDS(seu.live, "output/seu.live_norm_cellstate.rds")

seu.fix <- seu.fix.filtered
seu.live <- seu.live.filtered

Ig-classifaction

Seurat, protein-based

The supernatant of differentiated cells shows excretion of IgA and IgG antibodies, showing these cell-types should be within the population at 11 days. Using protein-measurements from antibodies against the major Ig-classes, and (in live-cell dataset) the Ig-genes, cells are clustered and represented in a UMAP.

## Run PCA on Ig proteins
seu.fix <- RunPCA(seu.fix, reduction.name = 'pcaIG', features = c("IgM","IgA","IgG", "IgD", "IgE"))

seu.fix <- RunUMAP(seu.fix,dims = 1:4, reduction = "pcaIG", reduction.name = "IGUMAP")
seu.fix <- FindNeighbors(seu.fix, reduction = "pcaIG",assay = "PROT", dims = 1:4, graph.name = "pcaIG_nn")
seu.fix <- FindClusters(seu.fix, graph.name = "pcaIG_nn",  algorithm = 3,resolution = 0.2, verbose = FALSE)
seu.fix[["clusters_pcaIG"]] <- Idents(object = seu.fix)


seu.fix <- RenameIdents(object = seu.fix, `0` = "IgM", `1` = "IgG", `2` = "IgA")
seu.fix[["clusters_pcaIG_named"]] <- Idents(object = seu.fix)

## ---- live cells

## Run PCA on Ig proteins and on Ig genes
seu.live <- RunPCA(seu.live, assay = "PROT",reduction.name = 'pcaIGPROT', features = c("IgM","IgA","IgG", "IgD", "IgE"))




seu.live <- RunUMAP(seu.live,dims = 1:4, reduction = "pcaIGPROT", reduction.name = "IGUMAP")
seu.live <- FindNeighbors(seu.live, reduction = "pcaIGPROT",assay = "PROT", dims = 1:4, graph.name = "pcaIG_nn")
seu.live <- FindClusters(seu.live, graph.name = "pcaIG_nn",  algorithm = 3,resolution = 0.2, verbose = FALSE)
seu.live[["clusters_pcaIG"]] <- Idents(object = seu.live)


seu.live <- RenameIdents(object = seu.live, `0` = "IgM", `1` = "IgG", `2` = "IgA")
seu.live[["clusters_pcaIG_named"]] <- Idents(object = seu.live)

### Plots Ig-protein based clustering
p.umap.Ig.prot.fix <- FeaturePlot(seu.fix, features = c("IgM","IgA","IgG"),slot = "scale.data", reduction = 'IGUMAP', max.cutoff = 2, 
                  cols = c("lightgrey","deeppink3"), ncol = 3, pt.size = 0.9)  &
  labs(x = "UMAP 1", y = "UMAP 2", color = "Scaled \ncounts") &
  theme_half_open()& 
  theme(legend.position = c(0.85,0.85), legend.key.size = unit(2, 'mm')) &
  theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          text = element_text(size=7), axis.ticks = element_blank(),
          axis.text=element_text(size=7),
          plot.title = element_text(size=7, face = "bold",hjust = 0.5))

p.clusters.umap.fix <- DimPlot(seu.fix, reduction = 'IGUMAP',group.by ="clusters_pcaIG_named" , label = TRUE, repel = TRUE, label.size = 2.5, cols = c("#009E73", "#D55E00","#F0E442")
)   &
  labs(x = "UMAP 1", y = "UMAP 2", color = "Clusters", title = "Fixed cells clusters based on \nIg-protein marker expression") &
  theme_half_open()& 
  theme(legend.position = "none", legend.key.size = unit(2, 'mm')) &
  theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          text = element_text(size=7), axis.ticks = element_blank(),
          axis.text=element_text(size=7),
          plot.title = element_text(size=6, face = "bold"))

p.umap.Ig.prot.live <- FeaturePlot(seu.live, features = c("IgM","IgA","IgG"),slot = "scale.data", reduction = 'IGUMAP', max.cutoff = 2, 
                  cols = c("lightgrey","deeppink3"), ncol = 3, pt.size = 0.5)  &
  labs(x = "UMAP 1", y = "UMAP 2", color = "") &
  theme_half_open()& 
  theme(legend.position = c(0.85,0.95), legend.key.size = unit(2, 'mm')) &
  theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          text = element_text(size=7), axis.ticks = element_blank(),
          axis.text=element_text(size=7),
          plot.title = element_text(size=7, face = "bold",hjust = 0.5))

p.umap.Ig.RNA.live <- FeaturePlot(seu.live, features = c("IGHM","IGHA1","IGHG1"),slot = "scale.data", reduction = 'IGUMAP', max.cutoff = 3, 
                  cols = c("lightgrey","dodgerblue3"), ncol = 3, pt.size =  0.5)  &
  labs(x = "UMAP 1", y = "UMAP 2", color = "") &
  theme_half_open()& 
  theme(legend.position = c(0.85,0.95), legend.key.size = unit(2, 'mm')) &
  theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          text = element_text(size=7), axis.ticks = element_blank(),
          axis.text=element_text(size=7),
          plot.title = element_text(size=7, face = "bold",hjust = 0.5))
          
p.clusters.umap.live <- DimPlot(seu.live, reduction = 'IGUMAP', label = TRUE, repel = TRUE, label.size = 2.5,pt.size =  0.5, cols = c("#009E73", "#D55E00","#F0E442","grey")
)   &
  labs(x = "UMAP 1", y = "UMAP 2", color = "Clusters", title = "Igsubtype classification based on \nProtein modality") &
  theme_half_open()& 
  theme(legend.position = "none", legend.key.size = unit(1, 'mm')) &
  theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          text = element_text(size=7), axis.ticks = element_blank(),
          axis.text=element_text(size=7),
          plot.title = element_text(size=6, face = "bold"))

p.donor.umap.live <- DimPlot(seu.live, reduction = 'IGUMAP', group.by = "donor",label = F, pt.size =  0.5,repel = TRUE, label.size = 2.5, cols = colors.donors)   &
  labs(x = "UMAP 1", y = "UMAP 2", color = "Donor", title = "Ig-subtypes classified across \nthree donors") &
  theme_half_open()& 
  theme(legend.position = c(0.85, 0.95), legend.key.size = unit(1, 'mm'), legend.text =element_text(size=6), legend.title =element_text(size=7)  ) &
  theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          text = element_text(size=7), axis.ticks = element_blank(),
          axis.text=element_text(size=7),
          plot.title = element_text(size=6, face = "bold"))

fix.meta <- seu.fix.filtered@meta.data %>%
  mutate(cell = rownames(seu.fix.filtered@meta.data),
         dataset = "fix") 
live.meta <- seu.live.filtered@meta.data%>%
  mutate(cell = rownames(seu.live.filtered@meta.data),
         dataset = "live")

MOFA analysis

library(MOFA2)

fix.DT.RNA <- seu.fix.filtered@assays$SCT@scale.data[rownames(seu.fix.filtered@assays$SCT@scale.data) %chin% VariableFeatures(seu.fix.filtered,assay = "SCT"),]
fix.var.genes <- rownames(fix.DT.RNA)

fix.DT.PROT <- seu.fix.filtered@assays$PROT@scale.data[rownames(seu.fix.filtered@assays$PROT@scale.data) %chin% VariableFeatures(seu.fix.filtered,assay = "PROT"),]
fix.var.prot <- rownames(fix.DT.PROT)

live.DT.RNA <- seu.live.filtered@assays$SCT@scale.data[rownames(seu.live.filtered@assays$SCT@scale.data) %chin% VariableFeatures(seu.live.filtered,assay = "SCT"),]
live.var.genes <- rownames(live.DT.RNA)

live.DT.PROT <- seu.live.filtered@assays$PROT@scale.data[rownames(seu.live.filtered@assays$PROT@scale.data) %chin% VariableFeatures(seu.live.filtered,assay = "PROT"),]
live.var.prot <- rownames(live.DT.PROT)

fix.meta <- seu.fix.filtered@meta.data %>%
  mutate(cell = rownames(seu.fix.filtered@meta.data),
         dataset = "fix") 
live.meta <- seu.live.filtered@meta.data%>%
  mutate(cell = rownames(seu.live.filtered@meta.data),
         dataset = "live")

## 
common.proteins <- c("TACI","CD70","IgD","CD86","CD5","PDL1","PDL2","CD21","IgM","CD11b","CD53","CD45","CD80","IgG","IgA","IgE","BAFFR")

fix.DT.RNA <- data.table(fix.DT.RNA) %>%
  mutate(gene = fix.var.genes)

fix.DT.PROT <- data.table(fix.DT.PROT) %>%
  mutate(gene = fix.var.prot)
fix.DT.PROT.unique <- fix.DT.PROT[!(fix.var.prot %chin% common.proteins),] %>%
  mutate(gene = paste0(gene, "_f"))

fix.DT.PROT.common <- fix.DT.PROT[(fix.var.prot %chin% common.proteins),]
fix.DT.PROT.rename <- bind_rows(fix.DT.PROT.unique,fix.DT.PROT.common)
#rownames(fix.DT.PROT.rename) <- fix.DT.PROT.rename$gene

live.DT.RNA <- data.table(live.DT.RNA) %>%
  mutate(gene = live.var.genes)

live.DT.PROT <- data.table(live.DT.PROT) %>%
  mutate(gene = live.var.prot)
live.DT.PROT.unique <- live.DT.PROT[!(live.var.prot %chin% common.proteins),] %>%
  mutate(gene = paste0(gene, "_l"))
live.DT.PROT.common <- live.DT.PROT[(live.var.prot %chin% common.proteins),]
live.DT.PROT.rename <- bind_rows(live.DT.PROT.unique,live.DT.PROT.common)
## merge data
meta.all <- merge.data.table(x = fix.meta, y=live.meta, all = TRUE)
meta.all$sample <- meta.all$cell
meta.all$group <- meta.all$donor

RNA.all <- merge.data.table(x = fix.DT.RNA, y=live.DT.RNA, all = TRUE)
rownames(RNA.all) <- RNA.all$gene
RNA.all <- RNA.all[,gene:=NULL]

PROT.all <- merge.data.table(x = fix.DT.PROT.rename, y=live.DT.PROT.rename, all = TRUE)

rownames(PROT.all) <- PROT.all$gene
PROT.all <- PROT.all[,gene:=NULL]

RNA.live.only <- RNA.all[rownames(RNA.all) %chin% VariableFeatures(seu.live.filtered, assay = "SCT"),]
genes.RNA.live.only <- rownames(RNA.all)[rownames(RNA.all) %chin% VariableFeatures(seu.live.filtered, assay = "SCT")]
RNA.live.only <- as.matrix(RNA.live.only)
#genes.RNA.live.only <- tolower(genes.RNA.live.only)
rownames(RNA.live.only) <- genes.RNA.live.only
#RNA.live.only[,c(which(meta.all$dataset== "fix"))] <- NA #?

## split data tables

### PROT

PROT.common <- PROT.all[rownames(PROT.all) %chin% common.proteins,]
list.common.prot <- rownames(PROT.all)[rownames(PROT.all) %chin% common.proteins]
PROT.common <- as.matrix(PROT.common, rownames.value    = list.common.prot)
#rownames(PROT.common) <- paste0(rownames(PROT.common),"_PROT")
#PROT.IG <- PROT.all[rownames(PROT.all) %chin% c("IgA", "IgM", "IgG", "IgD", "IgE"),]
#list.IG.prot <- rownames(PROT.all)[rownames(PROT.all) %chin% c("IgA", "IgM", "IgG", "IgD", "IgE")]
#PROT.IG <- as.matrix(PROT.IG, rownames.value   = list.IG.prot)

PROT.fix <- PROT.all[c(rownames(PROT.all) %chin% fix.DT.PROT.unique$gene),]
list.fix.prot <-rownames(PROT.all)[c(rownames(PROT.all) %chin% fix.DT.PROT.unique$gene)]
PROT.fix <- as.matrix(PROT.fix, rownames.value  = list.fix.prot)
#rownames(PROT.fix) <- paste0(rownames(PROT.fix),"_PROT")

PROT.live <- PROT.all[c(rownames(PROT.all) %chin% live.DT.PROT.unique$gene) ,]
list.live.prot <-  rownames(PROT.all)[c(rownames(PROT.all) %chin% live.DT.PROT.unique$gene) ]
PROT.live <- as.matrix(PROT.live, rownames.value    = list.live.prot)


### RNA 

RNA.common <- RNA.all[rownames(RNA.all) %chin% rownames(seu.fix.filtered[["SCT"]])[rownames(seu.fix.filtered[["SCT"]]) %chin% rownames(seu.live.filtered[["SCT"]])],]
list.common.RNA <- rownames(RNA.all)[rownames(RNA.all) %chin% rownames(seu.fix.filtered[["SCT"]])[rownames(seu.fix.filtered[["SCT"]]) %chin% rownames(seu.live.filtered[["SCT"]])]]
list.common.RNA <- tolower(list.common.RNA)
RNA.common <- as.matrix(RNA.common, rownames.value  = list.common.RNA)

#RNA.IG <- RNA.all[rownames(RNA.all) %chin% c("IgA", "IgM", "IgG", "IgD", "IgE"),]
#list.IG.RNA <- rownames(RNA.all)[rownames(RNA.all) %chin% c("IgA", "IgM", "IgG", "IgD", "IgE")]
#RNA.IG <- as.matrix(RNA.IG, rownames.value = list.IG.RNA)

RNA.fix <- RNA.all[!c(rownames(RNA.all) %chin% list.common.RNA) & rownames(RNA.all) %chin% rownames(seu.fix.filtered[["SCT"]]),]
list.fix.RNA <-rownames(RNA.all)[!c(rownames(RNA.all) %chin% list.common.RNA) & rownames(RNA.all) %chin% rownames(seu.fix.filtered[["SCT"]])]
list.fix.RNA <- tolower(list.fix.RNA)
RNA.fix <- as.matrix(RNA.fix, rownames.value    = list.fix.RNA)

RNA.live <- RNA.all[!c(rownames(RNA.all) %chin% list.common.RNA) & !c(rownames(RNA.all) %chin% list.fix.RNA) ,]
list.live.RNA <-  rownames(RNA.all)[!c(rownames(RNA.all) %chin% list.common.RNA) & !c(rownames(RNA.all) %chin% list.fix.RNA) ]
#list.live.RNA <- tolower(list.live.RNA)
RNA.live <- as.matrix(RNA.live, rownames.value  = list.live.RNA)

###
all.features <- c(rownames(RNA.live.only),rownames(PROT.common), rownames(PROT.live), rownames(PROT.fix))
duplicated.features <- all.features[duplicated(all.features)]

RNA.live.only <- data.table(RNA.live.only) %>%
  mutate(gene = rownames(RNA.live.only)) 

rownames.rna.live.only <- RNA.live.only$gene
RNA.live.only <- RNA.live.only %>%
  dplyr::select(-c("gene"))
RNA.live.only <- as.matrix(RNA.live.only)
#genes.RNA.live.only <- tolower(genes.RNA.live.only)
rownames(RNA.live.only) <- rownames.rna.live.only

all.features <- c(rownames(RNA.live.only),rownames(PROT.common), rownames(PROT.live), rownames(PROT.fix))
duplicated.features <- all.features[duplicated(all.features)]

myfiles <- list.files(path="output/", pattern = ".rds$")

if("MOFA_analysis_Donorgroup.rds" %in%  myfiles){mofa <- readRDS("output/MOFA_analysis_Donorgroup.rds")} else { #If so, read object, else do:
     
      mofa <- create_mofa(data = list(RNA = RNA.live.only,  PROT.common = PROT.common,PROT.live = PROT.live, PROT.fix = PROT.fix ), groups = meta.all$group)

      # Default settings used (try 15 factors, excludes all non-informative factors)
      data_opts <- get_default_data_options(mofa)
      model_opts <- get_default_model_options(mofa)
      train_opts <- get_default_training_options(mofa)
      train_opts$seed <- 42 # use same seed for reproducibility
      mofa <- prepare_mofa(
        object = mofa,
        data_options = data_opts,
        model_options = model_opts,
        training_options = train_opts
        )

    mofa <- run_mofa(mofa, outfile = "output/MOFA_analysis_Donorgroup.hdf5", use_basilisk = TRUE)
    mofa <- run_umap(mofa)
    samples_metadata(mofa) <- meta.all
    #saveRDS(mofa, file= "output/MOFA_analysis_Donorgroup.rds")
  
}
mofa

Trained MOFA with the following characteristics: 
 Number of views: 4 
 Views names: RNA PROT.common PROT.live PROT.fix 
 Number of features (per view): 3000 17 33 59 
 Number of groups: 3 
 Groups names: D1 D2 D3 
 Number of samples (per group): 433 871 867 
 Number of factors: 9 

MOFA description

p.mofaoverview.input <- plot_data_overview(mofa)

## Variance per factor
p.variance.perfactor.all <- plot_variance_explained(mofa, x="view", y="factor") +
    labs(title = "Variance explained by each factor per modality")  +
  RotatedAxis()

## variance total
p.variance.total <- plot_variance_explained(mofa, x="view", y="factor", plot_total = T)
# 
p.variance.total <- plot_variance_explained(mofa, x="group", y="factor", plot_total = T)
p.variance.total <- p.variance.total[[2]] +
  add.textsize +
  labs(title = "Total variance per modality") +
  geom_text(aes(label=round(R2,1)), vjust=1.6, color="white", size=2.5) +
  theme_classic()+
  RotatedAxis()

factors.selected <- paste0("Factor",1:5)

plot.rank.PROT <- plot_weights(mofa, 
  view = "PROT.common", 
  factors = factors.selected, 
  nfeatures = 10, 
  text_size =3
) +   
  labs(title = "Top 10  Common Protein loadings per factor") +
  theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank()
)

p.factors.Igclasses.PROT<- plot_factors(mofa, 
  factors = c(1:2),
  color_by = "IgM",dot_size = 1.2
) +
  labs(title="MOFA factor 1 & 2 distinguish Ig-classes \nProtein expression")+
  plot_factors(mofa, 
  factors =c(1:2),
  color_by = "IgA",dot_size = 1.2
)+
  plot_factors(mofa, 
  factors = c(1:2),
  color_by = "IgG",dot_size = 1.2
) &
    add.textsize &  
  labs(x = "Factor 1", y = "Factor 2", color = "Scaled \ncounts") &
  theme_half_open()& 
  theme(legend.position = c(0.85,0.95), legend.key.size = unit(2, 'mm')) &
  theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          text = element_text(size=7), axis.ticks = element_blank(),
          axis.text=element_text(size=7),
          plot.title = element_text(size=7, face = "bold",hjust = 0)) &
  scale_fill_gradient(low = "white", high = "deeppink3")

p.factors.Igclasses.mRNA <- plot_factors(mofa, 
  factors =c(1:2),
  color_by = "IGHM",dot_size = 1.2, show_missing = F,alpha = 0.5
) +labs(title="mRNA expression")+
  plot_factors(mofa, 
  factors = c(1:2),
  color_by = "IGHA1",dot_size = 1.2, show_missing = F,alpha = 0.5
) +
  plot_factors(mofa, 
  factors =c(1:2),
  color_by = "IGHG1",dot_size = 1.2, show_missing = F,alpha = 0.5
) &
    add.textsize &  
  labs(x = "Factor 1", y = "Factor 2", color = "Scaled \ncounts") &
  theme_half_open()& 
  theme(legend.position = c(0.85,0.95), legend.key.size = unit(2, 'mm')) &
  theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          text = element_text(size=7), axis.ticks = element_blank(),
          axis.text=element_text(size=7),
          plot.title = element_text(size=7, face = "bold",hjust = 0))  &
  scale_fill_gradient(low = "white", high = "dodgerblue3")

Clustering Ig-classes

clusters <- cluster_samples(mofa, k=3, factors=c(1,2) )

meta.all.cluster <-  meta.all[,] %>% 
  left_join(data.frame(sample = names(clusters$cluster), cluster = clusters$cluster))
meta.all.cluster$cluster <- as.factor(meta.all.cluster$cluster)
samples_metadata(mofa) <-meta.all.cluster

clustermeta <- mofa@samples_metadata %>%
  mutate(IgClass = ifelse(cluster == "1", "IgM", ifelse(cluster == "3", "IgA", "IgG")))

samples_metadata(mofa) <- clustermeta

Percentage_cluster <- clustermeta %>% dplyr::count(IgClass = factor(IgClass),group = factor(donor)) %>% 
  group_by(group)%>%
    mutate(pct = prop.table(n)) 

p.percentages.Igclasses.MOFAclusters <- ggplot(Percentage_cluster, aes(x = IgClass , y = pct, fill = group, label = scales::percent(pct,accuracy = 0.1))) + 
    geom_col(position = 'dodge') + 
    geom_text(position = position_dodge(width = .9),    # move to center of bars
              vjust = -0.5,    # nudge above top of bar
              size = 1.8) + 
    scale_y_continuous(labels = scales::percent) +
  cowplot::theme_cowplot() + 
  labs(title="Percentage Ig-class",x="", y = "Percentage (%)")+
  scale_x_discrete(labels = c("1" = "IgG", "2" ,"3")) +
scale_color_manual("Donor", values=c("#999999", "#E69F00", "#56B4E9"))+
scale_fill_manual("Donor",values=c("#999999", "#E69F00", "#56B4E9"))+
  add.textsize

p.MOFA.factors.dataset <- plot_factors(mofa, 
  factors = c(1:2),
  color_by = "dataset",
  alpha = 0.8,
  dot_size = 1.2
) +
  labs(title = "MOFA analysis integrates \nthe two datasets") +
scale_fill_manual("Dataset",values=c("lightgrey", "black")) &
    add.textsize &  
  labs(x = "Factor 1", y = "Factor 2", color = "Scaled \ncounts") &
  theme_half_open()& 
  theme(legend.position = c(0.85,0.95), legend.key.size = unit(2, 'mm')) &
  theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          text = element_text(size=7), axis.ticks = element_blank(),
          axis.text=element_text(size=7),
          plot.title = element_text(size=7, face = "bold",hjust = 0)) 
  
p.MOFA.factors.donors <- plot_factors(mofa, 
  factors = c(1:2),
  color_by = "donor",
  alpha = 0.8,
  dot_size = 1.2
) +labs(title = "MOFA factors are computed \nacross three donors") +
scale_fill_manual("Donor",values=colors.donors) &
    add.textsize &  
  labs(x = "Factor 1", y = "Factor 2", color = "Donor") &
  theme_half_open()& 
  theme(legend.position = c(0.85,0.95), legend.key.size = unit(2, 'mm')) &
  theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          text = element_text(size=7), axis.ticks = element_blank(),
          axis.text=element_text(size=7),
          plot.title = element_text(size=7, face = "bold",hjust = 0)) 


p.MOFA.factors.cluster <-  plot_factors(mofa, 
   factors =c(1:2),
  color_by = "IgClass",
  dot_size = 1.2
) +labs(title = "k-means clustering on Factor 1 & 2") +
scale_fill_manual("Cluster",values= colors.clusters) &
  labs(x = "Factor 1", y = "Factor 2", color = "Donor") &
  theme_half_open()& 
  theme(legend.position = c(0.85,0.95), legend.key.size = unit(2, 'mm')) &
  theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          text = element_text(size=10), axis.ticks = element_blank(),
          axis.text=element_text(size=10),
          plot.title = element_text(size=10, face = "bold",hjust = 0)) 

p.loadings.scatter <- plot_weights_scatter(mofa, factors = c(1,2), view = "PROT.common")+
  ggrepel::geom_text_repel(aes(x=x, y=y, label=feature), size=2.5, max.overlaps = 10)+
   coord_cartesian(xlim =c(-0.5,0.5), ylim = c(-0.5,0.5)) +

plot_weights_scatter(mofa, factors = c(1,2), view = "PROT.fix")+
  ggrepel::geom_text_repel(aes(x=x, y=y, label=feature), size=2.5, max.overlaps = 10)+
   coord_cartesian(xlim =c(-0.6,0.6), ylim = c(-0.6,0.6)) +

plot_weights_scatter(mofa, factors = c(1,2), view = "PROT.live")+
  ggrepel::geom_text_repel(aes(x=x, y=y, label=feature), size=2.5, max.overlaps = 10) +

plot_weights_scatter(mofa, factors = c(1,2), view = "RNA")+
  ggrepel::geom_text_repel(aes(x=x, y=y, label=feature), size=2, max.overlaps =20) +
   coord_cartesian(xlim =c(-0.002,0.002), ylim = c(-0.002,0.002)) &
  add.textsize

p.topweights.positive <- plot_top_weights(mofa, factors = c(1,2), sign = "positive", view = "RNA", nfeatures = 35) +
  plot_top_weights(mofa, factors = c(1,2), sign = "positive", view = "PROT.common") /
  plot_top_weights(mofa, factors =c(1,2), sign = "positive", view = "PROT.fix")/
  plot_top_weights(mofa, factors = c(1,2), sign = "positive", view = "PROT.live") &
  add.textsize

p.topweights.negative <- plot_top_weights(mofa, factors = c(1,2), sign = "negative", view = "RNA", nfeatures = 35) +
  plot_top_weights(mofa, factors = c(1,2), sign = "negative", view = "PROT.common") /
  plot_top_weights(mofa, factors = c(1,2), sign = "negative", view = "PROT.fix")/
  plot_top_weights(mofa, factors = c(1,2), sign = "negative", view = "PROT.live") &
  add.textsize

p.MOFA.factors.clusters <- plot_factors(mofa, 
  factors = c(1:2),
  color_by = "dataset",
  alpha = 0.5,
  dot_size = 1.2
) +
  scale_color_manual("Dataset", values=c("yellow", "blue"))+
  labs(title = "MOFA factors integrate two \n datasets across donors") +
scale_fill_manual("Donor",values=c("yellow", "blue")) +
  
plot_factors(mofa, 
  factors = c(1:2),
  color_by = "donor",
  alpha = 0.8,
  dot_size = 1.2
) +scale_color_manual("Donor", values=c("#999999", "#E69F00", "#56B4E9"))+
  labs(title = "Donors overlab (no batch/donor effects)") +
scale_fill_manual("Donor",values=c("#999999", "#E69F00", "#56B4E9")) +
  plot_factors(mofa, 
   factors =c(1:2),
  color_by = "IgClass",
  dot_size = 1.2
) +
  labs(title = "K-means clustering on Factor 1&2")+ scale_color_manual("Ig-class", values=c("#009E73", "#D55E00", "#CC79A7"))+
scale_fill_manual("Ig-class",values=c(c("#009E73", "#D55E00", "#CC79A7"))) & add.textsize

Add clusters & Factors to Seurat object

fix.meta2 <- fix.meta %>%
  mutate(sample_cell = rownames(fix.meta))

clustermeta2 <- clustermeta %>%
  mutate(sample_cell = sample) %>%
  dplyr::select(c(sample_cell,IgClass, clusterMOFA = cluster))

factorvalues <- get_factors(mofa,factors = c(1:2),as.data.frame = TRUE )

factorvalues <- factorvalues %>%
  dplyr::select(-group)%>%
  spread(factor,value)%>%
  mutate(sample_cell = sample) %>%
  dplyr::select(c(sample_cell,Factor1, Factor2))

fix.meta.new <- left_join(fix.meta2, clustermeta2)
fix.meta.new <- left_join(fix.meta.new,factorvalues)
rownames(fix.meta.new) <- fix.meta.new$sample_cell
seu.fix <- AddMetaData(seu.fix, fix.meta.new)
seu.fix <- SetIdent(seu.fix,value = "IgClass")
Idents(seu.fix) <- factor(x = Idents(seu.fix), levels = c("IgM", "IgG", "IgA"))

live.meta2 <- live.meta %>%
  mutate(sample_cell = rownames(live.meta))

clustermeta2 <- clustermeta %>%
  mutate(sample_cell = sample) %>%
  dplyr::select(c(sample_cell,IgClass))

factorvalues <- get_factors(mofa,factors = c(1:2),as.data.frame = TRUE )

factorvalues <- factorvalues %>%
  dplyr::select(-group)%>%
  spread(factor,value)%>%
  mutate(sample_cell = sample) %>%
  dplyr::select(c(sample_cell,Factor1, Factor2))

live.meta.new <- left_join(live.meta2, clustermeta2)
live.meta.new <- left_join(live.meta.new,factorvalues)
rownames(live.meta.new) <- live.meta.new$sample_cell
seu.live <- AddMetaData(seu.live, live.meta.new)
seu.live <- SetIdent(seu.live,value = "IgClass")
Idents(seu.live) <- factor(x = Idents(seu.live), levels = c("IgM", "IgG", "IgA"))

Figure 1

fig.1.row1 <- plot_grid(p.CCscore.MKI67.live, 
                        p.CCscore.pRb,
                        p.percentage.cyclinglive,
                     labels = panellabels[4:6], label_size = 10, ncol = 3, rel_widths = c(1.2,1.2,0.9))

fig.1.row2 <- plot_grid(p.scatter.CD27.IgD.live,
                        p.CD27IgD.dotplot.PROT.intra.markers,
                        p.CD27IgD.Vln.gene.markers,
                     labels = panellabels[7:9], label_size = 10, ncol = 3, rel_widths = c(1.2,1.2,0.9))


Fig.1.full <- plot_grid(fig.1.row1, fig.1.row2, labels = "", label_size = 10, ncol = 1, rel_heights = c(1,1))

# #ggsave(Fig.1.full, filename = "output/paper_figures/Fig1D-I.eps", width = 177, height = 120, units = "mm",  dpi = 300, useDingbats = FALSE)
# #ggsave(Fig.1.full, filename = "output/paper_figures/Fig1D-I.TIFF", width = 177, height = 120, units = "mm",  dpi = 300)

Fig.1.full

Figure 2

fig.2.row1 <- plot_grid( p.factors.Igclasses.PROT,p.MOFA.factors.dataset, labels = panellabels[c(1,2)], label_size = 10, rel_widths = c(2.1,0.9))

fig.2.row2 <- plot_grid( p.factors.Igclasses.mRNA,p.MOFA.factors.donors,labels = c("", panellabels[3]), label_size = 10,
          rel_widths = c(2.1,0.9))


fig.2.full <- plot_grid(fig.2.row1, fig.2.row2, labels = "", label_size = 10, ncol = 1, rel_heights = c(1,1))

# #ggsave(fig.2.full, filename = "output/paper_figures/Fig2.eps", width = 177, height = 120, units = "mm",  dpi = 300, useDingbats = FALSE)
# #ggsave(fig.2.full, filename = "output/paper_figures/Fig2.TIFF", width = 177, height = 120, units = "mm",  dpi = 300)

fig.2.full

Supplementary figures

Supplement CD27+IgD- subclasses (fig1)

p.percentages <- plot_grid(p.percentage.cycling.fix,
          p.percentage.class_switch_plasma.live,
          p.percentage.class_switch_plasma.fix,
          labels = panellabels[1:3], label_size = 10, ncol = 3, rel_widths = c(1,1))

p.dotplot.genes<- plot_grid(p.CD27IgD.dotplot.genesign,
          p.CD27IgD.dotplot.genesign.neg,
          labels = c(panellabels[4:5]), label_size = 10, ncol = 2, rel_heights =  c(1,1))


p.dotplot.proteins.surface <- plot_grid(p.CD27IgD.dotplot.PROTsign,
          p.CD27IgD.dotplot.PROTsign.neg,
          labels = c(panellabels[6:7]), label_size = 10, ncol = 2, rel_heights =  c(1,1))

p.dotplot.proteins.intra <- plot_grid(p.CD27IgD.dotplot.PROT.intra.sign,p.CD27IgD.dotplot.PROTsign.intra.neg,
          labels = c(panellabels[8:9]), label_size = 10, ncol = 2, rel_heights =  c(1,1))




p.supplement.percentages.dotplots <- plot_grid(p.percentages,
                                               p.dotplot.genes,
                                               p.dotplot.proteins.surface,
                                               p.dotplot.proteins.intra,
                                               labels = "", label_size = 10, ncol = 1, rel_heights = c(1.5,1.7, 1.25,1.15))

#ggsave(p.supplement.percentages.dotplots, filename = "output/paper_figures/Suppl_CC_CD27_dotplots.eps", width = 177, height = 220, units = "mm",  dpi = 300, useDingbats = FALSE)
#ggsave(p.supplement.percentages.dotplots, filename = "output/paper_figures/Suppl_CC_CD27_dotplots.TIFF", width = 177, height = 220, units = "mm",  dpi = 300)

p.supplement.percentages.dotplots

Supplementary figure Additional information supplementing figure 1. a. fixed data percentage cycling. b-c. Percentages CD27+IgD- cells in datasets. d-j Differentially expressed genes and proteins.

Supplement MOFA properties (fig2)

UMAP on shared proteins

library(umap)

common.prot <- left_join(live.DT.PROT.common,fix.DT.PROT.common)

umap.common.prot <- umap(t(common.prot[,-"gene"]), random_state=42)


meta_inorder <- data.frame( sample = colnames(common.prot[,-"gene"]))
meta_inorder <- left_join(meta_inorder, meta.all)
umap.labels <- meta_inorder$dataset

umap.common.prot.withmeta <- cbind(data.frame(umap.common.prot$layout),umap.labels)

reviewer_umap_commonProt <- ggplot(as.data.frame(umap.common.prot.withmeta), aes(x = X1, y = X2, fill = umap.labels ))+
  geom_point(size = 0.9, alpha = 0.8, shape = 21) +
  labs(title = "UMAP on common proteins \nshows full overlap of the two datasets") +
scale_fill_manual("Dataset",values=c("lightgrey", "black")) &
    add.textsize &  
  labs(x = "UMAP 1", y = "UMAP 2", color = "Scaled \ncounts") &
  theme_half_open()& 
  theme(legend.position = c(0.85,0.25), legend.key.size = unit(2, 'mm')) &
  theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          text = element_text(size=9), axis.ticks = element_blank(),
          axis.text=element_text(size=9),
          plot.title = element_text(size=9, face = "bold",hjust = 0)) 
reviewer_umap_commonProt

# 
# ggsave(reviewer_umap_commonProt, filename = "output/paper_figures/reviewer_umap_commonProt.png", width = 70, height = 80, units = "mm",  dpi = 300)
# ggsave(reviewer_umap_commonProt, filename = "output/paper_figures/reviewer_umap_commonProt.eps", width = 70, height = 80, units = "mm",  dpi = 300)

p.mofa.suppl <- plot_grid(
p.mofaoverview.input,p.variance.perfactor.all,reviewer_umap_commonProt,NULL,plot.rank.PROT,p.MOFA.factors.cluster,labels = c(panellabels[c(1:2)],panellabels[4],"",panellabels[c(3,5)]), 
          label_size = 10,ncol = 3,
          rel_widths = c(0.5,2,1), rel_heights = c(1,1))

 ggsave(p.mofa.suppl, filename = "output/paper_figures/Suppl_mofa.pdf", width = 483, height = 300, units = "mm",  dpi = 300, useDingbats = FALSE)
 ggsave(p.mofa.suppl, filename = "output/paper_figures/Suppl_mofa.png", width = 483, height = 300, units = "mm",  dpi = 300)

p.mofa.suppl

Supplement Seurat protein clustering (fig2)

fig.2.suppl.row1 <- plot_grid(p.clusters.umap.live,  p.umap.Ig.prot.live,labels = panellabels[c(1,3)], label_size = 10, rel_widths = c(1,3))

fig.2.suppl.row2 <- plot_grid(p.donor.umap.live, p.umap.Ig.RNA.live,labels = panellabels[2], label_size = 10,
          rel_widths = c(1,3))

p.umap.fix.all <- plot_grid(p.clusters.umap.fix, 
          p.umap.Ig.prot.fix, labels = panellabels[c(4,5)], 
          label_size = 10,
          rel_widths = c(1,3))


fig.2.suppl.full <- plot_grid(fig.2.suppl.row1, fig.2.suppl.row2,p.umap.fix.all, labels = "", label_size = 10, ncol = 1, rel_heights = c(1,1,1))

#ggsave(fig.2.suppl.full, filename = "output/paper_figures/Suppl_SeuratCluster.eps", width = 177, height = 180, units = "mm",  dpi = 300, useDingbats = FALSE)
#ggsave(fig.2.suppl.full, filename = "output/paper_figures/Suppl_SeuratCluster.TIFF", width = 177, height = 180, units = "mm",  dpi = 300)

fig.2.suppl.full

Supplementary figure Fixed cells Ig-classes

Supplement Ig-class % and ELISA

Import ELISA results

#pzfx_tables("data/supplementary/Figure 3/Fig3ELISA/GM6642 Ig ELISA.pzfx")

df.ELISA.IgA <- read_pzfx("data/supplementary/Figure 3/Fig3ELISA/GM6642 Ig ELISA.pzfx", table = "IgA production in time")

df.ELISA.IgA <- df.ELISA.IgA %>% gather(sampleID,signal,2:ncol(df.ELISA.IgA))%>%
  mutate(sample = ROWTITLE,Cytokine = "IgA")

df.ELISA.IgG <- read_pzfx("data/supplementary/Figure 3/Fig3ELISA/GM6642 Ig ELISA.pzfx", table = "IgG production in time")

df.ELISA.IgG <- df.ELISA.IgG %>% gather(sampleID,signal,2:ncol(df.ELISA.IgG))%>%
  mutate(sample = ROWTITLE,Cytokine = "IgG")

df.ELISA.IgM <- read_pzfx("data/supplementary/Figure 3/Fig3ELISA/GM6642 Ig ELISA.pzfx", table = "IgM production in time")

df.ELISA.IgM <- df.ELISA.IgM %>% gather(sampleID,signal,2:ncol(df.ELISA.IgM))%>%
  mutate(sample = ROWTITLE,Cytokine = "IgM")


df.ELISA.all <- rbind(df.ELISA.IgM, df.ELISA.IgA) %>%
  rbind(.,df.ELISA.IgG)%>%
  separate(sampleID,c("sampleID", "rep"), sep = "_") 

df.ELISA.all$sample <- factor(df.ELISA.all$sample, levels = c("Day 0" , "Day 1",  "Day 4",  "Day 7",  "Day 11"))

summarySE <- function(data=NULL, measurevar, groupvars=NULL, na.rm=FALSE,
                      conf.interval=.95, .drop=TRUE) {
    library(plyr)

    # New version of length which can handle NA's: if na.rm==T, don't count them
    length2 <- function (x, na.rm=FALSE) {
        if (na.rm) sum(!is.na(x))
        else       length(x)
    }

    # This does the summary. For each group's data frame, return a vector with
    # N, mean, and sd
    datac <- ddply(data, groupvars, .drop=.drop,
      .fun = function(xx, col) {
        c(N    = length2(xx[[col]], na.rm=na.rm),
          mean = mean   (xx[[col]], na.rm=na.rm),
          sd   = sd     (xx[[col]], na.rm=na.rm)
        )
      },
      measurevar
    )

    # Rename the "mean" column    
    datac <- rename(datac, c("mean" = measurevar))

    datac$se <- datac$sd / sqrt(datac$N)  # Calculate standard error of the mean

    # Confidence interval multiplier for standard error
    # Calculate t-statistic for confidence interval: 
    # e.g., if conf.interval is .95, use .975 (above/below), and use df=N-1
    ciMult <- qt(conf.interval/2 + .5, datac$N-1)
    datac$ci <- datac$se * ciMult

    return(datac)
}

tgc <- summarySE(df.ELISA.all, measurevar="signal", groupvars=c("sample", "Cytokine","sampleID"))


p.ELISA <- ggplot(tgc, aes(x=sample, y=signal, colour=sampleID, group=sampleID)) + geom_point( size=2, shape=21, fill="white") + # 21 is filled circle
    geom_errorbar(aes(ymin=signal-se, ymax=signal+se), colour="black", width=.1) +
    geom_line() +
    
    xlab("") +
    ylab("Concentration (ng/ml)") +
    scale_colour_hue(name="Donor")+
    theme_few()+
  facet_wrap(.~Cytokine, scales = "free_y")+
  add.textsize+ theme(axis.text.x = element_text(angle = 45, hjust = 1))

Suppl.IgG.validation <- plot_grid(p.percentages.Igclasses.MOFAclusters,p.ELISA, labels = panellabels[1:2], label_size = 10, ncol = 2, rel_widths = c(1,2))

#ggsave(Suppl.IgG.validation, filename = "output/paper_figures/Suppl.IgG.validation.eps", width = 177, height = 60, units = "mm",  dpi = 300, useDingbats = FALSE)
#ggsave(Suppl.IgG.validation, filename = "output/paper_figures/Suppl.IgG.validation.TIFF", width = 177, height = 60, units = "mm",  dpi = 300)

Suppl.IgG.validation

PCA without IG’s

## Run PCA without Ig
seu.fix <- RunPCA(seu.fix, reduction.name = 'pcanoIG', features = rownames(seu.fix)[!(rownames(seu.fix) %in% grep(rownames(seu.fix),pattern = "^Ig",value = TRUE,invert = FALSE))])

seu.fix <- RunUMAP(seu.fix,dims = 1:4, reduction = "pcanoIG", reduction.name = "noIGUMAP" )
seu.fix <- FindNeighbors(seu.fix, reduction = "pcanoIG",assay = "PROT", dims = 1:4, graph.name = "pcanoIG_nn")
seu.fix <- FindClusters(seu.fix, graph.name = "pcanoIG_nn",  algorithm = 3,resolution = 0.3, verbose = FALSE)
seu.fix[["clusters_pcanoIG"]] <- Idents(object = seu.fix)


p.umap.Ig.prot.fix2 <- FeaturePlot(seu.fix, features = c("IgM","IgA","IgG"),slot = "scale.data", reduction = 'noIGUMAP', max.cutoff = 2, 
                  cols = c("lightgrey","deeppink3"), ncol = 3, pt.size = 0.5)  &
  labs(x = "UMAP 1", y = "UMAP 2", color = "Scaled \ncounts") &
  theme_half_open()& 
  theme(legend.position = c(0.85,0.85), legend.key.size = unit(2, 'mm')) &
  theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          text = element_text(size=7), axis.ticks = element_blank(),
          axis.text=element_text(size=7),
          plot.title = element_text(size=7, face = "bold",hjust = 0.5))

p.clusters.umap.fix2 <- DimPlot(seu.fix, reduction = 'noIGUMAP',group.by = "clusters_pcaIG_named", label = TRUE, repel = TRUE, label.size = 2.5, cols = c("#009E73", "#D55E00","#F0E442"), pt.size = 0.5,
)   &
  labs(x = "UMAP 1", y = "UMAP 2", color = "Ig-based \nclusters", title = "Fixed cells PCA & UMAP \nwithout Ig-proteins") &
  theme_half_open()& 
  theme(legend.position = "none", legend.key.size = unit(2, 'mm')) &
  theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          text = element_text(size=7), axis.ticks = element_blank(),
          axis.text=element_text(size=7),
          plot.title = element_text(size=6, face = "bold"))



# seu.fix <- RenameIdents(object = seu.fix, `1` = "IgM", `2` = "IgG", `3` = "IgA")
# seu.fix[["clusters_pcaIG_named"]] <- Idents(object = seu.fix)

## ---- live cells

## Run PCA on Ig proteins and on Ig genes
seu.live <- RunPCA(seu.live, assay = "PROT",reduction.name = 'pcanoIGPROT', features = rownames(seu.live@assays$PROT)[!(rownames(seu.live@assays$PROT) %in% grep(rownames(seu.live@assays$PROT),pattern = "^Ig",value = TRUE,invert = FALSE))])


seu.live <- RunPCA(seu.live, assay = "SCT",reduction.name = 'pcanoIGRNA', features = rownames(seu.live@assays$SCT)[!(rownames(seu.live@assays$SCT) %in% grep(rownames(seu.live@assays$SCT),pattern = "^IG",value = TRUE,invert = FALSE))])


seu.live <- RunUMAP(seu.live,dims = 1:10, reduction = "pcanoIGPROT", reduction.name = "pcanoIGPROT_UMAP" )

p.clusters.umap.live_noIg_PROT <- DimPlot(seu.live, reduction = 'pcanoIGPROT_UMAP',group.by = "clusters_pcaIG_named", label = TRUE, repel = TRUE, label.size = 2.5,pt.size =  0.5, cols = c("#009E73", "#D55E00","#F0E442","grey")
)   &
  labs(x = "UMAP 1", y = "UMAP 2", color = "Clusters", title = "Live cells UMAP based on \nprotein modality PCA \nwithout Ig-proteins") &
  theme_half_open()& 
  theme(legend.position = "none", legend.key.size = unit(1, 'mm')) &
  theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          text = element_text(size=7), axis.ticks = element_blank(),
          axis.text=element_text(size=7),
          plot.title = element_text(size=6, face = "bold"))

p.umap.Ig.live_noIg_PROT <- FeaturePlot(seu.live, features = c("IgM","IgA","IgG"),slot = "scale.data", reduction = 'pcanoIGPROT_UMAP', max.cutoff = 2, 
                  cols = c("lightgrey","deeppink3"), ncol = 3, pt.size = 0.5)  &
  labs(x = "UMAP 1", y = "UMAP 2", color = "") &
  theme_half_open()& 
  theme(legend.position = c(0.85,0.95), legend.key.size = unit(2, 'mm')) &
  theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          text = element_text(size=7), axis.ticks = element_blank(),
          axis.text=element_text(size=7),
          plot.title = element_text(size=7, face = "bold",hjust = 0.5))

seu.live <- RunUMAP(seu.live,dims = 1:20, reduction = "pcanoIGRNA", reduction.name = "pcanoIGRNA_UMAP" )



p.clusters.umap.live_noIg_RNA <-DimPlot(seu.live, reduction = 'pcanoIGRNA_UMAP',group.by = "clusters_pcaIG_named", label = TRUE, repel = TRUE, label.size = 2.5,pt.size =  0.5, cols = c("#009E73", "#D55E00","#F0E442","grey")
)   &
  labs(x = "UMAP 1", y = "UMAP 2", color = "Clusters", title = "Live cells UMAP based on \nRNA modality PCA \nwithout Ig-genes") &
  theme_half_open()& 
  theme(legend.position = "none", legend.key.size = unit(1, 'mm')) &
  theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          text = element_text(size=7), axis.ticks = element_blank(),
          axis.text=element_text(size=7),
          plot.title = element_text(size=6, face = "bold"))

p.umap.Ig.live_noIg_RNA <- FeaturePlot(seu.live, features = c("IGHM","IGHA1","IGHG1"),slot = "scale.data", reduction = 'pcanoIGRNA_UMAP', max.cutoff = 3, 
                  cols = c("lightgrey","dodgerblue3"), ncol = 3, pt.size =  0.5)  &
  labs(x = "UMAP 1", y = "UMAP 2", color = "") &
  theme_half_open()& 
  theme(legend.position = c(0.85,0.95), legend.key.size = unit(2, 'mm')) &
  theme(axis.text.x = element_blank(),
          axis.text.y = element_blank(),
          text = element_text(size=7), axis.ticks = element_blank(),
          axis.text=element_text(size=7),
          plot.title = element_text(size=7, face = "bold",hjust = 0.5))
# 
# seu.live <- FindMultiModalNeighbors(
#   seu.live, reduction.list = list("pcanoIGRNA", "pcanoIGPROT"), 
#   dims.list = list(1:10, 1:10), modality.weight.name = "RNA.weight"
# )
# 
# seu.live <- RunUMAP(seu.live, nn.name = "weighted.nn", reduction.name = "wnn.umap", reduction.key = "wnnUMAP_")
# seu.live <- FindClusters(seu.live, graph.name = "wsnn", algorithm = 3, resolution = 1.2, verbose = FALSE)
# 
# seu.live[["clusters_pcanoIG"]] <- Idents(object = seu.live)
# 
# 
# p.umap.Ig.prot.live2 <- FeaturePlot(seu.live, features = c("IgM","IgA","IgG"),slot = "scale.data", reduction = 'wnn.umap', max.cutoff = 2, 
#                   cols = c("lightgrey","deeppink3"), ncol = 3, pt.size = 0.5)  &
#   labs(x = "UMAP 1", y = "UMAP 2", color = "") &
#   theme_half_open()& 
#   theme(legend.position = c(0.85,0.95), legend.key.size = unit(2, 'mm')) &
#   theme(axis.text.x = element_blank(),
#           axis.text.y = element_blank(),
#           text = element_text(size=7), axis.ticks = element_blank(),
#           axis.text=element_text(size=7),
#           plot.title = element_text(size=7, face = "bold",hjust = 0.5))
# 
# p.umap.Ig.RNA.live2 <- FeaturePlot(seu.live, features = c("IGHM","IGHA1","IGHG1"),slot = "scale.data", reduction = 'wnn.umap', max.cutoff = 3, 
#                   cols = c("lightgrey","dodgerblue3"), ncol = 3, pt.size =  0.5)  &
#   labs(x = "UMAP 1", y = "UMAP 2", color = "") &
#   theme_half_open()& 
#   theme(legend.position = c(0.85,0.95), legend.key.size = unit(2, 'mm')) &
#   theme(axis.text.x = element_blank(),
#           axis.text.y = element_blank(),
#           text = element_text(size=7), axis.ticks = element_blank(),
#           axis.text=element_text(size=7),
#           plot.title = element_text(size=7, face = "bold",hjust = 0.5))
#           
# p.clusters.umap.live2 <- DimPlot(seu.live, reduction = 'wnn.umap',group.by = "clusters_pcaIG_named", label = TRUE, repel = TRUE, label.size = 2.5,pt.size =  0.5, cols = c("#009E73", "#D55E00","#F0E442","grey")
# )   &
#   labs(x = "UMAP 1", y = "UMAP 2", color = "Clusters", title = "Live cells WNN & UMAP \nwithout Ig-proteins or genes") &
#   theme_half_open()& 
#   theme(legend.position = "none", legend.key.size = unit(1, 'mm')) &
#   theme(axis.text.x = element_blank(),
#           axis.text.y = element_blank(),
#           text = element_text(size=7), axis.ticks = element_blank(),
#           axis.text=element_text(size=7),
#           plot.title = element_text(size=6, face = "bold"))
# 
# p.donor.umap.live2 <- DimPlot(seu.live, reduction = 'wnn.umap', group.by = "donor",label = F, pt.size =  0.5,repel = TRUE, label.size = 2.5, cols = colors.donors)   &
#   labs(x = "UMAP 1", y = "UMAP 2", color = "Donor", title = "No-Ig clustering overlap \nacross 3 donors") &
#   theme_half_open()& 
#   theme(legend.position = c(0.85, 0.95), legend.key.size = unit(1, 'mm'), legend.text =element_text(size=6), legend.title =element_text(size=7)  ) &
#   theme(axis.text.x = element_blank(),
#           axis.text.y = element_blank(),
#           text = element_text(size=7), axis.ticks = element_blank(),
#           axis.text=element_text(size=7),
#           plot.title = element_text(size=6, face = "bold"))
# 
# seu.live <- RenameIdents(object = seu.live, `1` = "IgM", `2` = "IgG", `3` = "IgA")
# seu.live[["clusters_pcaIG_named"]] <- Idents(object = seu.live)

fig.2.suppl.row1_2 <- plot_grid(p.clusters.umap.fix2,  p.umap.Ig.prot.fix2,labels = panellabels[c(1,3)], label_size = 10, rel_widths = c(1,3))

fig.2.suppl.row2_2 <- plot_grid(p.clusters.umap.live_noIg_RNA, p.umap.Ig.live_noIg_RNA,labels = panellabels[2], label_size = 10,
          rel_widths = c(1,3))

p.umap.fix.all_2 <- plot_grid(p.clusters.umap.live_noIg_PROT, 
          p.umap.Ig.live_noIg_PROT, labels = panellabels[c(4,5)], 
          label_size = 10,
          rel_widths = c(1,3))


fig.2.suppl.full_2 <- plot_grid(fig.2.suppl.row1_2, fig.2.suppl.row2_2,p.umap.fix.all_2, labels = "", label_size = 10, ncol = 1, rel_heights = c(1,1,1))

 ggsave(fig.2.suppl.full_2, filename = "output/paper_figures/Suppl_SeuratCluster_noIg.pdf", width = 177, height = 180, units = "mm",  dpi = 300, useDingbats = FALSE)
#  # ggsave(fig.2.suppl.full_2, filename = "output/paper_figures/Suppl_SeuratCluster_noIg.TIFF", width = 177, height = 180, units = "mm",  dpi = 300)
ggsave(fig.2.suppl.full_2, filename = "output/paper_figures/Suppl_SeuratCluster_noIg.png", width = 177, height = 180, units = "mm",  dpi = 300)

fig.2.suppl.full_2

Save dataset

Save seurat objects with selected non-cycling & CD27+IgD- cells.

#saveRDS(seu.fix, "output/seu.fix_norm_plasmacells.rds")
#saveRDS(seu.live, "output/seu.live_norm_plasmacells.rds")

#saveRDS(mofa, file= "output/MOFA_analysis_Donorgroup_clustered.rds")

seu.live.filtered.RNA <- seu.live
DefaultAssay(seu.live.filtered.RNA) <- "RNA"
seu.live.filtered.RNA[["SCT"]] <- NULL

seu.live.filtered.RNA[["PROT"]] <- NULL

#saveRDS(seu.live.filtered.RNA, "output/seu.live_norm_plasmacells_RNA.rds")

Session information

sessionInfo()

R version 4.0.3 (2020-10-10)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19044)

Matrix products: default

locale:
[1] LC_COLLATE=English_Netherlands.1252  LC_CTYPE=English_Netherlands.1252   
[3] LC_MONETARY=English_Netherlands.1252 LC_NUMERIC=C                        
[5] LC_TIME=English_Netherlands.1252    

attached base packages:
[1] parallel  stats4    stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] plyr_1.8.6           umap_0.2.9.0         MOFA2_1.1.17        
 [4] pzfx_0.3.0           ggupset_0.3.0        RColorBrewer_1.1-2  
 [7] enrichplot_1.10.2    UCell_1.0.0          data.table_1.14.2   
[10] scales_1.1.1         cowplot_1.1.1        ggthemes_4.2.4      
[13] kableExtra_1.3.4     knitr_1.36           org.Hs.eg.db_3.12.0 
[16] AnnotationDbi_1.52.0 IRanges_2.24.1       S4Vectors_0.28.1    
[19] Biobase_2.50.0       BiocGenerics_0.36.1  forcats_0.5.1       
[22] stringr_1.4.0        dplyr_1.0.7          purrr_0.3.4         
[25] readr_2.1.0          tidyr_1.1.4          tibble_3.1.5        
[28] ggplot2_3.3.5        tidyverse_1.3.1      Matrix_1.3-4        
[31] SeuratObject_4.0.2   Seurat_4.0.2         workflowr_1.6.2     

loaded via a namespace (and not attached):
  [1] utf8_1.2.2            reticulate_1.22       tidyselect_1.1.1     
  [4] RSQLite_2.2.8         htmlwidgets_1.5.4     BiocParallel_1.24.1  
  [7] grid_4.0.3            Rtsne_0.15            scatterpie_0.1.7     
 [10] munsell_0.5.0         ragg_1.2.0            codetools_0.2-16     
 [13] ica_1.0-2             future_1.23.0         miniUI_0.1.1.1       
 [16] withr_2.4.3           colorspace_2.0-2      GOSemSim_2.16.1      
 [19] filelock_1.0.2        highr_0.9             rstudioapi_0.13      
 [22] ROCR_1.0-11           tensor_1.5            DOSE_3.16.0          
 [25] listenv_0.8.0         MatrixGenerics_1.2.1  labeling_0.4.2       
 [28] git2r_0.28.0          polyclip_1.10-0       pheatmap_1.0.12      
 [31] bit64_4.0.5           farver_2.1.0          rhdf5_2.34.0         
 [34] rprojroot_2.0.2       basilisk_1.2.1        parallelly_1.29.0    
 [37] vctrs_0.3.8           generics_0.1.1        xfun_0.26            
 [40] R6_2.5.1              graphlayouts_0.7.2    rhdf5filters_1.2.1   
 [43] DelayedArray_0.16.3   fgsea_1.16.0          spatstat.utils_2.2-0 
 [46] cachem_1.0.6          assertthat_0.2.1      promises_1.2.0.1     
 [49] ggraph_2.0.5          gtable_0.3.0          globals_0.14.0       
 [52] goftest_1.2-2         tidygraph_1.2.0       rlang_0.4.11         
 [55] systemfonts_1.0.3     splines_4.0.3         lazyeval_0.2.2       
 [58] spatstat.geom_2.2-2   broom_0.7.10          yaml_2.2.1           
 [61] reshape2_1.4.4        abind_1.4-5           modelr_0.1.8         
 [64] backports_1.3.0       httpuv_1.6.3          qvalue_2.22.0        
 [67] tools_4.0.3           ellipsis_0.3.2        spatstat.core_2.3-0  
 [70] jquerylib_0.1.4       ggridges_0.5.3        Rcpp_1.0.7           
 [73] basilisk.utils_1.2.2  openssl_1.4.5         rpart_4.1-15         
 [76] deldir_1.0-2          pbapply_1.5-0         viridis_0.6.2        
 [79] zoo_1.8-9             haven_2.4.3           ggrepel_0.9.1        
 [82] cluster_2.1.0         fs_1.5.0              magrittr_2.0.1       
 [85] RSpectra_0.16-0       scattermore_0.7       DO.db_2.9            
 [88] lmtest_0.9-38         reprex_2.0.1          RANN_2.6.1           
 [91] whisker_0.4           fitdistrplus_1.1-6    matrixStats_0.61.0   
 [94] hms_1.1.1             patchwork_1.1.1       mime_0.12            
 [97] evaluate_0.14         xtable_1.8-4          readxl_1.3.1         
[100] gridExtra_2.3         compiler_4.0.3        shadowtext_0.0.9     
[103] KernSmooth_2.23-17    crayon_1.4.2          htmltools_0.5.2      
[106] ggfun_0.0.4           mgcv_1.8-33           later_1.3.0          
[109] tzdb_0.2.0            lubridate_1.8.0       DBI_1.1.1            
[112] corrplot_0.92         tweenr_1.0.2          dbplyr_2.1.1         
[115] rappdirs_0.3.3        MASS_7.3-53           cli_3.1.0            
[118] igraph_1.2.6          pkgconfig_2.0.3       plotly_4.10.0        
[121] spatstat.sparse_2.0-0 xml2_1.3.2            svglite_2.0.0        
[124] bslib_0.3.1           webshot_0.5.2         rvest_1.0.2          
[127] digest_0.6.28         sctransform_0.3.2     RcppAnnoy_0.0.19     
[130] spatstat.data_2.1-0   fastmatch_1.1-3       rmarkdown_2.11       
[133] cellranger_1.1.0      leiden_0.3.9          uwot_0.1.10          
[136] shiny_1.7.1           lifecycle_1.0.1       nlme_3.1-149         
[139] jsonlite_1.7.2        Rhdf5lib_1.12.1       askpass_1.1          
[142] limma_3.46.0          viridisLite_0.4.0     fansi_0.5.0          
[145] pillar_1.6.4          lattice_0.20-41       fastmap_1.1.0        
[148] httr_1.4.2            survival_3.2-7        GO.db_3.12.1         
[151] glue_1.4.2            png_0.1-7             bit_4.0.4            
[154] HDF5Array_1.18.1      ggforce_0.3.3         stringi_1.7.5        
[157] sass_0.4.0            blob_1.2.2            textshaping_0.3.6    
[160] memoise_2.0.1         irlba_2.3.3           future.apply_1.8.1