Introduction

Here we will be looking at gene expression in B-cells. Here are the sample descriptions:

  • FLW3: healthy no covid

  • FLW5: healthy no covid

  • 1Recov3M: had covid post-3 months

  • 2Recov3M: had covid, same person,

  • 2Recov12M: had covid, same person 9 months later

  • MDC8: had covid (sample ~3 month post, was treated with steroids)

  • MMC7: had covid (sample ~3 month post, was treated with steroids).

suppressPackageStartupMessages({
  library("dplyr")
  library("ggplot2")
  library("plyr")
  library("Seurat")
  library("hdf5r")
  library("SingleCellExperiment")
  library("parallel")
  library("stringi")
  library("beeswarm")
  library("muscat")
  library("DESeq2")
  library("mitch")
  library("harmony")
  library("celldex")
  library("SingleR")
  library("kableExtra")
  library("gplots")
  library("eulerr")
})

Load data

Load the h5 matrices.

flw3 <- Read10X_h5("FLW3/outs/filtered_feature_bc_matrix.h5")
## Genome matrix has multiple modalities, returning a list of matrices for this genome
flw3 <- flw3[[1]]
colnames(flw3) <- paste("flw3",colnames(flw3))
flw3_metadata <- data.frame(colnames(flw3))

flw5 <-  Read10X_h5("FLW5/outs/filtered_feature_bc_matrix.h5")
## Genome matrix has multiple modalities, returning a list of matrices for this genome
flw5 <- flw5[[1]]
colnames(flw5) <- paste("flw5",colnames(flw5))
flw5_metadata <- data.frame(colnames(flw5))

X1Recov3M <- Read10X_h5("1Recov3M/outs/filtered_feature_bc_matrix.h5")
## Genome matrix has multiple modalities, returning a list of matrices for this genome
X1Recov3M <- X1Recov3M[[1]]
colnames(X1Recov3M) <- paste("X1Recov3M",colnames(X1Recov3M))
X1Recov3M_metadata <- data.frame(colnames(X1Recov3M))

X2Recov12M <- Read10X_h5("2Recov12M/outs/filtered_feature_bc_matrix.h5")
## Genome matrix has multiple modalities, returning a list of matrices for this genome
X2Recov12M <- X2Recov12M[[1]]
colnames(X2Recov12M) <- paste("X2Recov12M",colnames(X2Recov12M))
X2Recov12M_metadata <- data.frame(colnames(X2Recov12M))

X2Recov3M <- Read10X_h5("2Recov3M/outs/filtered_feature_bc_matrix.h5")
## Genome matrix has multiple modalities, returning a list of matrices for this genome
X2Recov3M <- X2Recov3M[[1]]
colnames(X2Recov3M) <- paste("X2Recov3M",colnames(X2Recov3M))
X2Recov3M_metadata <- data.frame(colnames(X2Recov3M))

MDC8 <- Read10X_h5("MDC8/outs/filtered_feature_bc_matrix.h5")
## Genome matrix has multiple modalities, returning a list of matrices for this genome
MDC8 <- MDC8[[1]]
colnames(MDC8) <- paste("MDC8",colnames(MDC8))
MDC8_metadata <- data.frame(colnames(MDC8))

MMC7 <- Read10X_h5("MMC7/outs/filtered_feature_bc_matrix.h5")
## Genome matrix has multiple modalities, returning a list of matrices for this genome
MMC7 <- MMC7[[1]]
colnames(MMC7) <- paste("MMC7",colnames(MMC7))
MMC7_metadata <- data.frame(colnames(MMC7))

Analyse the data for quality. Look at the depth of quantification per cell.

message("FLW3")
## FLW3
dim(flw3)
## [1] 36611  8438
summary(colSums(flw3))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##      49    1454    2184    3913    3830  103443
sum(flw3)
## [1] 33020506
message("FLW5")
## FLW5
dim(flw5)
## [1] 36611  3844
summary(colSums(flw5))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##      61    1822    2650    5628    4412  854439
sum(flw5)
## [1] 21634457
message("X1Recov3M")
## X1Recov3M
dim(X1Recov3M)
## [1] 36611  6624
summary(colSums(X1Recov3M))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##      11    1020    1642    2840    3082  105824
sum(X1Recov3M)
## [1] 18809738
message("X2Recov12M")
## X2Recov12M
dim(X2Recov12M)
## [1] 36611  7454
summary(colSums(X2Recov12M))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##      40    1545    2166    3898    3609  304679
sum(X2Recov12M)
## [1] 29057677
message("X2Recov3M")
## X2Recov3M
dim(X2Recov3M)
## [1] 36611  2819
summary(colSums(X2Recov3M))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##      44    2040    3073    8079    4905  338583
sum(X2Recov3M)
## [1] 22775187
message("MDC8")
## MDC8
dim(MDC8)
## [1] 36611 11020
summary(colSums(MDC8))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##      45    1481    1938    2656    3125   78241
sum(MDC8)
## [1] 29269764
message("MMC7")
## MMC7
dim(MMC7)
## [1] 36611  8928
summary(colSums(MMC7))
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##      73    1477    2266    3559    3759  107588
sum(MMC7)
## [1] 31777959

FLW3 has 8k cells and FLW5 has 4k. Median seq depth is also significantly higher for FLW5. FLW1 has 50% more reads as compared to FLW5.

Filter to remove cells with fewer than 100 reads.

flw3 <- flw3[,colSums(flw3)>=100]
flw5 <- flw5[,colSums(flw5)>=100]
X1Recov3M <- X1Recov3M[,colSums(X1Recov3M)>=100]
X2Recov12M <- X2Recov12M[,colSums(X2Recov12M)>=100]
X2Recov3M <- X2Recov3M[,colSums(X2Recov3M)>=100]
MDC8 <- MDC8[,colSums(MDC8)>=100]
MMC7 <- MMC7[,colSums(MMC7)>=100]

Rename cells by sample origin, and join.

comb <- cbind(flw3,flw5,X1Recov3M,X2Recov3M,X2Recov12M,MDC8,MMC7)

cellmetadata <- as.data.frame(colnames(comb))
colnames(cellmetadata) = "cell_id"
cellmetadata$library <-sapply(strsplit(cellmetadata[,1]," "),"[[",1)

comb <- CreateSeuratObject(counts = comb, project = "bcell", min.cells = 5, meta.data = cellmetadata)
## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')

Normalise scRNA-seq data.

comb <- NormalizeData(comb)
## Normalizing layer: counts
comb <- FindVariableFeatures(comb, selection.method = "vst", nfeatures = 2000)
## Finding variable features for layer counts
comb <- ScaleData(comb)
## Centering and scaling data matrix
comb <- RunPCA(comb, features = VariableFeatures(object = comb), npcs = 20)
## PC_ 1 
## Positive:  LYZ, CD36, PID1, HCK, RTN1, MS4A6A, CPPED1, IFI30, FTH1, TMTC2 
##     TCF7L2, G0S2, MEGF9, DYSF, LINC02432, S100Z, PLXNB2, NRG1, HK3, MS4A7 
##     PDGFC, CD14, ZNF804A, PECAM1, MYOF, FAM198B-AS1, CXCL8, TMTC1, SERPINA1, LPAR1 
## Negative:  CAMK4, LEF1, BACH2, BCL11B, INPP4B, IL7R, SKAP1, BCL2, ANK3, TC2N 
##     NELL2, PCED1B, THEMIS, NR3C2, CD96, ZEB1, TCF7, TSHZ2, TXK, CCND3 
##     PDE7A, ITK, MLLT3, ABLIM1, RIPOR2, LTB, RPL3, PRKCA, HIVEP2, RPL13 
## PC_ 2 
## Positive:  BANK1, AFF3, FCRL1, RALGPS2, MS4A1, PAX5, EBF1, MEF2C, IGHM, LINC00926 
##     COL19A1, OSBPL10, CD74, MARCH1, KHDRBS2, CD79A, BLK, HLA-DRA, ADAM28, COBLL1 
##     ARHGAP24, BCL11A, NIBAN3, HDAC9, CD22, IGHD, FCHSD2, PLEKHG1, CCSER1, GNG7 
## Negative:  B2M, BCL11B, IL32, ITK, TMSB4X, THEMIS, SKAP1, PRKCQ, SLFN12L, HLA-B 
##     INPP4B, CD3D, TRAC, CTSW, TC2N, CD96, HLA-C, HLA-A, CAMK4, IL7R 
##     PITPNC1, PDE3B, TRBC1, PFN1, CD6, SPOCK2, TRBC2, LINC00861, MGAT4A, CD7 
## PC_ 3 
## Positive:  B2M, TMSB4X, HLA-B, HLA-C, PFN1, CD79A, HLA-A, BANK1, MS4A1, RPS2 
##     TMSB10, IGHM, RPL41, FCRL1, PAX5, RPLP1, SH3BGRL3, EBF1, CFL1, RPL13A 
##     RPL3, CTSW, FAU, LINC00926, CD79B, COL19A1, RPS11, OSBPL10, RPL13, RPS15 
## Negative:  CD36, DISC1, FYB1, XIST, TMTC2, RTN1, PID1, CPPED1, LDLRAD4, JAML 
##     MS4A6A, SERINC5, PRKCA, MAML2, LEF1, MEGF9, HCK, FHIT, PLCL1, PECAM1 
##     FAM13A, S100Z, MBNL1, ZNF609, FOXP1, LYZ, IGF1R, NELL2, NRG1, LINC02432 
## PC_ 4 
## Positive:  BACH2, LEF1, CAMK4, ANK3, NELL2, INPP4B, BCL2, TSHZ2, SKAP1, NR3C2 
##     CCR7, BCL11B, IL7R, TC2N, ZEB1, ABLIM1, PATJ, DOCK9, PDE7A, TRABD2A 
##     TCF7, COL19A1, TESPA1, CD28, HIVEP2, CSGALNACT1, SATB1-AS1, PCED1B, RASGRF2, MLLT3 
## Negative:  CST3, ACTB, RTN1, LYZ, CPPED1, FCER1G, FTH1, IFI30, TCF7L2, TYROBP 
##     S100A6, HCK, PECAM1, ACTG1, FLNA, LST1, CSF1R, CD36, PFN1, SH3BGRL3 
##     PID1, MYOF, HDAC9, CFL1, FTL, MS4A7, ARPC1B, OAZ1, GAPDH, FCGR3A 
## PC_ 5 
## Positive:  TCF7L2, HCK, POU2F2, MS4A7, AC105402.3, CSF1R, CPPED1, LINC02432, RTN1, S100Z 
##     SPRED1, LST1, MARCH1, IFI30, MYOF, MEGF9, TMTC2, PID1, AC104809.2, TMTC1 
##     LYZ, FCGR3A, PECAM1, SSH2, HK3, AC009093.2, AC092546.1, EPS8, SERPINA1, CDKN1C 
## Negative:  LINC01478, TUBB1, LINC01374, CAVIN2, GP1BB, ITGA2B, ITGB3, CUX2, PTPRS, CLU 
##     RHEX, EPHB1, FAM160A1, AC023590.1, LINC00996, NRP1, SPARC, PPBP, SCN9A, NRGN 
##     COL24A1, PTCRA, PF4, PLXNA4, GP9, GNG11, MPIG6B, TREML1, SCN1A-AS1, PHEX
comb <- RunHarmony(comb,"library")
## Transposing data matrix
## Initializing state using k-means centroids initialization
## Harmony 1/10
## Harmony 2/10
## Harmony 3/10
## Harmony 4/10
## Harmony 5/10
## Harmony converged after 5 iterations
DimHeatmap(comb, dims = 1:6, cells = 500, balanced = TRUE)

ElbowPlot(comb)

#comb <- JackStraw(comb, num.replicate = 100)
comb <- FindNeighbors(comb, dims = 1:5)
## Computing nearest neighbor graph
## Computing SNN
comb <- FindClusters(comb, algorithm = 3, resolution = 0.35, verbose = FALSE)
comb <- RunUMAP(comb, dims = 1:5)
## 00:09:17 UMAP embedding parameters a = 0.9922 b = 1.112
## 00:09:17 Read 48988 rows and found 5 numeric columns
## 00:09:17 Using Annoy for neighbor search, n_neighbors = 30
## 00:09:17 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 00:09:21 Writing NN index file to temp file /tmp/RtmpjywqML/file27ed8f65cbc7ac
## 00:09:21 Searching Annoy index using 1 thread, search_k = 3000
## 00:09:35 Annoy recall = 100%
## 00:09:36 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 00:09:38 Initializing from normalized Laplacian + noise (using RSpectra)
## 00:09:39 Commencing optimization for 200 epochs, with 1782842 positive edges
## 00:09:52 Optimization finished
DimPlot(comb, reduction = "umap")

#comb <- FindNeighbors(comb, dims = 1:10)
#comb <- FindClusters(comb, algorithm = 3, resolution = 0.5, verbose = FALSE)
#comb <- RunUMAP(comb, dims = 1:10)
#DimPlot(comb, reduction = "umap")

Look at the prominent cell markers.

message("Naive CD4+ T") # 1,3,5,6,7
## Naive CD4+ T
VlnPlot(comb, features = c("IL7R", "CCR7"))

message("CD14+ Mono")
## CD14+ Mono
VlnPlot(comb, features = c("CD14", "LYZ"))

message("Memory CD4+")
## Memory CD4+
VlnPlot(comb, features = c("IL7R", "S100A4"))

message("B") # 2,7
## B
VlnPlot(comb, features = c("MS4A1"))

message("CD8+ T") #2
## CD8+ T
VlnPlot(comb, features = c("CD8A"))

message("FCGR3A+ Mono") #?
## FCGR3A+ Mono
VlnPlot(comb, features = c("FCGR3A", "MS4A7"))

message("NK") # 2,5,6,7?
## NK
VlnPlot(comb, features = c("GNLY", "NKG7"))

message("DC") # 7?
## DC
VlnPlot(comb, features = c("FCER1A", "CST3"))

message("Platelet") #6
## Platelet
VlnPlot(comb, features = c("PPBP"))

new.cluster.ids <- c("Naive CD4 T", "CD14+ Mono", "Memory CD4 T", "B", "CD8 T", "FCGR3A+ Mono",
    "NK", "DC", "Platelet")

FeaturePlot(comb, features = c("IL7R","CCR7","CD14","LYZ","S100A4","MS4A1", "CD8A" ,
  "FCGR3A", "MS4A7", "GNLY", "NKG7", "FCER1A", "CST3", "PPBP","CD3E" , "CD19"))

# b cell markers
FeaturePlot(comb, features = c("CD19","CD27","CD38","CD24"))

FeaturePlot(comb, features = c("CR2", "CD34", "MME", "MS4A1"))

FeaturePlot(comb, features = c("MZB1", "CXCR3", "FCRL5", "TBX21"))

FeaturePlot(comb, features = c("CCR6", "ITGAX","MX1","BST2"))

Cell Type Annotation

From Turner Lab: Create a reference using the monaco immune database. The Monaco reference consists of bulk RNA-seq samples of sorted immune cell populations from GSE107011 (Monaco et al. 2019). (from the CellDex documentation)

ref <- celldex::MonacoImmuneData()

DefaultAssay(comb) <- "RNA"
comb2 <- as.SingleCellExperiment(comb)

lc <- logcounts(comb2)

pred_imm_broad <- SingleR(test=comb2, ref=ref,
                          labels=ref$label.main)

head(pred_imm_broad)
## DataFrame with 6 rows and 4 columns
##                                                    scores       labels
##                                                  <matrix>  <character>
## flw3 AAACAGCCACATTAAC-1 0.1838240:0.1798454:0.0668172:...    Monocytes
## flw3 AAACATGCACCGGTAT-1 0.2133065:0.1682427:0.3923648:... CD4+ T cells
## flw3 AAACATGCATAATGTC-1 0.1241514:0.1089673:0.2266038:... CD4+ T cells
## flw3 AAACATGCATACTCCT-1 0.2994710:0.1875276:0.3675424:...      T cells
## flw3 AAACATGCATTATGGT-1 0.0908606:0.0467543:0.1950370:... CD8+ T cells
## flw3 AAACCAACAATTAAGG-1 0.1425373:0.1192283:0.2563216:...      T cells
##                         delta.next pruned.labels
##                          <numeric>   <character>
## flw3 AAACAGCCACATTAAC-1  0.0883235     Monocytes
## flw3 AAACATGCACCGGTAT-1  0.0653410  CD4+ T cells
## flw3 AAACATGCATAATGTC-1  0.0272195  CD4+ T cells
## flw3 AAACATGCATACTCCT-1  0.0302115       T cells
## flw3 AAACATGCATTATGGT-1  0.0705356  CD8+ T cells
## flw3 AAACCAACAATTAAGG-1  0.0359376       T cells
table(pred_imm_broad$pruned.labels)
## 
##         B cells       Basophils    CD4+ T cells    CD8+ T cells Dendritic cells 
##            5497              45            5624            1409            2576 
##       Monocytes     Neutrophils        NK cells     Progenitors         T cells 
##           14472              29            9120             240            9425
cellmetadata$label <- pred_imm_broad$pruned.labels

pred_imm_fine <- SingleR(test=comb2, ref=ref,
                          labels=ref$label.fine)
head(pred_imm_fine)
## DataFrame with 6 rows and 4 columns
##                                                    scores
##                                                  <matrix>
## flw3 AAACAGCCACATTAAC-1 0.0950415:0.2862279:0.0972819:...
## flw3 AAACATGCACCGGTAT-1 0.3375173:0.1610507:0.2956709:...
## flw3 AAACATGCATAATGTC-1 0.1923652:0.0942964:0.1856419:...
## flw3 AAACATGCATACTCCT-1 0.3508795:0.1867607:0.3342319:...
## flw3 AAACATGCATTATGGT-1 0.1960738:0.0724535:0.1618841:...
## flw3 AAACCAACAATTAAGG-1 0.2261702:0.1108775:0.2228686:...
##                                         labels delta.next
##                                    <character>  <numeric>
## flw3 AAACAGCCACATTAAC-1    Classical monocytes 0.06841161
## flw3 AAACATGCACCGGTAT-1      Naive CD4 T cells 0.00749510
## flw3 AAACATGCATAATGTC-1              Th1 cells 0.00723785
## flw3 AAACATGCATACTCCT-1 Central memory CD8 T.. 0.00731102
## flw3 AAACATGCATTATGGT-1      Naive CD8 T cells 0.08295119
## flw3 AAACCAACAATTAAGG-1         Vd2 gd T cells 0.00329823
##                                  pruned.labels
##                                    <character>
## flw3 AAACAGCCACATTAAC-1    Classical monocytes
## flw3 AAACATGCACCGGTAT-1      Naive CD4 T cells
## flw3 AAACATGCATAATGTC-1              Th1 cells
## flw3 AAACATGCATACTCCT-1 Central memory CD8 T..
## flw3 AAACATGCATTATGGT-1      Naive CD8 T cells
## flw3 AAACCAACAATTAAGG-1         Vd2 gd T cells
table(pred_imm_fine$pruned.labels)
## 
##    Central memory CD8 T cells           Classical monocytes 
##                           648                          8946 
##   Effector memory CD8 T cells             Exhausted B cells 
##                           262                           319 
##     Follicular helper T cells        Intermediate monocytes 
##                           582                          4734 
##         Low-density basophils       Low-density neutrophils 
##                            48                            12 
##                    MAIT cells       Myeloid dendritic cells 
##                           918                          1710 
##                 Naive B cells             Naive CD4 T cells 
##                          3936                          2873 
##             Naive CD8 T cells          Natural killer cells 
##                          1163                          8056 
##       Non classical monocytes   Non-switched memory B cells 
##                          1413                          1024 
##            Non-Vd2 gd T cells                  Plasmablasts 
##                          4084                            84 
##  Plasmacytoid dendritic cells              Progenitor cells 
##                           280                           225 
##       Switched memory B cells            T regulatory cells 
##                           832                           788 
## Terminal effector CD4 T cells Terminal effector CD8 T cells 
##                           849                           697 
##                     Th1 cells                Th1/Th17 cells 
##                          1146                           672 
##                    Th17 cells                     Th2 cells 
##                           550                           536 
##                Vd2 gd T cells 
##                          1300
cellmetadata$finelabel <- pred_imm_fine$pruned.labels

col_pal <- c('#e31a1c', '#ff7f00', "#999900", '#cc00ff', '#1f78b4', '#fdbf6f',
             '#33a02c', '#fb9a99', "#a6cee3", "#cc6699", "#b2df8a", "#99004d", "#66ff99",
             "#669999", "#006600", "#9966ff", "#cc9900", "#e6ccff", "#3399ff", "#ff66cc",
             "#ffcc66", "#003399")

annot_df <- data.frame(
  barcodes = rownames(pred_imm_broad),
  monaco_broad_annotation = pred_imm_broad$labels,
  monaco_broad_pruned_labels = pred_imm_broad$pruned.labels,
  monaco_fine_annotation = pred_imm_fine$labels,
  monaco_fine_pruned_labels = pred_imm_fine$pruned.labels
)

meta_inf <- comb@meta.data
meta_inf$cell_barcode <- colnames(comb)

meta_inf <- meta_inf %>% dplyr::left_join(y = annot_df,
                                          by = c("cell_barcode" = "barcodes"))
rownames(meta_inf) <- colnames(lc)

comb@meta.data <- meta_inf

DimPlot(comb, label=TRUE, group.by = "monaco_broad_annotation", reduction = "umap",
  cols = col_pal, pt.size = 0.5) + ggtitle("Annotation With the Monaco Reference Database")

Extract B cells

bcells <- comb[,which(meta_inf$monaco_broad_annotation == "B cells")]
bcells_metainf <- meta_inf[which(meta_inf$monaco_broad_annotation == "B cells"),]

# remove non bcells
bcells_metainf1 <- bcells_metainf[grep("B cells",bcells_metainf$monaco_fine_pruned_labels),]
bcells_metainf2 <- bcells_metainf[grep("Plasmablasts",bcells_metainf$monaco_fine_pruned_labels),]
bcells_metainf <- rbind(bcells_metainf1,bcells_metainf2)

bcells <- bcells[,which(colnames(bcells) %in% rownames(bcells_metainf))]

bcells <- FindVariableFeatures(bcells, selection.method = "vst", nfeatures = 2000)
## Finding variable features for layer counts
bcells <- RunPCA(bcells, features = VariableFeatures(object = bcells))
## Warning: The following 64 features requested have zero variance; running
## reduction without them: MZB1, INPP4B, PDE3B, GLDC, PTCHD1-AS, TRABD2A, RASGRF2,
## SSR4, SGCD, GLIS3, DTHD1, HESX1, ADGRB3, TMEM232, DOCK9, TNFAIP3, VMP1, TRBC1,
## NR3C2, PTPRM, AC007100.1, SPARCL1, AL589740.1, GPC5, RAPGEF4-AS1, CCDC149,
## GZMK, PCAT4, MARVELD3, LINC01146, RASGRP1, BAALC-AS1, AC089985.1, ADARB2,
## NMNAT3, C15orf62, NEXMIF, AL049828.1, LINC01937, LINC01182, ANKRD34B, CNTLN,
## SND1, AC090844.2, LINC00159, CSF1R, F11-AS1, NIPAL3, AC009120.4, AC004895.1,
## EPHX3, PNCK, AC114956.1, CORT, ZDHHC11B, EFNB2, AC008474.1, GINS1, PRKCA-AS1,
## TLR5, RNF213, ZFYVE9, ABCC4, MAP3K13
## Warning in PrepDR5(object = object, features = features, layer = layer, : The
## following features were not available: VCAN, PLXDC2, MAN1A1, NAMPT, FNDC3B,
## HSPA5, LRMDA, XBP1, PLCB1, AOAH, PRKCH, SLC8A1, LINC01505, CD247, PRDM1,
## SDF2L1, NUSAP1, RBM47, AAK1, NEAT1, DPYD, CENPE, DCC, ACSL1, GPC6, TNIK, MANF,
## NCALD, PCSK5, FKBP11, STAT4, TOX, RORA, PPP2R2B, DMXL2, SLCO3A1, MYDGF, SEC24D,
## GNAQ, SNX32, CDC14A, ITGB1, SHCBP1, CEP128, CSF3R, IQGAP2, ARHGAP26, RRBP1,
## CALR, VCAN-AS1, GTSE1, PDIA4, CREB5, DOCK4, SLC11A1, A2M, STXBP6, SLAMF7,
## PDE4D, TFEC, RAB39A, AC008574.1, GAB1, AC079336.5, SSR3, MCTP1, MYBL1, KIF4A,
## AGAP1, AC243829.2, AL442163.1, AL138716.1, HIST1H4C, ATXN1, WWC1, PLAAT2, SMC4,
## SPATS2, IPCEF1, LINC02384, MAML3, EIF4E3, AC044781.1, SEC11C, LINC02296, TPRG1,
## MAGI1, ERN1, MCM10, NLGN4Y, CADPS2, FYN, LCP2, TOX2, SAMHD1, LMAN1, LINC02456,
## GPRIN3, PDIA6, SORL1, LINC02363, LINC01821, GLT1D1, CLEC12A, PPP1R26-AS1,
## SNAP25-AS1, ESRRB, WDFY3, IRAK3, TRERF1, BX322234.1, IL6ST, SLC22A15, PTPRE,
## SRGN, AC002072.1, LRRC8C, ARLNC1, CHODL-AS1, ME1, EPB41L4A, AC091938.1,
## COL12A1, AL807761.4, AC037198.1, ZNF239, LRRC4C, PWRN1, CCR9, NDFIP1, ROR2,
## IGF2BP3, ANO5, DENND1B, TYMS, BTBD11, RCAN2, WNK3, RGS13, AL356108.1, SEMA3B,
## AC026353.1, AL442636.1, CFAP54, NFIA, ESR1, AC073332.1, RIN2, FKBP5,
## AC006581.2, FBXO5, INPP4A, EXT1, MAFTRR, ELL2, RUNX1, CERS6, KNL1, MIR181A1HG,
## MSC-AS1, NIBAN1, SEMA3C, MYO1F, ADAMTSL4-AS1, AC108169.1, AADACL2-AS1,
## MAMDC2-AS1, SPART-AS1, SAMD3, SYNE2, DAPK1, ESCO2, AL078602.1, AC011995.2,
## SPC25, Z99758.1, TCF7L1, LEPROTL1, IL20RB-AS1, LDLRAD3, CENPP, AL358777.1,
## DOCK5, DEPDC1B, MAGI3, SLFN5, HIVEP3, TAFA2, TRPA1, ANXA10, LINC02211, KLF12,
## SPN, LINC00334, ATF5, TGFBR3, PIK3R5, AP001977.1, KRT1, ADCY9, SPANXA2-OT1,
## IGF2BP2, MPP7, KIFC1, PPARG, CRELD2, VWC2, TLL2, AC079352.1, AC009974.2,
## AC027702.1, AC108925.1, TNC, TRAV8-3, CDC45, SEPTIN4-AS1, AC109811.1, ART4,
## MOCOS, UPK1A-AS1, ZNF541, AC026979.4, AP001831.1, RBMS1, AC005355.1,
## AC024558.2, AC104964.3, ANXA1, TNFAIP2, FRY, ARHGAP10, AC009975.1, NCAPG,
## EEPD1, CCL5, SULT2B1, NUCB2, ARL17B, AC011416.3, NIM1K, PTPRJ, PROX1,
## AC092611.1, TRIP13, TCERG1L, AC084809.1, AC090376.1, AL646090.2, IGHV3-73,
## AC012459.1, NHS, AC073050.1, CCDC113, SVIL, ADGRE2, DIO1, LINC01847,
## AL451166.1, TNFRSF10C, CLEC19A, AC103987.2, CDH2, BMX, LHFPL1, LINC00907,
## OPRM1, FANCC, SPOCK1, TBXAS1, AL360093.1, AC019117.2, TIAM1, CD226, AC073114.1,
## E2F1, CNGA1, AC007262.2, ARHGAP6, MVB12B, ZEB2, CYP27A1, LINC01163, PRICKLE2,
## PF4V1, TGFB2, LINC01169, CDCA5, PDGFD, FRMD4B, SUB1, RFX2, AC025263.1, PARP8,
## TXNDC11, PERP, ZCCHC24, NLRP3, KCNN3, HIST1H2AJ, AC008170.1, C1orf189, FCN1,
## DENND2A, KLRG1, SYTL3, TLR2, CDHR3, STK3, POU6F2, GAB2, PLA2G4A, COL18A1,
## HIST1H3C, AC020978.2, LINC02080, AC011509.2, GATA1, C8orf88, C14orf178, NRG3,
## TSPAN9, TAF1A-AS1, MUC13, PCDHGA12, OR1L8, LINC02625, AC079601.1, AL121772.3,
## CCN5, AC024382.1, CABCOCO1, SFMBT2, TENM1, POF1B, SYTL2, BRIP1, SDSL, FER1L5,
## SULF1, AC018845.1, PYHIN1, RBL2, IGLV1-51, C8orf37-AS1, DPYD-AS1, FAM222A-AS1,
## AC073475.1, ASPH, SLC26A7, TRAV17, NPIPA8, VSTM1, IL6R, CACNB4, GAS7,
## AL096854.1, ATP8B4, CD8A, MRC2, STK32B, JAKMIP2-AS1, AQP9, RNF175, CCDC88A,
## CFAP44-AS1, AL445928.2, CTBP2, UPP2, LINC02137, SMYD3, LINC01605, AC135803.1,
## CLCA4, AL592402.1, FRZB, LINC01878, FAM124B, AC023503.1, LINC01983, ODAPH,
## LINC01511, AC010280.1, LINC02526, AC005682.2, VPS37D, AC105206.3, AL512634.1,
## AL450322.1, ZNF214, AC090099.1, AC092111.1, CLEC1A, AC078962.2, IL22,
## AL355516.1, C13orf42, AL358334.4, LINC00911, AC087639.2, MESP1, KCNJ12,
## AC015911.2, AC115100.1, SLC14A2-AS1, AC138470.2, HIF3A, LINC01620, KCNS1,
## AL162457.1, IGLV3-9, AL023574.1, GDPD2, SLC1A7, SMG7-AS1, MINDY4B, LINC01276,
## AL591468.1, AC007938.3, GHET1, CARD17, AC121338.1, PAH, PCDH9-AS1, LINC00433,
## MCF2L-AS1, DSCAS, LINC01919, AP000265.1, AC245140.3, SELENOS, AL109840.2,
## MICAL2, SBF2, SOBP, LINC02008, KIF13A, FGD4, CADPS, HIST1H2BB, NMT2, KIAA0825,
## TANC2, TSHZ3, BUB1, ICOS, AL353151.2, LINC01229, BEX5, AL139294.1, CDS1, FOLR3,
## UBE2J1, NHSL1, NLRP12, SLC44A3, STOX2, AP000593.4, AC074051.5, LINC02860, FGD6,
## BUB1B, AC105180.2, LRGUK, AC010609.1, C1orf21, AC104459.1, AC079760.2,
## AP003110.1, CCDC7, H2AFX, IL4I1, AC023480.1, TRAIP, LYPLAL1-DT, JAKMIP2, TRPS1,
## MCC, AC092164.1, ASTL, VOPP1, TNR, SCOC-AS1, BMPR1A, AC026341.1, CLDN11,
## AL024498.1, LINC02676, IGHV5-10-1, DLGAP1-AS5, ZNF534, AC002472.1, MN1,
## APOBEC3H, ASMT, ST8SIA6, WDR64, OSCAR, LINC02328, TEX41, AC003044.1, NEDD4,
## SLC44A1, TSPAN5, TENT5A, KIF23, CAV3, ALPL, NCAPG2, C16orf71, FBXO15, WNT2B,
## HIPK2, SPCS2, PCNX1, SNTB1, PAM, TVP23C, ZBTB16, MYO3A, CEP70, FANCI, L2HGDH,
## KCNMB2-AS1, RGL1, DNM3, SIPA1L2, AC005833.2, ZYG11A, VIPR2, AP003071.2,
## RTN4RL1, PTPN4, AGTPBP1, CRADD, AL365295.1, GASK1A, ALDH1L2, AUXG01000058.1,
## PIP4P2, GLRA2, AC109927.1, C8orf31, AL355303.1, ID2, CEP85L, AC125618.1, NEIL3,
## AL022724.3, AC019257.1, BASP1-AS1, AL354733.2, CTNNA3, MARCO, DACT3, CFAP61,
## AC007610.1, PPP4R4, CMTM2, LPCAT2, POLE2, ESYT2, ETS1, MIR4435-2HG, AC068643.1,
## LINC01933, PPM1L, PLEKHA5, CYTOR, VIM, CYSLTR2, AC066613.2, LMAN2, UTY, SCLT1,
## TNFSF8, AL359644.1, ERC1, AC015923.1, ARHGAP22, AC097460.1, CMTM4, CLEC7A,
## CPVL, CLIP4, SPAG16, AL365255.1, RGS1, ANKRD7, ATP8A2, PRAM1, TUBA1B,
## LGALSL-DT, LINC02234, CHN1, ZRANB2-AS2, TMEM252-DT, KLRC2, LINC00278, TRAT1,
## AC108718.1, DIAPH2-AS1, CIT, SGMS2, ZBP1, SEMA4D, C6orf58, HACD1, SCN2A, MCF2,
## MAP7, MLEC, FRMD3, ASAP3, AL512288.1, AC108210.1, ZNF391, PFKFB1, HLF, CD3G,
## ME3, CCSER2, HLCS, AC020916.1, LINC00298, AL357054.5, AL357141.1, TSPAN11,
## AC007611.1, AC100863.1, KLHL34, LINC00571, ISOC2, SGSM1, CPQ, MYO5B, PLAUR,
## GAS2, AL591518.1, DACH2, EVA1C, C20orf194, FOS, PLS1, SSC5D, AC107068.2,
## AL358937.1, PAK1, HYOU1, TIMP3, KIF2C, PRLR, HNRNPLL, TRG-AS1, TLN2, C4orf36,
## CERS6-AS1, RLN2, AL354989.1, CPEB1, DHRS3, NAV1, SLC24A4, MTHFD1L, SYNE1,
## AL357873.1, COL16A1, ERICH3-AS1, AC095030.1, LINC02609, AL590666.4, SERPINC1,
## TEX35, NEK2, WNT3A, STON1-GTF2A1L, LINC01793, AC016747.2, IGKV5-2, KLHL41,
## AC108066.1, GRM7-AS3, ZNF662, CADM2-AS2, IGSF10, AC092944.3, TRPC3, AC104083.1,
## AC026726.1, LINC01942, AL365226.1, CFAP206, TRDN-AS1, ITPRID1, TRGV9,
## AC005076.1, AC067930.1, AL445526.1, FAM201A, AL353742.1, AKR1C4, LINC00707,
## AL354916.1, ASAH2, PLEKHS1, AL139123.1, AC104383.2, CBLIF, NECTIN1-AS1, EMP1,
## GPD1, LINC02389, AC091214.1, AC089999.2, CCDC169, TNFSF11, DZIP1, AL358332.1,
## C15orf54, AC103740.1, CCDC33, AC068870.2, KIF7, AC022819.1, CERS3-AS1,
## AC090907.1, AC120498.4, AC012676.1, AC099518.5, AC015908.2, FBXW10, TTLL6,
## CACNG1, KIF19, AC091691.2, AC011446.1, AC008747.1, AC006213.1, AL354993.2,
## AP000355.1, APOL4, AL683807.2, DIPK2B, RAB40A, RNASEH2B-AS1, TMEM170B, HDAC4,
## ZEB1-AS1, LINC02798, OSTN-AS1, TTTY10, LRIG1, ARSG, ATP2B4, TMEM163, JAZF1-AS1,
## WWC2, NETO2, AC122697.1, DNAJC1, MANEA, WDR17, AC079174.2, AC008870.4,
## LINC01692, NR6A1, CA3, DST, SLC2A13, CCDC144A, P4HB, SLC5A9, PTGFR, AL391811.1,
## LINC01816, SULT1C4, LINC01173, CNTN4-AS1, THOC7-AS1, SPTSSB, LINC02505,
## AC096759.1, AC107067.2, AC116616.1, DNAH5, LINC02060, LINC01622, TFAP2A,
## POPDC3, LINC00326, AL139393.1, AC005100.1, PCOLCE-AS1, AC073314.1, LINC00681,
## CDH17, AL157829.1, NRARP, PITX3, TAC3, AVPR1A, AC010201.1, LINC00423, TRAV26-1,
## FLRT2, AL355836.2, IGHV4-28, AC090617.1, TMIGD1, AC022903.2, AP001029.1, MEP1B,
## CNN1, AC010463.3, DLL3, EHD2, LINC00945, BRWD1-AS2, LINC01424, AP001468.1,
## TCN2, ARSF, TLR8-AS1, FMR1NB, AL133334.1, AC004817.5, LINC02133, KLRK1, DGKG,
## PSD3, CD2, ARNTL2, KIF9-AS1, AC004158.1, LARP1B, GATA3, RAB31, SSR1, NUGGC,
## MIR646HG, LINC01625, AC025569.1, PPEF1, RFX8, MIR193BHG, MYBL2, AL596218.1,
## PSTPIP2, UACA, ATL1, ALG14, MS4A4E, AL035427.1, WAKMAR2, ITPKB, AC139769.2,
## LMNB1, SYBU, ARRB1, CLIC5, FRRS1, AC119674.1, BCAT1, NMRK1, MBNL1-AS1, RAB3C,
## AC018695.2, LINC02246, RAB20, AC097515.1, MAP4K3-DT, VCL, CCDC30, AMPH, MELK,
## AC021086.1, HIST1H2AL, UGGT2, SAT1, DDAH1, FSD1L, MCM4, RPS6KA6, AC104170.1,
## STK38, GCNT4, ANO10, PCBP3, TRDC, DENND2C, TK1, LINC01353, SLC9A2, RNF180,
## AL158801.2, AC136431.1, ITGA2, TECTA, KANK2, SMPDL3A, CDC6, ZSCAN31,
## AC096536.3, CRMP1, AMBN, EGR3, CLEC6A, TMPRSS12, AC122685.1, FERMT2,
## AC027458.1, PLAC1, SLC8A1-AS1, PXDNL, HM13, NCOA1, HDX, TACR1, HIST1H1B,
## NT5DC2, CDCA8, BRCA1, FBXO32, ZNF438, LINC01118, LINC00299, SAMSN1, CD59, Z
## PC_ 1 
## Positive:  FYB1, PITPNC1, BCL11B, PRKCA, ITK, THEMIS, PHACTR2, PRKCQ, SLFN12L, IL32 
##     TC2N, SERINC5, CAMK4, CD36, ATP10A, JAML, TXK, PLCL1, DISC1, MGAT4A 
##     IL7R, PAG1, PID1, LINC00861, CD6, RTN1, TCF7L2, CD96, TRAC, ARHGEF3 
## Negative:  IGHM, IGKC, STEAP1B, SOX5, AL355076.2, FCRL5, PCDH9, KLHL14, SLC38A11, IGLC2 
##     GPM6A, SSPN, AL592429.2, JCHAIN, RHEX, KCNQ5, IGLC3, AC083837.1, LINC01811, MACROD2 
##     IFNG-AS1, AL078459.1, AC007368.1, RGS7, PARM1, GEN1, NCOA3, AC108879.1, AC008014.1, IGLC1 
## PC_ 2 
## Positive:  AL355076.2, SSPN, IGHA2, SOX5, IGHA1, JCHAIN, RHEX, AC008691.1, CSGALNACT1, EML6 
##     DNAH8, MYO1D, TNFRSF17, TEX9, PVT1, AL592429.2, IGHG1, LINC01320, GEN1, TP63 
##     FUT8, IGHG3, AC007368.1, ITM2C, DERL3, XIST, TXNDC5, SAMD12, DOCK10, FA2H 
## Negative:  IGHM, PCDH9, AC108879.1, STEAP1B, SLC38A11, SYN3, AKAP6, THRB, LINC01811, NETO1 
##     AL139020.1, FOXP1, MACROD2, RAPGEF5, SDK1, KLHL14, LINC01340, LINC01374, PLD5, ACSM3 
##     CARMIL1, LINC02161, RIMBP2, TTC28, SKAP1, TTTY14, LINC01572, RGS7, LINC02550, CHL1 
## PC_ 3 
## Positive:  SLC38A11, STEAP1B, RGS7, AKAP6, AC108879.1, FOXP1, AC007368.1, LEF1, KLHL14, GPM6A 
##     IGHM, NETO1, CAMK4, MBNL1, BCL11B, BCL2, MAML2, PCDH9, IGF1R, SYN3 
##     SKAP1, PRKCA, CCR7, SERINC5, LINC02550, IL7R, MLLT3, PRKN, TC2N, HIVEP2 
## Negative:  SSPN, IGHA1, IGHA2, TXNDC5, EML6, AC008691.1, JCHAIN, SOX5, LINC01320, TNFRSF17 
##     TP63, DNAH8, IGHG1, IGHG3, TEX9, DERL3, IGHGP, AL355076.2, CRIP1, PPP1R9A 
##     IGHG2, KCNMA1, IGLC1, EYA2, ITM2C, IGHG4, FA2H, MYO1D, SOX5-AS1, ACOXL 
## PC_ 4 
## Positive:  TXNDC5, AC108879.1, IGHA1, IGHA2, JCHAIN, AL589693.1, LINC01811, IFNG-AS1, PARM1, TNFRSF17 
##     AKAP6, DERL3, SYN3, TP63, PPP1R9A, CD38, SSPN, MYO1D, PCDH9, IGKC 
##     THRB, AL139020.1, KLHL14, IRF4, RIMBP2, TSHR, CARMIL1, IGLC1, RAPGEF5, LINC01320 
## Negative:  AC007368.1, GPM6A, SOX5, GALNTL6, RGS7, AL592429.2, AL355076.2, AC083837.1, GEN1, ANK2 
##     SYT1, FCRL5, AC092546.1, MAST4, RHEX, CSGALNACT1, SOX5-AS1, TENM4, PRKD1, STEAP1B 
##     ADAMTS6, ZNF804A, AL450352.1, NTNG1, PLPP3, MIR3681HG, AL078459.1, DLGAP1, KCNJ3, AC093879.1 
## PC_ 5 
## Positive:  FCRL5, SOX5, PCDH9, MACROD2, AL079338.1, SYT1, TTC28, LINC01013, SOX5-AS1, MPP6 
##     MIR3681HG, NETO1, AL139020.1, GEN1, AC008691.1, IGHG3, LINC01374, ANK2, ROR1, RAPGEF5 
##     AC108879.1, SYN3, SKAP1, AC092546.1, PLPP3, ADAM23, PPP1R9A, MOXD1, PARM1, NRCAM 
## Negative:  SLC38A11, RGS7, STEAP1B, AC007368.1, GPM6A, JCHAIN, GALNTL6, PTPRG, RHEX, IGHA1 
##     IGHA2, AL355076.2, KLHL14, MAST4, EML6, STK33, TENM4, DERL3, MYO1D, NTNG1 
##     TNFRSF17, AL078459.1, IFNG-AS1, TXNDC5, ADAMTS6, AC083837.1, PRKD1, BCL2, ZNF804A, ARHGAP20
## Warning: Number of dimensions changing from 20 to 50
DimHeatmap(bcells, dims = 1:2, cells = 500, balanced = TRUE)

DimHeatmap(bcells, dims = 3:4, cells = 500, balanced = TRUE)

ElbowPlot(bcells)

#comb <- JackStraw(comb, num.replicate = 100)
bcells <- FindNeighbors(bcells, dims = 1:4)
## Computing nearest neighbor graph
## Computing SNN
bcells <- FindClusters(bcells, algorithm = 3, resolution = 0.3, verbose = FALSE)

bcells <- RunUMAP(bcells, dims = 1:4)
## 00:12:59 UMAP embedding parameters a = 0.9922 b = 1.112
## 00:12:59 Read 5314 rows and found 4 numeric columns
## 00:12:59 Using Annoy for neighbor search, n_neighbors = 30
## 00:12:59 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 00:13:00 Writing NN index file to temp file /tmp/RtmpjywqML/file27ed8f4f9d08ab
## 00:13:00 Searching Annoy index using 1 thread, search_k = 3000
## 00:13:01 Annoy recall = 100%
## 00:13:02 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 00:13:03 Initializing from normalized Laplacian + noise (using RSpectra)
## 00:13:04 Commencing optimization for 500 epochs, with 183396 positive edges
## 00:13:07 Optimization finished
DimPlot(bcells, reduction = "umap", label=TRUE,)

DimPlot(bcells, group.by="monaco_fine_annotation" , reduction = "umap", label=TRUE,)

bcells.markers <- FindAllMarkers(bcells, only.pos = TRUE)
## Calculating cluster 0
## Calculating cluster 1
## Calculating cluster 2
## Calculating cluster 3
## Calculating cluster 4
## Calculating cluster 5
## Calculating cluster 6
bcells.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1)
## # A tibble: 4,810 × 7
## # Groups:   cluster [7]
##       p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene     
##       <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>    
##  1 7.63e-20       1.06 0.441 0.381  2.27e-15 0       SNX9     
##  2 5.15e-11       1.27 0.208 0.148  1.53e- 6 0       YBX3     
##  3 9.16e-11       1.01 0.321 0.281  2.72e- 6 0       TMEM123  
##  4 4.37e-10       1.86 0.083 0.041  1.30e- 5 0       ARRDC4   
##  5 7.54e- 9       1.02 0.203 0.15   2.24e- 4 0       FOS      
##  6 1.25e- 8       1.05 0.145 0.094  3.72e- 4 0       USP9Y    
##  7 5.83e- 7       1.03 0.114 0.074  1.73e- 2 0       DDX3Y    
##  8 1.75e- 6       1.07 0.077 0.044  5.19e- 2 0       LINC00278
##  9 2.28e- 5       1.52 0.036 0.017  6.77e- 1 0       TTTY10   
## 10 4.53e- 5       1.08 0.166 0.134  1   e+ 0 0       RPGR     
## # ℹ 4,800 more rows
bcells.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1) %>%
    slice_head(n = 10) %>%
    ungroup() -> top10

DoHeatmap(bcells, features = top10$gene) + NoLegend()
## Warning in DoHeatmap(bcells, features = top10$gene): The following features
## were omitted as they were not found in the scale.data slot for the RNA assay:
## PRDM1, SDF2L1, SLAMF7, VIM, HIST1H4C, SAMHD1, MT-ND2, MT-ND1, PDE4D, NEK6,
## C12orf74, HIPK2, PTPRJ, CD247, PRKCH, RPGR, TTTY10, LINC00278, DDX3Y, USP9Y,
## FOS, ARRDC4, TMEM123, YBX3, SNX9

More featureplots.

# b cell markers
FeaturePlot(bcells, features = c("CD19","CD27","CD38","CD24"))

FeaturePlot(bcells, features = c("CR2", "CD34", "MME", "MS4A1"))

FeaturePlot(bcells, features = c("MZB1", "CXCR3", "FCRL5", "TBX21"))

FeaturePlot(bcells, features = c("CCR6", "ITGAX","MX1","BST2"))

# IGH genes
genes <- rownames(bcells)[grep("^IGH",rownames(bcells))]
genes
##  [1] "IGHEP2"     "IGHMBP2"    "IGHA2"      "IGHE"       "IGHG4"     
##  [6] "IGHG2"      "IGHGP"      "IGHA1"      "IGHEP1"     "IGHG1"     
## [11] "IGHG3"      "IGHD"       "IGHM"       "IGHJ6"      "IGHJ3P"    
## [16] "IGHJ5"      "IGHJ4"      "IGHV6-1"    "IGHV1-2"    "IGHV1-3"   
## [21] "IGHV4-4"    "IGHV7-4-1"  "IGHV2-5"    "IGHV3-7"    "IGHV3-64D" 
## [26] "IGHV5-10-1" "IGHV3-11"   "IGHV3-13"   "IGHV3-15"   "IGHV1-18"  
## [31] "IGHV3-20"   "IGHV3-21"   "IGHV3-23"   "IGHV1-24"   "IGHV2-26"  
## [36] "IGHV4-28"   "IGHV3-30"   "IGHV4-31"   "IGHV3-29"   "IGHV3-33"  
## [41] "IGHV4-34"   "IGHV4-39"   "IGHV3-43"   "IGHV1-46"   "IGHV3-48"  
## [46] "IGHV3-49"   "IGHV5-51"   "IGHV3-53"   "IGHV1-58"   "IGHV4-59"  
## [51] "IGHV3-64"   "IGHV3-66"   "IGHV1-69"   "IGHV2-70D"  "IGHV1-69-2"
## [56] "IGHV1-69D"  "IGHV2-70"   "IGHV3-72"   "IGHV3-73"   "IGHV3-74"  
## [61] "IGHV5-78"
lapply(seq(1,61,4), function(i) {
  j=i+3
  mygenes <- genes[i:j]
  FeaturePlot(bcells, features = mygenes )
})
## Warning: All cells have the same value (0) of "IGHEP2"
## Warning: All cells have the same value (0) of "IGHV3-29"
## Warning: All cells have the same value (0) of "IGHV1-58"
## Warning: All cells have the same value (0) of "IGHV3-66"
## Warning: All cells have the same value (0) of "IGHV1-69-2"
## Warning: The following requested variables were not found: NA
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

## 
## [[16]]

Analyse differences between samples in B cells

bcells@meta.data$covid <- grepl("M",bcells@meta.data$library)
Idents(bcells) <- "covid"

covid.de <- FindMarkers(bcells, ident.1 = "TRUE", ident.2 = "FALSE", verbose = TRUE)
head(covid.de,20)
##                    p_val avg_log2FC pct.1 pct.2     p_val_adj
## XIST       9.935340e-178 -1.4634068 0.287 0.700 2.952286e-173
## MALAT1     5.327151e-154 -0.4579590 0.986 0.999 1.582963e-149
## ARL17B     1.454177e-101 -2.8307874 0.037 0.222  4.321086e-97
## UTY         1.369802e-95 11.6504731 0.233 0.000  4.070368e-91
## DUSP1       4.019999e-95  2.4150213 0.403 0.126  1.194543e-90
## CXCR4       9.452891e-73  1.2802545 0.606 0.388  2.808926e-68
## PELI1       4.672331e-71  1.6305845 0.573 0.376  1.388383e-66
## RIPOR2      1.105943e-69 -0.5539199 0.771 0.920  3.286309e-65
## TSIX        4.940117e-66 -1.4471311 0.137 0.340  1.467956e-61
## HLA-DQA2    1.147568e-64  2.2683470 0.304 0.089  3.409998e-60
## PLCG2       7.114690e-64 -0.5670829 0.784 0.913  2.114130e-59
## GPM6A       4.241425e-62 -1.8532312 0.110 0.291  1.260339e-57
## TSC22D3     4.574745e-59  1.1956814 0.633 0.497  1.359386e-54
## RGS7        3.119841e-58 -2.0267471 0.069 0.221  9.270606e-54
## USP9Y       4.793243e-58  8.2360536 0.151 0.001  1.424312e-53
## AC007368.1  9.109379e-56 -2.6039285 0.031 0.148  2.706852e-51
## CSGALNACT1  5.953374e-53 -1.0873090 0.315 0.514  1.769045e-48
## GNLY        2.278711e-52  2.5926535 0.267 0.086  6.771189e-48
## ANK3        5.239092e-52 -1.0182389 0.230 0.430  1.556796e-47
## LINC02550   6.525279e-52 -2.4642119 0.024 0.127  1.938987e-47

Remove chrX and chrY genes

Shell command:

zcat gencode.v44.annotation.gtf.gz | grep chrX | cut -f9 | sed ‘s/gene_name //’ | grep @ | cut -d ‘"’ -f2 | uniq | sort -u > chrX_genes.txt

zcat gencode.v44.annotation.gtf.gz | grep chrY | cut -f9 | sed ‘s/gene_name //’ | grep @ | cut -d ‘"’ -f2 | uniq | sort -u > chrY_genes.txt

chrx <- readLines("ref/chrX_genes.txt")
chry <- readLines("ref/chrY_genes.txt")

covid.de <- covid.de[which(! rownames(covid.de) %in% chrx),]
covid.de <- covid.de[which(! rownames(covid.de) %in% chry),]

head(subset(covid.de,avg_log2FC>0),20) %>%
  kbl(caption="Upregulated in B cells after covid") %>%
  kable_paper("hover", full_width = F)
Upregulated in B cells after covid
p_val avg_log2FC pct.1 pct.2 p_val_adj
DUSP1 0 2.4150213 0.403 0.126 0
CXCR4 0 1.2802545 0.606 0.388 0
PELI1 0 1.6305845 0.573 0.376 0
HLA-DQA2 0 2.2683470 0.304 0.089 0
GNLY 0 2.5926535 0.267 0.086 0
EZR 0 0.9578868 0.709 0.631 0
YPEL5 0 1.7193243 0.374 0.198 0
KLF6 0 1.1403361 0.610 0.494 0
FOS 0 2.6020863 0.206 0.059 0
KLF2 0 1.2160893 0.484 0.329 0
CD69 0 1.0360232 0.519 0.381 0
EIF2AK3 0 1.0161844 0.548 0.416 0
NFKBIA 0 2.0214365 0.245 0.108 0
YBX3 0 2.0804468 0.200 0.072 0
JUNB 0 1.3445145 0.353 0.212 0
AC105402.3 0 0.7817414 0.475 0.309 0
JUN 0 1.2190982 0.417 0.284 0
HLA-C 0 0.6188876 0.677 0.611 0
RPS2 0 0.4887913 0.807 0.773 0
LYN 0 0.6114082 0.812 0.792 0 
head(subset(covid.de,avg_log2FC<0),20) %>%
  kbl(caption="Downregulated in B cells after covid") %>%
  kable_paper("hover", full_width = F)
Downregulated in B cells after covid
p_val avg_log2FC pct.1 pct.2 p_val_adj
MALAT1 0 -0.4579590 0.986 0.999 0
ARL17B 0 -2.8307874 0.037 0.222 0
RIPOR2 0 -0.5539199 0.771 0.920 0
PLCG2 0 -0.5670829 0.784 0.913 0
GPM6A 0 -1.8532312 0.110 0.291 0
RGS7 0 -2.0267471 0.069 0.221 0
AC007368.1 0 -2.6039285 0.031 0.148 0
CSGALNACT1 0 -1.0873090 0.315 0.514 0
ANK3 0 -1.0182389 0.230 0.430 0
LINC02550 0 -2.4642119 0.024 0.127 0
MBNL1 0 -0.4558884 0.835 0.924 0
MT-CO1 0 -0.2156893 0.989 0.992 0
GALNTL6 0 -2.1541408 0.035 0.134 0
DOCK10 0 -0.5854147 0.506 0.687 0
RNF17 0 -6.2601348 0.001 0.047 0
AC004687.1 0 -1.0713013 0.159 0.313 0
CPED1 0 -0.9486663 0.205 0.367 0
AL355076.2 0 -1.2214266 0.114 0.244 0
RALGPS2 0 -0.4363306 0.835 0.915 0
PPBP 0 -0.7972448 0.091 0.214 0 
VlnPlot(bcells, features <- c('DUSP1','CXCR4','PELI1','HLA-DQA2'),
  idents = c("FALSE", "TRUE"), group.by = "covid", ncol = 2)

VlnPlot(bcells, features <- c('ARL17B','GPM6A','RGS7','AC007368.1'),
  idents = c("FALSE", "TRUE"), group.by = "covid", ncol = 2)

Volcano

par(mar=c(5.1, 4.1, 4.1, 2.1))

HEADER="Effect of COVID infection on B cells"
plot(covid.de$avg_log2FC, -log10(covid.de$p_val),xlab="log2FC",ylab="-log10(p val)",pch=19,cex=0.8,main=HEADER)
sig <- subset(covid.de, p_val_adj<0.05)
points(sig$avg_log2FC, -log10(sig$p_val),pch=19,cex=0.8,col="red")
up=nrow(subset(sig,avg_log2FC>0))
dn=nrow(subset(sig,avg_log2FC<0))
nsig=nrow(sig)
ntot=nrow(covid.de)
SUBHEADER=paste("Total=",ntot,"genes;",nsig,"@5% FDR;",up,"up;",dn,"down")
mtext(SUBHEADER)

Heatmap of top genes.

top <- rownames(head(covid.de,15))

mx <- bcells[["RNA"]]["data"]
dim(mx)
## [1] 29715  5314
mx <- mx[which(rownames(mx) %in% top),]
dim(mx)
## [1]   15 5314
colfunc <- colorRampPalette(c("blue", "white", "red"))
colsidecols <- gsub("3","orange",gsub("2","gray", as.character(as.numeric(bcells@meta.data$covid)+2)))

heatmap.2(as.matrix(mx),col=colfunc(25), scale="row",trace="none",cexRow=0.8,
  cexCol=0.01,dendrogram="none",ColSideColors=colsidecols)

Pseudobulk

pseudo_bcells <- AggregateExpression(bcells, assays = "RNA", return.seurat = T, group.by = "library")
## Centering and scaling data matrix
Idents(pseudo_bcells) <- c("FALSE","FALSE","TRUE","TRUE","TRUE","TRUE","TRUE")

pseudo_bcells@meta.data$covid <- c("FALSE","FALSE","TRUE","TRUE","TRUE","TRUE","TRUE")

mx2 <- pseudo_bcells[["RNA"]]["counts"]
mx2 <- as.matrix(mx2)
dim(mx2)
## [1] 29715     7
mx2f <- mx2[which(rowMeans(mx2)>=10),]
dim(mx2f)
## [1] 11338     7
dds <- DESeqDataSetFromMatrix(countData = mx2f , colData = pseudo_bcells@meta.data , design = ~ covid )
## converting counts to integer mode
## Warning in DESeqDataSet(se, design = design, ignoreRank): some variables in
## design formula are characters, converting to factors
res <- DESeq(dds)
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
z<- results(res)
vsd <- vst(dds, blind=FALSE)
zz<-cbind(as.data.frame(z),assay(vsd))
dge<-as.data.frame(zz[order(zz$pvalue),])

head(subset(dge,log2FoldChange>0),20) %>%
  kbl(caption="Upregulated in B cells after covid") %>%
  kable_paper("hover", full_width = F)
Upregulated in B cells after covid
baseMean log2FoldChange lfcSE stat pvalue padj flw3 flw5 MDC8 MMC7 X1Recov3M X2Recov12M X2Recov3M
DUSP1 395.55815 2.4764963 0.4346142 5.698148 0.0000000 0.0000646 6.582319 7.257682 8.965373 9.461725 9.737484 8.243784 8.599466
TXNDC5 146.47800 1.9423747 0.4065961 4.777161 0.0000018 0.0037958 6.491252 6.066685 8.118962 7.266855 7.280582 8.310718 7.540042
JUNB 341.96398 1.3261093 0.2885687 4.595472 0.0000043 0.0057615 7.407362 7.814588 8.543616 9.198303 8.916368 8.430589 8.785369
ELL2 212.44250 1.8718481 0.4221863 4.433701 0.0000093 0.0109868 6.549910 6.845662 8.814079 8.328032 8.473961 7.761967 7.328673
HERPUD1 130.47015 1.0173990 0.2345095 4.338412 0.0000144 0.0152618 6.810331 6.615941 7.494263 7.354909 7.725897 7.588041 7.483969
DUSP5 22.14794 2.2378217 0.5379643 4.159796 0.0000319 0.0261564 4.906877 4.894999 5.771838 5.750654 6.285509 5.501499 5.804975
SEC61G 47.31752 1.4250258 0.3474793 4.101039 0.0000411 0.0313616 5.637449 5.649305 6.816338 6.476638 6.205853 6.416690 6.313656
YBX3 165.62580 1.9000438 0.4664356 4.073539 0.0000463 0.0329530 6.047595 6.758334 8.371819 8.466255 7.752942 7.298886 7.234394
RAB20 18.95957 3.3387606 0.8298575 4.023294 0.0000574 0.0382898 4.695479 4.450625 5.472012 5.896833 6.431670 5.128758 5.438513
IRS2 64.41087 1.9038183 0.4882451 3.899308 0.0000965 0.0510705 5.809811 5.366227 7.022173 7.043492 7.242693 6.060923 6.405831
SOCS3 12.03247 3.2040194 0.8253029 3.882235 0.0001035 0.0510705 4.642907 4.084235 5.334428 5.490062 5.708983 5.151229 5.477264
UCP2 149.97664 0.9329086 0.2408735 3.873024 0.0001075 0.0510705 6.762061 7.114117 7.736774 7.858991 7.641535 7.578439 7.586402
ISG15 33.85000 1.3780489 0.3636906 3.789070 0.0001512 0.0630285 5.327362 5.575581 6.322463 6.126733 5.932125 6.145795 6.151232
CXCR4 793.17877 1.4488624 0.3828790 3.784126 0.0001542 0.0630285 8.626848 8.576801 10.427065 10.483624 9.206415 9.310576 9.929458
RBM38 103.52866 1.6178946 0.4284766 3.775923 0.0001594 0.0630285 6.317717 5.933403 7.641388 7.423080 7.725897 6.523616 7.001582
MKI67 10.91401 4.6273760 1.2446548 3.717799 0.0002010 0.0766184 4.438889 4.084235 5.174966 5.312983 4.084235 5.801716 5.746448
PWP1 118.47655 1.0260538 0.2778518 3.692810 0.0002218 0.0816417 6.590301 6.648892 7.521359 7.550298 7.582379 7.023085 7.346795
NFKBIA 252.42803 2.1665559 0.5885538 3.681152 0.0002322 0.0826184 6.574289 6.788063 8.002564 8.969015 9.472207 7.310576 7.328673
NLRP3 28.72989 3.4498879 0.9427388 3.659431 0.0002528 0.0870445 4.788321 4.601014 5.231107 5.553400 7.353433 5.469770 5.438513
LINC01505 64.13954 1.0828392 0.2970429 3.645397 0.0002670 0.0890628 5.947510 6.066685 6.831054 6.575422 7.036207 6.635355 6.573041 
head(subset(dge,log2FoldChange<0),20) %>%
  kbl(caption="Downregulated in B cells after covid") %>%
  kable_paper("hover", full_width = F)
Downregulated in B cells after covid
baseMean log2FoldChange lfcSE stat pvalue padj flw3 flw5 MDC8 MMC7 X1Recov3M X2Recov12M X2Recov3M
SOX5 705.82052 -1.2561359 0.1791376 -7.012130 0.0000000 0.0000000 10.236080 10.305568 9.057044 9.051942 9.309653 9.038877 8.905264
LINC02550 46.74057 -2.0523302 0.3915757 -5.241210 0.0000002 0.0005677 7.330121 6.773279 5.738364 5.670222 5.261825 5.605844 6.129420
AC007368.1 86.96251 -2.7089394 0.5250238 -5.159650 0.0000002 0.0006603 7.597419 8.307365 5.897105 5.189549 5.579256 6.481863 6.214539
MSR1 15.23551 -2.0517712 0.4397875 -4.665370 0.0000031 0.0048193 5.866769 6.066685 5.115178 5.047080 5.054271 5.057700 4.997913
GALNTL6 125.45101 -2.1652405 0.4646314 -4.660125 0.0000032 0.0048193 7.954641 8.543098 6.384925 5.776287 7.280582 6.732747 6.234973
RGS7 151.60592 -1.6141461 0.3737955 -4.318260 0.0000157 0.0152618 8.319308 8.327319 6.902372 6.321537 7.582379 6.732747 7.194840
ARL17B 61.84379 -2.6122834 0.6264255 -4.170142 0.0000304 0.0261564 8.212793 6.254573 5.231107 6.182824 5.261825 5.739998 5.939965
AL450352.1 10.73580 -2.8385730 0.7179583 -3.953674 0.0000770 0.0483266 5.973292 5.613043 4.084235 4.410478 5.054271 4.923894 4.932090
MRVI1-AS1 37.74215 -1.2004468 0.3063022 -3.919158 0.0000889 0.0510705 6.456547 6.665073 5.866981 5.850046 5.708983 5.727258 5.887747
ARHGAP20 16.07999 -1.9132769 0.4895500 -3.908235 0.0000930 0.0510705 5.734585 6.187279 5.174966 5.047080 5.054271 4.923894 5.312819
EGLN3 23.03402 -1.8734757 0.4844369 -3.867327 0.0001100 0.0510705 6.447729 6.066685 5.703910 5.047080 5.432463 5.294949 4.997913
ADAMTS6 121.23221 -0.9915976 0.2588111 -3.831357 0.0001274 0.0566836 7.782105 7.785947 7.245664 6.963762 6.990982 6.999063 6.634525
LINC00996 21.75474 -1.5986022 0.4500602 -3.551974 0.0003824 0.1133780 6.140333 6.254573 5.115178 5.312983 4.776488 5.403337 5.716062
AC092279.1 17.83620 -1.5181387 0.4472592 -3.394315 0.0006880 0.1360084 5.960468 6.015059 5.174966 4.993816 5.054271 5.332381 5.514670
ZNF804A 154.94694 -0.9733200 0.2941799 -3.308588 0.0009377 0.1756087 7.730227 8.361585 7.319795 7.033771 7.488750 7.422440 6.990019
SAXO2 9.91328 -2.6592247 0.8089926 -3.287082 0.0010123 0.1801078 5.070553 6.210100 4.455508 4.936773 4.084235 4.923894 4.779756
AL355076.2 158.20850 -1.1058951 0.3386378 -3.265716 0.0010919 0.1855372 8.094995 8.261426 6.888399 6.753726 7.203756 7.656539 7.467532
SUPT3H 303.12176 -0.6456667 0.1977592 -3.264914 0.0010950 0.1855372 8.823673 8.756612 8.027868 8.339749 8.085953 8.156539 8.346245
CSGALNACT1 678.83374 -0.9260578 0.2877283 -3.218514 0.0012886 0.1977923 10.145142 9.948853 8.860459 8.696148 9.095203 9.424037 9.554550
AL079338.1 42.41566 -1.1769347 0.3677302 -3.200539 0.0013717 0.2009894 6.629524 6.743225 6.364430 5.776287 5.432463 5.825538 5.914133 
VlnPlot(bcells, features <- c('DUSP1','TXNDC5','JUNB','YBX3'),
  idents = c("FALSE", "TRUE"), group.by = "covid", ncol = 2)

VlnPlot(bcells, features <- c('SOX5','LINC02550','MSR1','GALNTL6'),
  idents = c("FALSE", "TRUE"), group.by = "covid", ncol = 2)

More heat

top <- rownames(head(dge,20))

rpm <- apply(mx2f,2,function(x) {x/sum(x)*1000000} )

hm <- rpm[rownames(rpm) %in% top,]

colsidecols <- gsub("TRUE","orange",gsub("FALSE","gray",pseudo_bcells@meta.data$covid))

heatmap.2(hm,col=colfunc(25), scale="row",trace="none",cexRow=1.2,
  cexCol=1.2,dendrogram="none",ColSideColors=colsidecols,mar=c(8,10),
  main="B cell genes after COVID")

table(bcells@meta.data$seurat_clusters)
## 
##    0    1    2    3    4    5    6 
## 1343 1039 1022  701  655  447  107
table(bcells@meta.data$monaco_fine_annotation)
## 
##           Exhausted B cells               Naive B cells 
##                         278                        3296 
## Non-switched memory B cells                Plasmablasts 
##                         912                          79 
##     Switched memory B cells 
##                         749
table(paste(bcells@meta.data$seurat_clusters,bcells@meta.data$monaco_fine_annotation))
## 
##           0 Exhausted B cells               0 Naive B cells 
##                            23                          1213 
## 0 Non-switched memory B cells                0 Plasmablasts 
##                            79                             1 
##     0 Switched memory B cells           1 Exhausted B cells 
##                            27                             3 
##               1 Naive B cells 1 Non-switched memory B cells 
##                           991                            44 
##     1 Switched memory B cells           2 Exhausted B cells 
##                             1                            36 
##               2 Naive B cells 2 Non-switched memory B cells 
##                           725                           135 
##     2 Switched memory B cells           3 Exhausted B cells 
##                           126                            42 
##               3 Naive B cells 3 Non-switched memory B cells 
##                           110                           496 
##     3 Switched memory B cells           4 Exhausted B cells 
##                            53                           105 
##               4 Naive B cells 4 Non-switched memory B cells 
##                            35                            94 
##                4 Plasmablasts     4 Switched memory B cells 
##                             8                           413 
##           5 Exhausted B cells               5 Naive B cells 
##                            53                           221 
## 5 Non-switched memory B cells                5 Plasmablasts 
##                            63                             5 
##     5 Switched memory B cells           6 Exhausted B cells 
##                           105                            16 
##               6 Naive B cells 6 Non-switched memory B cells 
##                             1                             1 
##                6 Plasmablasts     6 Switched memory B cells 
##                            65                            24
bcellseuratclusters <- lapply(unique(bcells@meta.data$seurat_clusters),function(x) {
  rownames(bcells@meta.data[which(bcells@meta.data$seurat_clusters == x),])
})

names(bcellseuratclusters) <- paste("S",unique(bcells@meta.data$seurat_clusters),sep="")

bcellsubtypes <- lapply(unique(bcells@meta.data$monaco_fine_annotation),function(x) {
  rownames(bcells@meta.data[which(bcells@meta.data$monaco_fine_annotation == x),])
})

names(bcellsubtypes) <- unique(bcells@meta.data$monaco_fine_annotation)

v1 <- c(bcellseuratclusters,bcellsubtypes)

plot(euler(v1),quantities = list(cex = 1.0), labels = list(cex = 1.5))

Now just zoom in on marker genes in FLW samples

flw_metainf <- meta_inf[grep("flw",meta_inf$cell_id),]
flw_metainf <- flw_metainf[which(flw_metainf$monaco_broad_annotation == "B cells"),]
flw <- comb[,which(colnames(comb) %in% rownames(flw_metainf))]


# remove non bcells
flw_metainf1 <- flw_metainf[grep("B cells",flw_metainf$monaco_fine_pruned_labels),]
flw_metainf2 <- flw_metainf[grep("Plasmablasts",flw_metainf$monaco_fine_pruned_labels),]
flw_metainf <- rbind(flw_metainf1,flw_metainf2)

flw <- flw[,which(colnames(flw) %in% rownames(flw_metainf))]

flw <- FindVariableFeatures(flw, selection.method = "vst", nfeatures = 2000)
## Finding variable features for layer counts
flw <- RunPCA(flw, features = VariableFeatures(object = bcells))
## Warning: The following 139 features requested have zero variance; running
## reduction without them: IGHG3, MZB1, TP63, INPP4B, PDE3B, DERL3, ZNF804A, GLDC,
## PITPNC1, PTCHD1-AS, TRABD2A, RASGRF2, DMXL2, DNAH8, SSR4, SGCD, RRBP1, KCNIP4,
## DLG2, MAGI2, GLIS3, NEGR1, AC243829.2, DTHD1, HESX1, ATXN1, IPCEF1, MERTK,
## LMAN1, BX322234.1, NCAM1, FAAH2, ADGRB3, TMEM232, AL353660.1, DOCK9,
## AC093916.1, TNFAIP3, VMP1, DAPK1, AL078602.1, IGLV3-1, TRBC1, AC068286.2,
## NR3C2, SLFN5, PTPRM, AC007100.1, SPANXA2-OT1, IGF2BP2, SPARCL1, ANXA1,
## IGHV3-73, CLEC19A, AL592429.2, PRICKLE2, AL589740.1, AC078881.1, PRMT2, GPC5,
## RAPGEF4-AS1, TBC1D4, ASPH, CCDC149, GZMK, CKAP2L, STPG2, LINC01983, PCAT4,
## MARVELD3, AC007938.3, LINC01146, RASGRP1, SUGCT, BAALC-AS1, AC089985.1, MN1,
## SLC44A1, CAV3, ADARB2, AGTPBP1, NMNAT3, C15orf62, CYTOR, VIM, NEXMIF,
## AL049828.1, LINC01937, LINC01182, HNRNPLL, AC097654.1, ANKRD34B, EMP1, CNTLN,
## SND1, NETO2, AC090844.2, AC107067.2, LINC00423, EHD2, LINC00159, PPEF1,
## MIR193BHG, GIMAP7, MAP4K3-DT, CSF1R, F11-AS1, CRMP1, NIPAL3, PHEX, IGF1,
## AC009120.4, PRH1, HIST1H2AH, AC004895.1, EPHX3, PNCK, AC114956.1, LINC01297,
## CORT, ZDHHC11B, EFNB2, ADGRA3, EEF1AKMT1, GBP2, AC008474.1, GINS1, VEGFA,
## PRKCA-AS1, TLR5, AC090945.1, NPAS2, AC021231.1, GNLY, AL139317.5, RNF213,
## ZFYVE9, ABCC4, MAP3K13
## Warning in PrepDR5(object = object, features = features, layer = layer, : The
## following features were not available: VCAN, PLXDC2, MAN1A1, NAMPT, FNDC3B,
## HSPA5, LRMDA, XBP1, PLCB1, AOAH, PRKCH, SLC8A1, LINC01505, CD247, PRDM1,
## SDF2L1, NUSAP1, RBM47, AAK1, NEAT1, DPYD, CENPE, DCC, ACSL1, GPC6, TNIK, MANF,
## NCALD, PCSK5, FKBP11, STAT4, TOX, RORA, PPP2R2B, SLCO3A1, MYDGF, SEC24D, GNAQ,
## SNX32, CDC14A, ITGB1, SHCBP1, CEP128, CSF3R, IQGAP2, ARHGAP26, CALR, VCAN-AS1,
## GTSE1, PDIA4, CREB5, DOCK4, SLC11A1, A2M, STXBP6, SLAMF7, PDE4D, TFEC, RAB39A,
## AC008574.1, GAB1, AC079336.5, SSR3, MCTP1, MYBL1, KIF4A, AGAP1, AL442163.1,
## AL138716.1, HIST1H4C, WWC1, PLAAT2, SMC4, SPATS2, LINC02384, MAML3, EIF4E3,
## AC044781.1, SEC11C, LINC02296, TPRG1, MAGI1, ERN1, MCM10, NLGN4Y, CADPS2, FYN,
## LCP2, TOX2, SAMHD1, LINC02456, GPRIN3, PDIA6, SORL1, LINC02363, LINC01821,
## GLT1D1, CLEC12A, PPP1R26-AS1, SNAP25-AS1, ESRRB, WDFY3, IRAK3, TRERF1, IL6ST,
## SLC22A15, PTPRE, SRGN, AC002072.1, LRRC8C, ARLNC1, CHODL-AS1, ME1, EPB41L4A,
## AC091938.1, COL12A1, AL807761.4, AC037198.1, ZNF239, LRRC4C, PWRN1, CCR9,
## NDFIP1, ROR2, IGF2BP3, ANO5, DENND1B, TYMS, BTBD11, RCAN2, WNK3, RGS13,
## AL356108.1, SEMA3B, AC026353.1, AL442636.1, CFAP54, NFIA, ESR1, AC073332.1,
## RIN2, FKBP5, AC006581.2, FBXO5, INPP4A, EXT1, MAFTRR, ELL2, RUNX1, CERS6, KNL1,
## MIR181A1HG, MSC-AS1, NIBAN1, SEMA3C, MYO1F, ADAMTSL4-AS1, AC108169.1,
## AADACL2-AS1, MAMDC2-AS1, SPART-AS1, SAMD3, SYNE2, ESCO2, AC011995.2, SPC25,
## Z99758.1, TCF7L1, LEPROTL1, IL20RB-AS1, LDLRAD3, CENPP, AL358777.1, DOCK5,
## DEPDC1B, MAGI3, HIVEP3, TAFA2, TRPA1, ANXA10, LINC02211, KLF12, SPN, LINC00334,
## ATF5, TGFBR3, PIK3R5, AP001977.1, KRT1, ADCY9, MPP7, KIFC1, PPARG, CRELD2,
## VWC2, TLL2, AC079352.1, AC009974.2, AC027702.1, AC108925.1, TNC, TRAV8-3,
## CDC45, SEPTIN4-AS1, AC109811.1, ART4, MOCOS, UPK1A-AS1, ZNF541, AC026979.4,
## AP001831.1, RBMS1, AC005355.1, AC024558.2, AC104964.3, TNFAIP2, FRY, ARHGAP10,
## AC009975.1, NCAPG, EEPD1, CCL5, SULT2B1, NUCB2, ARL17B, AC011416.3, NIM1K,
## PTPRJ, PROX1, AC092611.1, TRIP13, TCERG1L, AC084809.1, AC090376.1, AL646090.2,
## AC012459.1, NHS, AC073050.1, CCDC113, SVIL, ADGRE2, DIO1, LINC01847,
## AL451166.1, TNFRSF10C, AC103987.2, CDH2, BMX, LHFPL1, LINC00907, OPRM1, FANCC,
## SPOCK1, TBXAS1, AL360093.1, AC019117.2, TIAM1, CD226, AC073114.1, E2F1, CNGA1,
## AC007262.2, ARHGAP6, MVB12B, ZEB2, CYP27A1, LINC01163, PF4V1, TGFB2, LINC01169,
## CDCA5, PDGFD, FRMD4B, SUB1, RFX2, AC025263.1, PARP8, TXNDC11, PERP, ZCCHC24,
## NLRP3, KCNN3, HIST1H2AJ, AC008170.1, C1orf189, FCN1, DENND2A, KLRG1, SYTL3,
## TLR2, CDHR3, STK3, POU6F2, GAB2, PLA2G4A, COL18A1, HIST1H3C, AC020978.2,
## LINC02080, AC011509.2, GATA1, C8orf88, C14orf178, NRG3, TSPAN9, TAF1A-AS1,
## MUC13, PCDHGA12, OR1L8, LINC02625, AC079601.1, AL121772.3, CCN5, AC024382.1,
## CABCOCO1, SFMBT2, TENM1, POF1B, SYTL2, BRIP1, SDSL, FER1L5, SULF1, AC018845.1,
## PYHIN1, RBL2, IGLV1-51, C8orf37-AS1, DPYD-AS1, FAM222A-AS1, AC073475.1,
## SLC26A7, TRAV17, NPIPA8, VSTM1, IL6R, CACNB4, GAS7, AL096854.1, ATP8B4, CD8A,
## MRC2, STK32B, JAKMIP2-AS1, AQP9, RNF175, CCDC88A, CFAP44-AS1, AL445928.2,
## CTBP2, UPP2, LINC02137, SMYD3, LINC01605, AC135803.1, CLCA4, AL592402.1, FRZB,
## LINC01878, FAM124B, AC023503.1, ODAPH, LINC01511, AC010280.1, LINC02526,
## AC005682.2, VPS37D, AC105206.3, AL512634.1, AL450322.1, ZNF214, AC090099.1,
## AC092111.1, CLEC1A, AC078962.2, IL22, AL355516.1, C13orf42, AL358334.4,
## LINC00911, AC087639.2, MESP1, KCNJ12, AC015911.2, AC115100.1, SLC14A2-AS1,
## AC138470.2, HIF3A, LINC01620, KCNS1, AL162457.1, IGLV3-9, AL023574.1, GDPD2,
## SLC1A7, SMG7-AS1, MINDY4B, LINC01276, AL591468.1, GHET1, CARD17, AC121338.1,
## PAH, PCDH9-AS1, LINC00433, MCF2L-AS1, DSCAS, LINC01919, AP000265.1, AC245140.3,
## SELENOS, AL109840.2, MICAL2, SBF2, SOBP, LINC02008, KIF13A, FGD4, CADPS,
## HIST1H2BB, NMT2, KIAA0825, TANC2, TSHZ3, BUB1, ICOS, AL353151.2, LINC01229,
## BEX5, AL139294.1, CDS1, FOLR3, UBE2J1, NHSL1, NLRP12, SLC44A3, STOX2,
## AP000593.4, AC074051.5, LINC02860, FGD6, BUB1B, AC105180.2, LRGUK, AC010609.1,
## C1orf21, AC104459.1, AC079760.2, AP003110.1, CCDC7, H2AFX, IL4I1, AC023480.1,
## TRAIP, LYPLAL1-DT, JAKMIP2, TRPS1, MCC, AC092164.1, ASTL, VOPP1, TNR, SCOC-AS1,
## BMPR1A, AC026341.1, CLDN11, AL024498.1, LINC02676, IGHV5-10-1, DLGAP1-AS5,
## ZNF534, AC002472.1, APOBEC3H, ASMT, ST8SIA6, WDR64, OSCAR, LINC02328, TEX41,
## AC003044.1, NEDD4, TSPAN5, TENT5A, KIF23, ALPL, NCAPG2, C16orf71, FBXO15,
## WNT2B, HIPK2, SPCS2, PCNX1, SNTB1, PAM, TVP23C, ZBTB16, MYO3A, CEP70, FANCI,
## L2HGDH, KCNMB2-AS1, RGL1, DNM3, SIPA1L2, AC005833.2, ZYG11A, VIPR2, AP003071.2,
## RTN4RL1, PTPN4, CRADD, AL365295.1, GASK1A, ALDH1L2, AUXG01000058.1, PIP4P2,
## GLRA2, AC109927.1, C8orf31, AL355303.1, ID2, CEP85L, AC125618.1, NEIL3,
## AL022724.3, AC019257.1, BASP1-AS1, AL354733.2, CTNNA3, MARCO, DACT3, CFAP61,
## AC007610.1, PPP4R4, CMTM2, LPCAT2, POLE2, ESYT2, ETS1, MIR4435-2HG, AC068643.1,
## LINC01933, PPM1L, PLEKHA5, CYSLTR2, AC066613.2, LMAN2, UTY, SCLT1, TNFSF8,
## AL359644.1, ERC1, AC015923.1, ARHGAP22, AC097460.1, CMTM4, CLEC7A, CPVL, CLIP4,
## SPAG16, AL365255.1, RGS1, ANKRD7, ATP8A2, PRAM1, TUBA1B, LGALSL-DT, LINC02234,
## CHN1, ZRANB2-AS2, TMEM252-DT, KLRC2, LINC00278, TRAT1, AC108718.1, DIAPH2-AS1,
## CIT, SGMS2, ZBP1, SEMA4D, C6orf58, HACD1, SCN2A, MCF2, MAP7, MLEC, FRMD3,
## ASAP3, AL512288.1, AC108210.1, ZNF391, PFKFB1, HLF, CD3G, ME3, CCSER2, HLCS,
## AC020916.1, LINC00298, AL357054.5, AL357141.1, TSPAN11, AC007611.1, AC100863.1,
## KLHL34, LINC00571, ISOC2, SGSM1, CPQ, MYO5B, PLAUR, GAS2, AL591518.1, DACH2,
## EVA1C, C20orf194, FOS, PLS1, SSC5D, AC107068.2, AL358937.1, PAK1, HYOU1, TIMP3,
## KIF2C, PRLR, TRG-AS1, TLN2, C4orf36, CERS6-AS1, RLN2, AL354989.1, CPEB1, DHRS3,
## NAV1, SLC24A4, MTHFD1L, SYNE1, AL357873.1, COL16A1, ERICH3-AS1, AC095030.1,
## LINC02609, AL590666.4, SERPINC1, TEX35, NEK2, WNT3A, STON1-GTF2A1L, LINC01793,
## AC016747.2, IGKV5-2, KLHL41, AC108066.1, GRM7-AS3, ZNF662, CADM2-AS2, IGSF10,
## AC092944.3, TRPC3, AC104083.1, AC026726.1, LINC01942, AL365226.1, CFAP206,
## TRDN-AS1, ITPRID1, TRGV9, AC005076.1, AC067930.1, AL445526.1, FAM201A,
## AL353742.1, AKR1C4, LINC00707, AL354916.1, ASAH2, PLEKHS1, AL139123.1,
## AC104383.2, CBLIF, NECTIN1-AS1, GPD1, LINC02389, AC091214.1, AC089999.2,
## CCDC169, TNFSF11, DZIP1, AL358332.1, C15orf54, AC103740.1, CCDC33, AC068870.2,
## KIF7, AC022819.1, CERS3-AS1, AC090907.1, AC120498.4, AC012676.1, AC099518.5,
## AC015908.2, FBXW10, TTLL6, CACNG1, KIF19, AC091691.2, AC011446.1, AC008747.1,
## AC006213.1, AL354993.2, AP000355.1, APOL4, AL683807.2, DIPK2B, RAB40A,
## RNASEH2B-AS1, TMEM170B, HDAC4, ZEB1-AS1, LINC02798, OSTN-AS1, TTTY10, LRIG1,
## ARSG, ATP2B4, TMEM163, JAZF1-AS1, WWC2, AC122697.1, DNAJC1, MANEA, WDR17,
## AC079174.2, AC008870.4, LINC01692, NR6A1, CA3, DST, SLC2A13, CCDC144A, P4HB,
## SLC5A9, PTGFR, AL391811.1, LINC01816, SULT1C4, LINC01173, CNTN4-AS1, THOC7-AS1,
## SPTSSB, LINC02505, AC096759.1, AC116616.1, DNAH5, LINC02060, LINC01622, TFAP2A,
## POPDC3, LINC00326, AL139393.1, AC005100.1, PCOLCE-AS1, AC073314.1, LINC00681,
## CDH17, AL157829.1, NRARP, PITX3, TAC3, AVPR1A, AC010201.1, TRAV26-1, FLRT2,
## AL355836.2, IGHV4-28, AC090617.1, TMIGD1, AC022903.2, AP001029.1, MEP1B, CNN1,
## AC010463.3, DLL3, LINC00945, BRWD1-AS2, LINC01424, AP001468.1, TCN2, ARSF,
## TLR8-AS1, FMR1NB, AL133334.1, AC004817.5, LINC02133, KLRK1, DGKG, PSD3, CD2,
## ARNTL2, KIF9-AS1, AC004158.1, LARP1B, GATA3, RAB31, SSR1, NUGGC, MIR646HG,
## LINC01625, AC025569.1, RFX8, MYBL2, AL596218.1, PSTPIP2, UACA, ATL1, ALG14,
## MS4A4E, AL035427.1, WAKMAR2, ITPKB, AC139769.2, LMNB1, SYBU, ARRB1, CLIC5,
## FRRS1, AC119674.1, BCAT1, NMRK1, MBNL1-AS1, RAB3C, AC018695.2, LINC02246,
## RAB20, AC097515.1, VCL, CCDC30, AMPH, MELK, AC021086.1, HIST1H2AL, UGGT2, SAT1,
## DDAH1, FSD1L, MCM4, RPS6KA6, AC104170.1, STK38, GCNT4, ANO10, PCBP3, TRDC,
## DENND2C, TK1, LINC01353, SLC9A2, RNF180, AL158801.2, AC136431.1, ITGA2, TECTA,
## KANK2, SMPDL3A, CDC6, ZSCAN31, AC096536.3, AMBN, EGR3, CLEC6A, TMPRSS12,
## AC122685.1, FERMT2, AC027458.1, PLAC1, SLC8A1-AS1, PXDNL, HM13, NCOA1, HDX,
## TACR1, HIST1H1B, NT5DC2, CDCA8, BRCA1, FBXO32, ZNF438, LINC01118, LINC00299,
## SAMSN1, CD59, ZNF658, CEP112, AL162493.1, ARHGAP18, POT1-AS1, AL138828.1,
## SLC2A9, HPGD, KIAA0319, CCNH, CPEB4, AP001021.3, XRCC4, BMP8B, ZC2HC1B, ACVR2A,
## LINC02391, BNC2, KALRN, HIST1H2AG, GSTCD, CCL4, LINC00393, AC087482.1, SUSD5,
## AC003666.1, PIP4K2A, FBXW7, FOXN3, ETV4, AC104823.1, STPG2-AS1, KCTD16, LEPR,
## RAP1GA
## PC_ 1 
## Positive:  FYB1, ITK, BCL11B, PRKCA, IL32, CD36, THEMIS, PID1, TCF7L2, RTN1 
##     PRKCQ, JAML, SLFN12L, PHACTR2, DISC1, MS4A6A, ATP10A, CPPED1, SKAP1, TMTC2 
##     TC2N, PZP, UBASH3B, MGAT4A, TXK, CD6, SERINC5, PAG1, PLCL1, TRAC 
## Negative:  GPM6A, SOX5, AC007368.1, AL355076.2, IGHM, STEAP1B, RGS7, SLC38A11, GALNTL6, FCRL5 
##     RHEX, IGKC, AC083837.1, PCDH9, LINC02550, SSPN, CSGALNACT1, KLHL14, KCNQ5, GEN1 
##     MACROD2, AL078459.1, AL079338.1, IGLC2, FUT8, XIST, RIPOR2, SYT1, PARM1, ADAMTS6 
## PC_ 2 
## Positive:  SOX5, AL355076.2, SSPN, RHEX, GALNTL6, AC008691.1, AC007368.1, SYT1, GEN1, AC092546.1 
##     IGHG1, GPM6A, ANK2, SOX5-AS1, CSGALNACT1, FCRL5, EML6, IGHA2, MIR3681HG, TEX9 
##     IGHA1, PLPP3, AL450352.1, IGHGP, TENM4, JCHAIN, CRIP1, MAST4, OSTN, ITM2C 
## Negative:  SLC38A11, AKAP6, STEAP1B, AC108879.1, SYN3, LINC02550, PCDH9, IGHM, IFNG-AS1, KLHL14 
##     MACROD2, NETO1, PLD5, SATB1-AS1, CARMIL1, RGS7, CCR7, SKAP1, LINC01811, LINC01340 
##     FOXP1, PRKN, AL139020.1, AL079338.1, LEF1, PARM1, PCAT1, CCND3, IGF1R, MLLT3 
## PC_ 3 
## Positive:  PCDH9, MACROD2, SYN3, AL079338.1, AC108879.1, IGHM, NETO1, AL139020.1, FCRL5, SLC38A11 
##     LINC01013, AKAP6, IGKC, RAPGEF5, LINC01374, TTC28, STEAP1B, LINC02550, IFNG-AS1, AC105402.3 
##     LINC01340, PARM1, THRB, MOXD1, MPP6, IGLC1, KLHL14, AL133346.1, SYT1, TXNDC5 
## Negative:  AC007368.1, AL355076.2, LEF1, CAMK4, GALNTL6, CSGALNACT1, PRKCA, BCL11B, MAST4, ANK3 
##     NELL2, FYB1, SERINC5, TCF7, TC2N, IL7R, TXK, PLCL1, GPM6A, FHIT 
##     RHEX, TESPA1, TAFA1, ADAMTS6, THEMIS, ITK, TSHZ2, ITGA6, DOCK10, MLLT3 
## PC_ 4 
## Positive:  AC007368.1, RGS7, GPM6A, SLC38A11, GALNTL6, STEAP1B, MAST4, TENM4, STK33, KLHL14 
##     AC083837.1, ACSM3, PTPRG, FOXP1, NTNG1, ARHGAP20, DLGAP1, DCLK2, IGHM, KCNJ3 
##     AC105402.3, SEMA3D, PRKD1, LINC01572, TMEM182, MSR1, IGF1R, MYO3B, AC005699.1, NAALADL2-AS2 
## Negative:  SOX5, FCRL5, SSPN, PCDH9, AL079338.1, PARM1, SYT1, MACROD2, AC008691.1, TEX9 
##     LINC01013, AL589693.1, MIR3681HG, TTC28, MPP6, NETO1, IGHA2, SOX5-AS1, SYN3, IFNG-AS1 
##     NAV2, AC099560.1, AC108879.1, PPP1R9A, MOXD1, IGHG1, IGHGP, EPHA4, PLPP3, IGKC 
## PC_ 5 
## Positive:  SOX5, FCRL5, AL079338.1, AC007368.1, SYT1, ANK2, GALNTL6, PCDH9, SOX5-AS1, MACROD2 
##     NETO1, LINC01013, AC092546.1, MPP6, PLPP3, AC093879.1, GEN1, MIR3681HG, SMOC2, GPM6A 
##     TTC28, SKAP1, KCNJ3, FOXP1, SYN3, PCSK6, ROR1, FYB1, AL450352.1, BCL11B 
## Negative:  SSPN, AL355076.2, PARM1, IFNG-AS1, TEX9, IGHA2, JCHAIN, RASSF6, AC008691.1, NAV2 
##     AL589693.1, IGHG1, AL078459.1, IGHA1, EML6, LINC01811, SNED1, RHEX, KLHL14, MYO1D 
##     EPHA4, IGHGP, SLCO5A1, AC008014.1, EYA2, SAMD12, AKAP6, PPP1R9A, IGKC, LINC01320
## Warning: Number of dimensions changing from 20 to 50
DimHeatmap(flw, dims = 1:2, cells = 500, balanced = TRUE)

DimHeatmap(flw, dims = 3:4, cells = 500, balanced = TRUE)

ElbowPlot(flw)

#comb <- JackStraw(comb, num.replicate = 100)
flw <- FindNeighbors(flw, dims = 1:4)
## Computing nearest neighbor graph
## Computing SNN
flw <- FindClusters(flw, algorithm = 3, resolution = 0.3, verbose = FALSE)

flw <- RunUMAP(flw, dims = 1:4)
## 00:28:16 UMAP embedding parameters a = 0.9922 b = 1.112
## 00:28:16 Read 1558 rows and found 4 numeric columns
## 00:28:16 Using Annoy for neighbor search, n_neighbors = 30
## 00:28:16 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 00:28:16 Writing NN index file to temp file /tmp/RtmpjywqML/file27ed8f2890f821
## 00:28:16 Searching Annoy index using 1 thread, search_k = 3000
## 00:28:17 Annoy recall = 100%
## 00:28:18 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 00:28:19 Initializing from normalized Laplacian + noise (using RSpectra)
## 00:28:19 Commencing optimization for 500 epochs, with 53734 positive edges
## 00:28:21 Optimization finished
DimPlot(flw, reduction = "umap", label=TRUE,)

DimPlot(flw, group.by="monaco_fine_annotation" , reduction = "umap", label=TRUE,)

flw.markers <- FindAllMarkers(flw, only.pos = TRUE)
## Calculating cluster 0
## Calculating cluster 1
## Calculating cluster 2
## Calculating cluster 3
## Calculating cluster 4
flw.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1)
## # A tibble: 2,640 × 7
## # Groups:   cluster [5]
##        p_val avg_log2FC pct.1 pct.2 p_val_adj cluster gene    
##        <dbl>      <dbl> <dbl> <dbl>     <dbl> <fct>   <chr>   
##  1 4.99e-106       1.57 0.984 0.798 1.48e-101 0       BACH2   
##  2 1.12e- 74       1.88 0.839 0.549 3.32e- 70 0       COL19A1 
##  3 1.55e- 67       2.36 0.625 0.247 4.62e- 63 0       IL4R    
##  4 4.66e- 63       2.22 0.647 0.279 1.39e- 58 0       LIX1-AS1
##  5 6.95e- 63       1.36 0.861 0.685 2.07e- 58 0       SNX29   
##  6 7.04e- 54       1.27 0.865 0.774 2.09e- 49 0       FCRL1   
##  7 1.75e- 49       1.19 0.885 0.764 5.19e- 45 0       ADK     
##  8 5.28e- 47       1.17 0.885 0.802 1.57e- 42 0       MEF2C   
##  9 6.61e- 41       1.63 0.637 0.396 1.96e- 36 0       IGHD    
## 10 2.61e- 36       1.10 0.788 0.665 7.76e- 32 0       ITPR1   
## # ℹ 2,630 more rows
flw.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1) %>%
    slice_head(n = 10) %>%
    ungroup() -> top10

DoHeatmap(flw, features = top10$gene) + NoLegend()
## Warning in DoHeatmap(flw, features = top10$gene): The following features were
## omitted as they were not found in the scale.data slot for the RNA assay:
## AC021055.1, HYAL3, AC066580.1, LINC01143, AL161640.1, AC105275.2, LINC01768,
## LINC01739, AL137798.2, MT-ND2, PTPRJ, CD247, PRKCH, FAM49A, C12orf74, ITPR1,
## ADK, SNX29

More B cell heat maps.

# b cell markers
FeaturePlot(flw, features = c("CD19","CD27","CD38","CD24"))

FeaturePlot(flw, features = c("CR2", "CD34", "MME", "MS4A1"))

FeaturePlot(flw, features = c("MZB1", "CXCR3", "FCRL5", "TBX21"))

FeaturePlot(flw, features = c("CCR6", "ITGAX","MX1","BST2"))

# IGH genes
genes <- rownames(bcells)[grep("^IGH",rownames(bcells))]
genes
##  [1] "IGHEP2"     "IGHMBP2"    "IGHA2"      "IGHE"       "IGHG4"     
##  [6] "IGHG2"      "IGHGP"      "IGHA1"      "IGHEP1"     "IGHG1"     
## [11] "IGHG3"      "IGHD"       "IGHM"       "IGHJ6"      "IGHJ3P"    
## [16] "IGHJ5"      "IGHJ4"      "IGHV6-1"    "IGHV1-2"    "IGHV1-3"   
## [21] "IGHV4-4"    "IGHV7-4-1"  "IGHV2-5"    "IGHV3-7"    "IGHV3-64D" 
## [26] "IGHV5-10-1" "IGHV3-11"   "IGHV3-13"   "IGHV3-15"   "IGHV1-18"  
## [31] "IGHV3-20"   "IGHV3-21"   "IGHV3-23"   "IGHV1-24"   "IGHV2-26"  
## [36] "IGHV4-28"   "IGHV3-30"   "IGHV4-31"   "IGHV3-29"   "IGHV3-33"  
## [41] "IGHV4-34"   "IGHV4-39"   "IGHV3-43"   "IGHV1-46"   "IGHV3-48"  
## [46] "IGHV3-49"   "IGHV5-51"   "IGHV3-53"   "IGHV1-58"   "IGHV4-59"  
## [51] "IGHV3-64"   "IGHV3-66"   "IGHV1-69"   "IGHV2-70D"  "IGHV1-69-2"
## [56] "IGHV1-69D"  "IGHV2-70"   "IGHV3-72"   "IGHV3-73"   "IGHV3-74"  
## [61] "IGHV5-78"
lapply(seq(1,61,4), function(i) {
  j=i+3
  mygenes <- genes[i:j]
  FeaturePlot(flw, features = mygenes )
})
## Warning: All cells have the same value (0) of "IGHEP2"
## Warning: All cells have the same value (0) of "IGHV7-4-1"
## Warning: All cells have the same value (0) of "IGHV3-64D"
## Warning: All cells have the same value (0) of "IGHV5-10-1"
## Warning: All cells have the same value (0) of "IGHV1-24"
## Warning: All cells have the same value (0) of "IGHV4-28"
## Warning: All cells have the same value (0) of "IGHV3-29"
## Warning: All cells have the same value (0) of "IGHV3-49"
## Warning: All cells have the same value (0) of "IGHV1-58"
## Warning: All cells have the same value (0) of "IGHV3-64"
## Warning: All cells have the same value (0) of "IGHV3-66"
## Warning: All cells have the same value (0) of "IGHV1-69-2"
## Warning: The following requested variables were not found: NA
## Warning: All cells have the same value (0) of "IGHV5-78"
## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

## 
## [[5]]

## 
## [[6]]

## 
## [[7]]

## 
## [[8]]

## 
## [[9]]

## 
## [[10]]

## 
## [[11]]

## 
## [[12]]

## 
## [[13]]

## 
## [[14]]

## 
## [[15]]

## 
## [[16]]

Session information

sessionInfo()
## R version 4.4.0 (2024-04-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 22.04.4 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.10.0 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.10.0
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## time zone: Australia/Melbourne
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] eulerr_7.0.2                gplots_3.1.3.1             
##  [3] kableExtra_1.4.0            SingleR_2.6.0              
##  [5] celldex_1.14.0              harmony_1.2.0              
##  [7] Rcpp_1.0.12                 mitch_1.16.0               
##  [9] DESeq2_1.44.0               muscat_1.18.0              
## [11] beeswarm_0.4.0              stringi_1.8.4              
## [13] SingleCellExperiment_1.26.0 SummarizedExperiment_1.34.0
## [15] Biobase_2.64.0              GenomicRanges_1.56.0       
## [17] GenomeInfoDb_1.40.0         IRanges_2.38.0             
## [19] S4Vectors_0.42.0            BiocGenerics_0.50.0        
## [21] MatrixGenerics_1.16.0       matrixStats_1.3.0          
## [23] hdf5r_1.3.10                Seurat_5.0.3               
## [25] SeuratObject_5.0.2          sp_2.1-4                   
## [27] plyr_1.8.9                  ggplot2_3.5.1              
## [29] dplyr_1.1.4                
## 
## loaded via a namespace (and not attached):
##   [1] progress_1.2.3            goftest_1.2-3            
##   [3] Biostrings_2.72.0         HDF5Array_1.32.0         
##   [5] vctrs_0.6.5               spatstat.random_3.2-3    
##   [7] digest_0.6.35             png_0.1-8                
##   [9] corpcor_1.6.10            shape_1.4.6.1            
##  [11] gypsum_1.0.1              ggrepel_0.9.5            
##  [13] echarts4r_0.4.5           deldir_2.0-4             
##  [15] parallelly_1.37.1         MASS_7.3-60.2            
##  [17] reshape2_1.4.4            httpuv_1.6.15            
##  [19] foreach_1.5.2             withr_3.0.0              
##  [21] ggrastr_1.0.2             xfun_0.43                
##  [23] survival_3.6-4            memoise_2.0.1.9000       
##  [25] ggbeeswarm_0.7.2          systemfonts_1.0.6        
##  [27] zoo_1.8-12                GlobalOptions_0.1.2      
##  [29] gtools_3.9.5              pbapply_1.7-2            
##  [31] prettyunits_1.2.0         GGally_2.2.1             
##  [33] KEGGREST_1.44.0           promises_1.3.0           
##  [35] httr_1.4.7                globals_0.16.3           
##  [37] fitdistrplus_1.1-11       rhdf5filters_1.16.0      
##  [39] rhdf5_2.48.0              rstudioapi_0.16.0        
##  [41] UCSC.utils_1.0.0          miniUI_0.1.1.1           
##  [43] generics_0.1.3            curl_5.2.1               
##  [45] zlibbioc_1.50.0           ScaledMatrix_1.12.0      
##  [47] polylabelr_0.2.0          polyclip_1.10-6          
##  [49] GenomeInfoDbData_1.2.12   ExperimentHub_2.12.0     
##  [51] SparseArray_1.4.3         xtable_1.8-4             
##  [53] stringr_1.5.1             doParallel_1.0.17        
##  [55] evaluate_0.23             S4Arrays_1.4.0           
##  [57] BiocFileCache_2.12.0      hms_1.1.3                
##  [59] irlba_2.3.5.1             colorspace_2.1-0         
##  [61] filelock_1.0.3            ROCR_1.0-11              
##  [63] reticulate_1.36.1         spatstat.data_3.0-4      
##  [65] magrittr_2.0.3            lmtest_0.9-40            
##  [67] later_1.3.2               viridis_0.6.5            
##  [69] lattice_0.22-6            spatstat.geom_3.2-9      
##  [71] future.apply_1.11.2       scattermore_1.2          
##  [73] scuttle_1.14.0            cowplot_1.1.3            
##  [75] RcppAnnoy_0.0.22          pillar_1.9.0             
##  [77] nlme_3.1-164              iterators_1.0.14         
##  [79] caTools_1.18.2            compiler_4.4.0           
##  [81] beachmat_2.20.0           RSpectra_0.16-1          
##  [83] tensor_1.5                minqa_1.2.6              
##  [85] crayon_1.5.2              abind_1.4-5              
##  [87] scater_1.32.0             blme_1.0-5               
##  [89] locfit_1.5-9.9            bit_4.0.5                
##  [91] codetools_0.2-20          BiocSingular_1.20.0      
##  [93] bslib_0.7.0               alabaster.ranges_1.4.0   
##  [95] GetoptLong_1.0.5          plotly_4.10.4            
##  [97] remaCor_0.0.18            mime_0.12                
##  [99] splines_4.4.0             circlize_0.4.16          
## [101] fastDummies_1.7.3         dbplyr_2.5.0             
## [103] sparseMatrixStats_1.16.0  knitr_1.46               
## [105] blob_1.2.4                utf8_1.2.4               
## [107] clue_0.3-65               BiocVersion_3.19.1       
## [109] lme4_1.1-35.3             listenv_0.9.1            
## [111] DelayedMatrixStats_1.26.0 Rdpack_2.6               
## [113] tibble_3.2.1              Matrix_1.7-0             
## [115] statmod_1.5.0             svglite_2.1.3            
## [117] fANCOVA_0.6-1             pkgconfig_2.0.3          
## [119] tools_4.4.0               cachem_1.0.8             
## [121] RhpcBLASctl_0.23-42       rbibutils_2.2.16         
## [123] RSQLite_2.3.6             viridisLite_0.4.2        
## [125] DBI_1.2.2                 numDeriv_2016.8-1.1      
## [127] fastmap_1.1.1             rmarkdown_2.26           
## [129] scales_1.3.0              grid_4.4.0               
## [131] ica_1.0-3                 broom_1.0.5              
## [133] AnnotationHub_3.12.0      sass_0.4.9               
## [135] patchwork_1.2.0           BiocManager_1.30.23      
## [137] ggstats_0.6.0             dotCall64_1.1-1          
## [139] RANN_2.6.1                alabaster.schemas_1.4.0  
## [141] farver_2.1.1              aod_1.3.3                
## [143] mgcv_1.9-1                yaml_2.3.8               
## [145] cli_3.6.2                 purrr_1.0.2              
## [147] leiden_0.4.3.1            lifecycle_1.0.4          
## [149] uwot_0.2.2                glmmTMB_1.1.9            
## [151] mvtnorm_1.2-4             backports_1.4.1          
## [153] BiocParallel_1.38.0       gtable_0.3.5             
## [155] rjson_0.2.21              ggridges_0.5.6           
## [157] progressr_0.14.0          limma_3.60.0             
## [159] jsonlite_1.8.8            edgeR_4.2.0              
## [161] RcppHNSW_0.6.0            bitops_1.0-7             
## [163] bit64_4.0.5               Rtsne_0.17               
## [165] alabaster.matrix_1.4.0    spatstat.utils_3.0-4     
## [167] BiocNeighbors_1.22.0      highr_0.10               
## [169] jquerylib_0.1.4           alabaster.se_1.4.0       
## [171] pbkrtest_0.5.2            lazyeval_0.2.2           
## [173] alabaster.base_1.4.1      shiny_1.8.1.1            
## [175] htmltools_0.5.8.1         sctransform_0.4.1        
## [177] rappdirs_0.3.3            glue_1.7.0               
## [179] spam_2.10-0               httr2_1.0.1              
## [181] XVector_0.44.0            gridExtra_2.3            
## [183] EnvStats_2.8.1            boot_1.3-30              
## [185] igraph_2.0.3              variancePartition_1.34.0 
## [187] TMB_1.9.11                R6_2.5.1                 
## [189] tidyr_1.3.1               labeling_0.4.3           
## [191] cluster_2.1.6             Rhdf5lib_1.26.0          
## [193] nloptr_2.0.3              DelayedArray_0.30.1      
## [195] tidyselect_1.2.1          vipor_0.4.7              
## [197] xml2_1.3.6                AnnotationDbi_1.66.0     
## [199] future_1.33.2             rsvd_1.0.5               
## [201] munsell_0.5.1             KernSmooth_2.23-22       
## [203] data.table_1.15.4         htmlwidgets_1.6.4        
## [205] ComplexHeatmap_2.20.0     RColorBrewer_1.1-3       
## [207] rlang_1.1.3               spatstat.sparse_3.0-3    
## [209] spatstat.explore_3.2-7    lmerTest_3.1-3           
## [211] fansi_1.0.6