Macrophage scRNA-seq project: differential expression

Source: TBA

Introduction

macrophage scRNAseq project michelle Wong 2023-07-05

Edits by Mark Ziemann Sept 2023

Purpose of the edits is to output all detected genes from the Seurat differential expression objects for downstream enrichment analysis using the mitch package.

suppressPackageStartupMessages({
  library("Seurat")
  library("SeuratWrappers")
  library("ggplot2")
  library("dplyr")
  library("tidyr")
  library("harmony")
  library("mitch")
  library("kableExtra")
})

Loading processed data

toload<-list.files("counts")
toloadMDM<-grep("^O238|^PV238|^P239", toload, value = TRUE)
processeddata<-list()

for (i in toloadMDM){
sample<-gsub("_filtered_feature_bc_matrix", "", i)
data<-Read10X(paste("counts", i, sep = "/"))
processeddata[[i]]<-CreateSeuratObject(counts = data$`Gene Expression`, project = sample)
processeddata[[i]][["ADT"]]<-CreateAssayObject(counts = data$`Antibody Capture`)
hash<-processeddata[[i]]@assays$ADT@data
hash<-hash[rowMeans(hash)>1,]
print(rownames(hash))
processeddata[[i]][["HTO"]]<-CreateAssayObject(hash)
processeddata[[i]][["HTO"]] <- NormalizeData(processeddata[[i]][["HTO"]], assay = "HTO", normalization.method = "CLR")
processeddata[[i]]<-HTODemux(processeddata[[i]], positive.quantile = 0.95)
processeddata[[i]]$percent.mt<-PercentageFeatureSet(processeddata[[i]], pattern = "^MT-")
processeddata[[i]]<-NormalizeData(processeddata[[i]])
processeddata[[i]]<-FindVariableFeatures(processeddata[[i]], selection.method = "vst", nfeatures = 2000)
processeddata[[i]]<-ScaleData(processeddata[[i]]) 
processeddata[[i]]<-RunPCA(processeddata[[i]])
processeddata[[i]]<-RunUMAP(processeddata[[i]], dims = 1:15)
processeddata[[i]]<-RenameCells(processeddata[[i]], add.cell.id = sample)
processeddata[[i]]<-subset(processeddata[[i]], subset = percent.mt<10 & nCount_RNA <150000 & nFeature_RNA>2000 & HTO_classification.global=="Singlet")
}

## 10X data contains more than one type and is being returned as a list containing matrices of each type.

## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')

## [1] "Mock-MDM"    "GFP-pos-MDM" "GFP-neg-MDM"

## Warning: The following arguments are not used: assay

## Normalizing across features

## Cutoff for Mock-MDM : 120 reads

## Cutoff for GFP-pos-MDM : 200 reads

## Cutoff for GFP-neg-MDM : 5215 reads

## Centering and scaling data matrix

## PC_ 1 
## Positive:  GAL, CAPG, CSTB, SPP1, MT1G, MT2A, MYL9, MMP12, FABP3, MT1H 
##     S100A6, FABP4, CRABP2, UCHL1, IGFL2, TUBB3, PPP1R1A, UTS2, FABP6, ACAT2 
##     WFDC2, MT1X, GREM1, AANAT, MAEL, NRAP, NQO1, CEACAM3, HCG14, ANXA2 
## Negative:  EXOC4, NEAT1, ARL15, DPYD, FHIT, LRMDA, ATG7, FTX, VPS13B, ZEB2 
##     MALAT1, RAD51B, SPIDR, RASAL2, MBD5, GMDS-DT, DOCK3, DOCK4, ELMO1, MAML2 
##     ST18, RERE, ZSWIM6, PAN3, TBC1D5, ETV6, VTI1A, FRMD4B, MEF2A, ARHGAP15 
## PC_ 2 
## Positive:  TM4SF19, CSTB, FDX1, DBI, S100A6, BCL2A1, C15orf48, CSF1, CCL22, ANXA2 
##     CAPG, RGCC, PSMA7, GPC4, SPP1, S100A10, NCAPH, FABP3, LGALS3, ACAT2 
##     CHI3L1, ELOC, MGST1, TXN, TXNRD1, FBP1, FNIP2, CRABP2, MREG, RETREG1 
## Negative:  TGFBI, EPB41L3, RNASE1, CEBPD, SIPA1L1, CD14, HLA-DRA, VMO1, RGS2, HLA-DPB1 
##     TCF4, SSBP2, ARL4C, LILRB2, TIMP1, FOS, HLA-DRB5, GPR155, FPR3, HLA-DQA1 
##     PDE4B, C1QC, MAFB, ST8SIA4, HIF1A, BTG1, MPEG1, HLA-DRB1, LUCAT1, ARMC2 
## PC_ 3 
## Positive:  SRGN, GAPDH, CD74, PFN1, HSP90AA1, HLA-DPA1, HSP90B1, LGMN, CTSC, MS4A7 
##     HLA-DRA, BLVRB, CTSL, HLA-DRB1, PLA2G7, LDHB, ANXA1, MMP9, TXN, GLIPR1 
##     FUCA1, MARCKS, TUBB, TUBA1C, MIF, ACTB, LDHA, HLA-DRB5, S100A10, PSME2 
## Negative:  CT69, AC116534.1, AL035665.2, KCNA2, OSBP2, AC073359.2, AC067751.1, ANKUB1, AL136317.2, LY86-AS1 
##     LINC02320, NOS1AP, TMEM117, CDCP1, AC008771.1, LINC00342, AC078845.1, AC068413.1, AL157886.1, SLC25A48 
##     AC023194.3, LINC01374, AC079142.1, STX18-AS1, CLVS1, AC093334.1, ZC3H12C, KCP, ANKRD36C, GDA 
## PC_ 4 
## Positive:  env, gag, pol, tat, nef, vif, envelope, vpu, eGFP, CCL3 
##     ACTG1, TPM4, vpr, UGCG, ACTB, MARCKS, SPOCD1, ANK2, MACC1, G0S2 
##     PPARG, TCTEX1D2, PCSK6, TNFRSF9, SNCA, CYTOR, TTC39B, SPRED1, TNFSF15, PSAT1 
## Negative:  PTGDS, LINC02244, CLU, APOC1, BX664727.3, ACP5, SYNGR1, CFD, NCAPH, FOLR2 
##     FXYD2, AL136317.2, CPE, RARRES1, PEBP4, AKR7A2, EPHB1, LY86-AS1, S100A8, RAB6B 
##     RCBTB2, COX5B, TMEM114, MT1X, CAMP, MT2A, SPON2, MT1G, MT1E, CSRP2 
## PC_ 5 
## Positive:  SLAMF9, AK8, CP, UTS2, AC078850.1, CTSL, MMP19, PLEKHA5, CLMP, SEL1L2 
##     INSM1, KALRN, BCAT1, SCD, HNRNPLL, CCDC26, GAL, PHLDA1, GPR183, DPP4 
##     SLC16A10, TGM2, OLFML2B, TBC1D8, PPP1R1A, GUCY1A2, MTSS1, VMO1, LINC00862, NMB 
## Negative:  CENPM, TYMS, MYBL2, PCLAF, TK1, MNDA, NUSAP1, NCAPG, SPC25, CENPF 
##     DIAPH3, BIRC5, MKI67, CEP55, ASPM, ANLN, CKS2, HMMR, CDK1, AURKB 
##     RRM2, HDAC9, CENPK, CENPE, TUBA1B, GTSE1, TPX2, DLGAP5, TTK, MT2A

## 14:07:41 UMAP embedding parameters a = 0.9922 b = 1.112

## 14:07:41 Read 2718 rows and found 15 numeric columns

## 14:07:41 Using Annoy for neighbor search, n_neighbors = 30

## 14:07:41 Building Annoy index with metric = cosine, n_trees = 50

## 0%   10   20   30   40   50   60   70   80   90   100%

## [----|----|----|----|----|----|----|----|----|----|

## **************************************************|
## 14:07:42 Writing NN index file to temp file /tmp/RtmpMjxT9O/file355b872c059b5
## 14:07:42 Searching Annoy index using 1 thread, search_k = 3000
## 14:07:42 Annoy recall = 100%
## 14:07:43 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 14:07:43 Initializing from normalized Laplacian + noise (using irlba)
## 14:07:43 Commencing optimization for 500 epochs, with 113234 positive edges
## 14:07:45 Optimization finished
## 10X data contains more than one type and is being returned as a list containing matrices of each type.

## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')

## [1] "Mock-MDM"    "GFP-pos-MDM" "GFP-neg-MDM" "Mock-ALV"    "GFP-pos-ALV"
## [6] "GFP-neg-ALV"

## Warning: The following arguments are not used: assay

## Normalizing across features
## Cutoff for Mock-MDM : 22 reads
## Cutoff for GFP-pos-MDM : 40 reads
## Cutoff for GFP-neg-MDM : 21 reads
## Cutoff for Mock-ALV : 53 reads
## Cutoff for GFP-pos-ALV : 42 reads
## Cutoff for GFP-neg-ALV : 1660 reads
## Centering and scaling data matrix
## PC_ 1 
## Positive:  MALAT1, DOCK3, RASAL2, ARL15, NEAT1, PLXDC2, ASAP1, DPYD, MITF, ATG7 
##     TMEM117, LRMDA, EXOC4, JMJD1C, ZEB2, FTX, CHD9, LRCH1, MEF2A, LPP 
##     ZNF438, DENND4C, MAML2, JARID2, RAPGEF1, FRMD4B, TPRG1, FMNL2, TNS3, ATXN1 
## Negative:  FABP4, NUPR1, STMN1, nef, HAMP, HLA-DQB1, S100A9, SLAMF9, CAMP, GAPDH 
##     HLA-DQA1, TMEM37, LILRA1, TIMP1, BLVRB, GYPC, C1QC, UTS2, RBP1, FXYD2 
##     AC078850.1, DNAJC5B, AC026369.3, IL18BP, VSIG4, S100A8, AL035446.1, MS4A6E, FAH, IL4I1 
## PC_ 2 
## Positive:  MARCH1, VAMP5, RCBTB2, ARHGAP15, FGL2, TRPS1, SAMSN1, SLC9A9, HLA-DPB1, SPRED1 
##     XYLT1, ABCA1, HLA-DPA1, RTN1, HLA-DRA, ME1, TGFBI, KCNMA1, ZNF609, TCF7L2 
##     TSPAN33, CAMK1D, LYZ, TET2, MRC1, FCHSD2, SIPA1L1, ST8SIA4, ATP8B4, SLC8A1 
## Negative:  FDX1, CYSTM1, TM4SF19, CSTB, CSF1, ATP6V1D, CD164, GM2A, PRDX6, GAL 
##     TGM2, ACAT2, SCD, UPP1, CYTOR, BCAT1, ATP6V1G1, FAH, ELOC, LGALS1 
##     TPM4, C15orf48, TXN, GOLGA7B, RHOF, TNFRSF12A, ADCY3, GAPLINC, CCL22, TFRC 
## PC_ 3 
## Positive:  NCAPH, CLU, LINC02244, PTGDS, CRABP2, SYNGR1, LINC01010, TRIM54, DUSP2, AC015660.2 
##     LINC01800, CRIP1, RGS20, CCL22, GPC4, C2orf92, SERPINI1, AC067751.1, STAC, FAIM2 
##     SERTAD2, SPON2, C11orf45, GAL, MT1E, TNFRSF21, MEIKIN, USP2-AS1, CT69, MT2A 
## Negative:  MARCKS, LGMN, AIF1, FUCA1, CTSZ, IFI30, IL18BP, FPR3, HLA-DQA1, HLA-DQB1 
##     BLVRB, CTSL, MS4A7, CD163, PLIN2, C1QC, CD36, FCGR3A, FABP4, IL4I1 
##     CADM1, SLC8A1, GYPC, TMEM37, MS4A6A, TPM4, CTSC, STMN1, C1QA, SMS 
## PC_ 4 
## Positive:  pol, LINC01091, tat, TCTEX1D2, env, nef, gag, FMN1, PRKCE, vif 
##     STOX2, C4orf45, AP000812.1, SLC25A48, ANK2, envelope, ME3, LINC00278, STX18-AS1, AC098829.1 
##     SLC51A, TBC1D8, KCP, vpr, PDE2A, vpu, ANKUB1, FRMD3, BCL2, AC119674.1 
## Negative:  HLA-DRB1, CD74, MMP9, IFI6, LYZ, HLA-DPA1, PTGDS, HLA-DRA, SPP1, ALOX5AP 
##     CFD, CA2, MX1, CTSB, ISG15, TUBA1A, S100A4, MGST1, CYP1B1, IFIT3 
##     HLA-DPB1, VAMP5, CTSZ, S100A6, ATP6V1G1, CHI3L1, RGCC, APLP2, LINC02244, FGL2 
## PC_ 5 
## Positive:  TYMS, PCLAF, TK1, MYBL2, MKI67, RRM2, CLSPN, CENPM, FAM111B, DIAPH3 
##     CDK1, CENPK, ESCO2, NCAPG, CEP55, BIRC5, ATAD2, HELLS, PRC1, CENPF 
##     NUSAP1, DTL, HMMR, CCNA2, KIF11, MCM7, ANLN, MCM10, CIT, TOP2A 
## Negative:  CTSB, SAT1, IFI30, RGCC, IGF2R, C15orf48, ATP6V1G1, PLEK, MMP19, GLRX 
##     CTSZ, env, gag, DCSTAMP, RNF13, NR1H3, CSTB, CD74, tat, nef 
##     PHLDA1, pol, HLA-DRA, LGMN, CD36, RHOF, PLIN2, ATP1B1, HLA-DQA1, HLA-DRB1 
## 14:08:07 UMAP embedding parameters a = 0.9922 b = 1.112
## 14:08:07 Read 9002 rows and found 15 numeric columns
## 14:08:07 Using Annoy for neighbor search, n_neighbors = 30
## 14:08:07 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 14:08:08 Writing NN index file to temp file /tmp/RtmpMjxT9O/file355b86bb6809c
## 14:08:08 Searching Annoy index using 1 thread, search_k = 3000
## 14:08:10 Annoy recall = 100%
## 14:08:10 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 14:08:11 Initializing from normalized Laplacian + noise (using irlba)
## 14:08:11 Commencing optimization for 500 epochs, with 380162 positive edges
## 14:08:17 Optimization finished
## 10X data contains more than one type and is being returned as a list containing matrices of each type.

## Warning: Feature names cannot have underscores ('_'), replacing with dashes
## ('-')

## [1] "Mock-MDM"    "GFP-pos-MDM" "GFP-neg-MDM"

## Warning: The following arguments are not used: assay

## Normalizing across features
## Cutoff for Mock-MDM : 171 reads
## Cutoff for GFP-pos-MDM : 140 reads
## Cutoff for GFP-neg-MDM : 1631 reads
## Centering and scaling data matrix
## PC_ 1 
## Positive:  MT1G, pol, FABP4, nef, MT1H, GAL, UTS2, AC105402.3, MT2A, gag 
##     MMP12, GREM1, MT1X, RAB33A, AC098479.1, AL449043.1, AC091973.1, LINC00460, RUNX1T1, FAM138B 
##     AC090398.2, RAMP3, PRSS35, AL162171.1, AC236972.3, PTN, AC011195.2, LIFR-AS1, AL353764.1, C8A 
## Negative:  ARL15, NEAT1, EXOC4, FTX, FHIT, PDE4D, DPYD, DOCK3, PLXDC2, ARHGAP15 
##     RAD51B, TRIO, MAML2, ZEB2, DENND4C, RASAL2, LRMDA, ZFAND3, FNDC3B, JMJD1C 
##     ATG7, VPS13B, COP1, FRMD4B, TMEM117, JARID2, UBE2E2, MITF, TPRG1, ATXN1 
## PC_ 2 
## Positive:  TGFBI, CD163, HLA-DPB1, ST8SIA4, C1QA, MPEG1, HLA-DQB1, HLA-DRA, HLA-DQA1, MS4A6A 
##     HLA-DQA2, RNASE1, HLA-DPA1, HLA-DRB5, PDE4B, C1QC, TCF4, FPR3, ARL4C, LILRB2 
##     CEBPD, SIPA1L1, SEMA4A, CSF1R, EPB41L3, CD14, CD74, HLA-DOA, SSBP2, MS4A7 
## Negative:  TM4SF19, SPP1, BCL2A1, ANXA2, DBI, GPC4, C15orf48, MGST1, RGCC, ANO5 
##     SNX10, TXN, TXNRD1, ACAT2, CHI3L1, CRABP2, GDF15, PRDX6, EMP1, FNIP2 
##     FBP1, RETREG1, SCD, HES2, S100A9, COX5B, TREM2, NMB, CCL22, MREG 
## PC_ 3 
## Positive:  SRGN, HLA-DRB1, LY96, CD74, HSP90AA1, RPS20, ACTB, HSP90B1, HLA-DPA1, CFD 
##     MARCKS, ANXA1, CTSL, LYZ, BLVRB, AIF1, HLA-DRA, GLIPR1, TXN, TREM2 
##     COX5B, LGMN, DBI, BCL2A1, HLA-DRB5, HLA-DPB1, IGSF6, CTSB, SELENOP, ADAMDEC1 
## Negative:  MIR34AHG, SH3RF3, STOX2, TTC28, PKD1L1, pol, KCNMB2-AS1, TPRG1, AC024084.1, XIST 
##     SLC4A8, LINC02320, KCP, APBB2, RAD51B, STX18-AS1, ST18, PRKCE, TMEM117, FTX 
##     AC079142.1, AC098829.1, LY86-AS1, PDE4D, AL357153.2, AL365295.1, OSBP2, ANK2, SLC16A1-AS1, CDKAL1 
## PC_ 4 
## Positive:  RPS4Y1, PTGDS, TMEM176B, NCAPH, LINC00278, TMEM176A, CLU, RARRES1, TRMT1, CRABP2 
##     UTY, EIF1AY, C11orf45, CCL22, USP9Y, BX664727.3, TTTY14, CHI3L1, S100A4, LINC02244 
##     CRIP1, MNDA, COX5B, AC015660.2, GCLC, CLEC12A, CCND2, RASGRF1, ACAT2, AL136317.2 
## Negative:  SLAMF9, gag, env, XIST, tat, pol, nef, G0S2, LGMN, MARCKS 
##     vif, CLEC4A, CTSB, C5orf17, HES4, SUPV3L1, CCL3, RAB38, MDN1, envelope 
##     STX18-AS1, SDS, PLA2G12A, AMZ1, SNCA, UTS2, BEX2, CTSL, eGFP, NIBAN1 
## PC_ 5 
## Positive:  KCNMA1, XIST, APOC1, RAB6B, PTGDS, MEIKIN, CLU, PCDH19, MMP9, C5orf17 
##     FOLR2, RCBTB2, MT1X, MS4A7, FCGR2B, NMB, FLRT2, BX664727.3, TREM2, LILRA1 
##     PLEKHA5, MT1G, SRGN, STX18-AS1, PEBP4, MT1E, MT1F, KLRD1, MT2A, SYNGR1 
## Negative:  pol, env, gag, tat, nef, RPS4Y1, vif, EIF1AY, TTC39B, LINC00278 
##     envelope, TTTY14, HS3ST2, UTY, USP9Y, ACTB, TNFSF14, UGCG, TMEM176A, AC008050.1 
##     HLA-DQA1, SNCA, ACTG1, vpu, HLA-DQB1, eGFP, TNFSF15, CHI3L1, vpr, TMEM176B 
## 14:08:46 UMAP embedding parameters a = 0.9922 b = 1.112
## 14:08:46 Read 11953 rows and found 15 numeric columns
## 14:08:46 Using Annoy for neighbor search, n_neighbors = 30
## 14:08:46 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 14:08:47 Writing NN index file to temp file /tmp/RtmpMjxT9O/file355b81510150d
## 14:08:47 Searching Annoy index using 1 thread, search_k = 3000
## 14:08:50 Annoy recall = 100%
## 14:08:50 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 14:08:51 Initializing from normalized Laplacian + noise (using irlba)
## 14:08:51 Commencing optimization for 200 epochs, with 515920 positive edges
## 14:08:54 Optimization finished

Integration and metadata

processeddata

## $O238_filtered_feature_bc_matrix
## An object of class Seurat 
## 36619 features across 1244 samples within 3 assays 
## Active assay: RNA (36613 features, 2000 variable features)
##  2 other assays present: ADT, HTO
##  2 dimensional reductions calculated: pca, umap
## 
## $P239_filtered_feature_bc_matrix
## An object of class Seurat 
## 36625 features across 5536 samples within 3 assays 
## Active assay: RNA (36613 features, 2000 variable features)
##  2 other assays present: ADT, HTO
##  2 dimensional reductions calculated: pca, umap
## 
## $PV238_filtered_feature_bc_matrix
## An object of class Seurat 
## 36619 features across 4914 samples within 3 assays 
## Active assay: RNA (36613 features, 2000 variable features)
##  2 other assays present: ADT, HTO
##  2 dimensional reductions calculated: pca, umap

#demultiplex P and V (donors 1/2)
meta<-processeddata[[3]]@meta.data
split<- read.csv("scSplit_result.csv") #note, this was not run on as many cells? may need to redo

library(stringr)
split[c('cell.bc', 'SNP')]<- str_split_fixed(split$Barcode.Cluster,"-1",2)
split$cell.bc <- paste("PV238", split$cell.bc, sep="_")
split$cell.bc <- paste(split$cell.bc, "-1", sep="")
split<-split[,2:3]

meta$cell.bc<-rownames(meta)
meta2<-merge(meta,split, by = "cell.bc", all.x = TRUE)
processeddata[[3]]@meta.data$SNP<-meta2$SNP

rm(data)
rm(hash)
rm(meta)
rm(meta2)
rm(split)

table(processeddata[[3]]$orig.ident)

## 
## PV238 
##  4914

Combine samples

combined <- merge(processeddata[[1]], processeddata[[2]])
for(i in 3:length(processeddata)){combined<-merge(combined, processeddata[[i]])}

combined

## An object of class Seurat 
## 36625 features across 11694 samples within 3 assays 
## Active assay: RNA (36613 features, 0 variable features)
##  2 other assays present: ADT, HTO

#demultiplex MDM and ALV for P239

combined$macro_type<-case_when(grepl("ALV$", combined$HTO_classification)~"ALV",
                               grepl("MDM$", combined$HTO_classification)~"MDM")

table(combined$macro_type)

## 
##  ALV  MDM 
## 3300 8394

combined<-subset(combined, macro_type == "MDM")

combined

## An object of class Seurat 
## 36625 features across 8394 samples within 3 assays 
## Active assay: RNA (36613 features, 0 variable features)
##  2 other assays present: ADT, HTO

####integration
DefaultAssay(combined) <- "RNA"
combined <- FindVariableFeatures(combined, selection.method = "vst", nfeatures = 2000)

features_var<-VariableFeatures(combined)

str(features_var)

##  chr [1:2000] "MT1G" "MT1H" "MT1X" "CAMP" "pol" "MT1M" "MT2A" "UTS2" "env" ...

features_var <- grep('pol|^gag|^env|^tat|^vif|^nef|^envelope|^vpu|^eGFP', features_var, value = TRUE, invert = TRUE)
VariableFeatures <- features_var

combined <- ScaleData(combined, features = features_var)

## Centering and scaling data matrix

combined <- RunPCA(combined, npcs = 15, features = features_var)

## PC_ 1 
## Positive:  S100A10, TXN, PRDX6, COX5B, ACP5, C15orf48, FABP3, BCL2A1, PSME2, TUBA1B 
##     CALM3, MGST1, ACTB, SPP1, CHI3L1, ACAT2, TUBA1A, ACTG1, MMP9, CTSL 
##     CRABP2, ANXA1, GSTM4, IFI30, TUBA1C, FAM89B, UCHL1, HAMP, LILRA1, FBP1 
## Negative:  ARL15, FTX, EXOC4, NEAT1, DPYD, PLXDC2, VPS13B, FHIT, MALAT1, RAD51B 
##     ZEB2, JMJD1C, ZFAND3, LRMDA, MBD5, TRIO, SPIDR, RASAL2, TCF12, DOCK3 
##     COP1, DOCK4, SBF2, ZNF438, ZSWIM6, ARHGAP15, ATG7, MAML2, ELMO1, TMEM117 
## PC_ 2 
## Positive:  TM4SF19, ANO5, GPC4, GM2A, TXNRD1, FNIP2, BCL2A1, PSD3, RETREG1, TXN 
##     SLC28A3, RGS20, SPP1, CALM3, TGM2, TCTEX1D2, LINC01010, BCAT1, EPB41L1, CCL22 
##     FABP3, CCDC26, SPSB1, NIBAN1, RGCC, LINC01135, AC092353.2, MGST1, FKBP15, SCD 
## Negative:  CD74, HLA-DRA, HLA-DPB1, HLA-DPA1, AIF1, TGFBI, HLA-DRB1, C1QA, HLA-DQB1, HLA-DQA1 
##     FPR3, MS4A6A, C1QC, MS4A7, LYZ, CD14, CD163, CEBPD, MPEG1, TIMP1 
##     FOS, MAFB, CSF1R, GLUL, ST8SIA4, EPB41L3, IFI30, TCF4, FCN1, HLA-DOA 
## PC_ 3 
## Positive:  ACTB, CALR, FBP1, ACTG1, SLC35F1, TUBA1B, TIMP3, MGLL, PLEK, HSP90B1 
##     LINC01091, GSN, IL1RN, LPL, IGSF6, GCLC, RGCC, LHFPL2, TMEM176B, INSIG1 
##     MADD, PDIA4, CSF1, LINC00278, HLA-DRB1, HSPA5, TNS3, CHI3L1, RPS4Y1, MREG 
## Negative:  CCPG1, HES2, NUPR1, CARD16, BTG1, NMB, PSAT1, CLGN, S100P, PLBD1 
##     STMN1, G0S2, PLEKHA5, SUPV3L1, CLEC4E, CXCR4, NIBAN1, IFI6, DUSP1, SEL1L3 
##     XIST, BEX2, RAB6B, ME1, CLEC4A, SLAMF9, TBL1X, VWA5A, PDE4D, QPCT 
## PC_ 4 
## Positive:  PTGDS, NCAPH, HLA-C, RNASE6, CLU, BX664727.3, JUNB, CCL22, SYNGR1, COX5B 
##     ADRA2B, AKR7A2, CSF3R, CAMP, CRABP2, CCND2, RCBTB2, AL136317.2, TSPAN33, LINC02244 
##     ZFP36L1, SSBP3, FXYD2, TMEM176B, PEBP4, CPE, BASP1, RASSF4, LY86-AS1, NIPAL2 
## Negative:  MARCKS, UGCG, TPM4, CLEC4A, SNCA, SUPV3L1, CD36, B4GALT5, CCL3, RAB38 
##     BCAT1, GDF15, NIBAN1, PLAU, PSAT1, ATP13A3, PPARG, PCSK6, CLGN, MDN1 
##     CARD16, BEX2, SLAMF9, CTSL, TM4SF19, ANXA1, G0S2, PAG1, ANK2, LGMN 
## PC_ 5 
## Positive:  RPS4Y1, CLEC12A, LYZ, MNDA, MT-CO3, MT-ATP6, LINC00278, UTY, MT-ND1, MT-CYB 
##     CHI3L1, USP9Y, TTTY14, MT-ND2, COX5B, ALOX5AP, TNFSF13B, PRSS41, MT-ATP8, S100P 
##     IGSF6, BCL2A1, P2RY13, CLGN, CAMK1D, PSAT1, NUPR1, LINC01500, LINC01515, FBP1 
## Negative:  MMP19, AC078850.1, TFRC, CTSL, LGMN, SLC7A8, GAL, AK8, SEL1L2, UTS2 
##     PHLDA1, IVNS1ABP, SDS, SCD, INSIG1, GPR183, FABP4, SLAMF9, MMP9, IL18BP 
##     CSF1, BCAT1, HNRNPLL, TBC1D8, TGM2, MTSS1, SRGN, CP, FCGR3A, CD300LB

combined <- RunHarmony(combined, group.by.vars="orig.ident", assay.use = "RNA", reduction = "pca", dims.use = 1:15, plot_convergence = TRUE)

## Harmony 1/10

## Harmony 2/10

## Harmony 3/10

## Harmony 4/10

## Harmony 5/10

## Harmony converged after 5 iterations

## Warning: Invalid name supplied, making object name syntactically valid. New
## object name is Seurat..ProjectDim.RNA.harmony; see ?make.names for more details
## on syntax validity

combined <- RunUMAP(combined, reduction = "harmony", reduction.name = "umap.rna", reduction.key = "rnaUMAP_", dims = 1:15, assay = "RNA")

## 14:09:07 UMAP embedding parameters a = 0.9922 b = 1.112

## 14:09:07 Read 8394 rows and found 15 numeric columns

## 14:09:07 Using Annoy for neighbor search, n_neighbors = 30

## 14:09:07 Building Annoy index with metric = cosine, n_trees = 50

## 0%   10   20   30   40   50   60   70   80   90   100%

## [----|----|----|----|----|----|----|----|----|----|

## **************************************************|
## 14:09:08 Writing NN index file to temp file /tmp/RtmpMjxT9O/file355b8491e26a1
## 14:09:08 Searching Annoy index using 1 thread, search_k = 3000
## 14:09:10 Annoy recall = 100%
## 14:09:10 Commencing smooth kNN distance calibration using 1 thread with target n_neighbors = 30
## 14:09:11 Initializing from normalized Laplacian + noise (using irlba)
## 14:09:11 Commencing optimization for 500 epochs, with 357442 positive edges
## 14:09:17 Optimization finished

combined <- FindNeighbors(combined, reduction = "harmony", dims = 1:15)

## Computing nearest neighbor graph
## Computing SNN

combined <- FindClusters(combined, resolution = 0.4)

## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 8394
## Number of edges: 264586
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8207
## Number of communities: 6
## Elapsed time: 0 seconds

DimPlot(combined, group.by = "HTO_classification")

Metadata

combined$donor<-case_when(grepl("\tSNG-1",combined$SNP)~"P238",
                          grepl("\tSNG-2", combined$SNP)~"V238",
                          grepl("O238", combined$orig.ident)~"O238",
                          grepl("P239", combined$orig.ident)~"P239")
table(combined$donor)

## 
## O238 P238 P239 V238 
## 1244 1815 2236 2901

combined$HIV_protein<-combined$HTO_classification
HIV_genes<-c("pol", "gag", "env", "tat", "vif","nef","envelope", "vpu","eGFP")
VlnPlot(combined, HIV_genes, group.by = "HIV_protein")

#adding a HIV_genes collection
RNA_counts<-combined@assays$RNA@counts
str(RNA_counts)

## Formal class 'dgCMatrix' [package "Matrix"] with 6 slots
##   ..@ i       : int [1:42731828] 44 48 59 60 62 69 73 77 83 86 ...
##   ..@ p       : int [1:8395] 0 4365 11361 16864 21790 27005 33785 39966 46888 53465 ...
##   ..@ Dim     : int [1:2] 36613 8394
##   ..@ Dimnames:List of 2
##   .. ..$ : chr [1:36613] "MIR1302-2HG" "FAM138A" "OR4F5" "AL627309.1" ...
##   .. ..$ : chr [1:8394] "O238_AAACCCAAGTGCAACG-1" "O238_AAACGCTCACTGGACC-1" "O238_AAACGCTCATGTTACG-1" "O238_AAACGCTGTCGCTGCA-1" ...
##   ..@ x       : num [1:42731828] 9 2 6 2 7 1 4 2 2 5 ...
##   ..@ factors : list()

RNA_counts_HIV<-RNA_counts[rownames(RNA_counts) %in% HIV_genes, ]

rownames(RNA_counts_HIV)

## [1] "env"      "gag"      "pol"      "vif"      "tat"      "vpu"      "envelope"
## [8] "nef"      "eGFP"

RNA_counts_HIV<-as.data.frame(RNA_counts_HIV)
RNA_counts_HIV<-t(RNA_counts_HIV)
RNA_counts_HIV<-as.data.frame(RNA_counts_HIV)
RNA_counts_HIV$HIV<-rowSums(RNA_counts_HIV)

combined$HIV_counts<-RNA_counts_HIV$HIV
meta<-combined@meta.data
#meta


VlnPlot(combined, "HIV_counts", group.by = "HIV_protein", y.max = 20)

## Warning: Removed 1938 rows containing non-finite values (`stat_ydensity()`).

## Warning: Removed 1938 rows containing missing values (`geom_point()`).

combined$HIV_RNA<-case_when(combined$HIV_counts>5~"positive",
                            combined$HIV_counts<=5~"negative")


combined$HIV<-case_when(grepl("GFP-pos-MDM",combined$HIV_protein)~"productive", 
                        grepl("GFP-neg-MDM", combined$HIV_protein)&grepl("positive", combined$HIV_RNA)~"latent",
                        grepl("GFP-neg-MDM", combined$HIV_protein)~"bystander",
                        grepl("Mock-MDM", combined$HIV_protein)~"mock-infected"
)

table(combined$HIV, combined$donor)

##                
##                 O238 P238 P239 V238
##   bystander      467  771 1413 1281
##   latent         207  161   30  125
##   mock-infected  185  422  447  762
##   productive     385  461  346  733

Annotation

DimPlot(combined)

combined.markers <- FindAllMarkers(combined,features=features_var, only.pos = TRUE,
  min.pct = 0.25, logfc.threshold = 0.25)

## Calculating cluster 0

## Calculating cluster 1

## Calculating cluster 2

## Calculating cluster 3

## Calculating cluster 4

## Calculating cluster 5

combined.markers2<-combined.markers %>%
  group_by(cluster) %>%
  slice_max(n = 10, order_by = avg_log2FC)

write.csv(combined.markers, "cluster_markers.csv")

saveRDS(combined, "MDM_object.RDS")

Analysis

Idents(combined)<-combined$HIV_protein
levels(combined)

## [1] "GFP-neg-MDM" "GFP-pos-MDM" "Mock-MDM"

levels(combined)<-c("Mock-MDM", "GFP-neg-MDM", "GFP-pos-MDM")

pdf("HIV_sort_group.pdf", width = 9, height = 27)
VlnPlot(combined, "HIV_counts")
VlnPlot(combined, HIV_genes)
dev.off()

## png 
##   2

VlnPlot(combined, "HIV_counts", log = TRUE)

VlnPlot(combined, "HIV_counts")

VlnPlot(combined, HIV_genes)

Idents(combined)<-combined$HIV
levels(combined)

## [1] "bystander"     "productive"    "mock-infected" "latent"

levels(combined) <- c("mock-infected", "bystander", "latent", "productive")
VlnPlot(combined, "HIV_counts")

VlnPlot(combined, "HIV_counts", log = TRUE)

VlnPlot(combined, HIV_genes)

VlnPlot(combined, HIV_genes, log = TRUE)

FeaturePlot(combined, "HIV_counts", min.cutoff = 5, max.cutoff = 2500)

DimPlot(combined, cells.highlight=grep("latent",combined$HIV), sizes.highlight = 0.5)+NoLegend()

DimPlot(combined, cells.highlight=grep("productive",combined$HIV), sizes.highlight = 0.5)+NoLegend()

DimPlot(combined, group.by = "seurat_clusters")

DEG analysis

In the following, “ident1” is considered the case, and “ident2” is considered the control.

Here are the comparisons:

1. Latent vs bystander

2. Latent vs productive

3. Productive vs all

4. Productive vs mock

5. Latent vs mock

DE1: latent vs bystander

Idents(combined)<-combined$HIV
levels(combined) <- c("mock-infected", "bystander", "latent", "productive")
saveRDS(combined, "MDM_object.RDS")

latent_v_bystander <- FindMarkers(combined, ident.1 = "latent", ident.2 = "bystander",
  logfc.threshold=0, only.pos=FALSE, group.by=NULL, subset.ident=NULL,
  assay=NULL, reduction=NULL, features=NULL,test.use="wilcox",
  min.pct=0.1, min.diff.pct=-Inf,verbose=TRUE,max.cells.per.ident=Inf,
  random.seed=1,latent.vars.=NULL,min.cells.feature=3, pseudocount.use=1)

dim(latent_v_bystander)

## [1] 10939     5

head(latent_v_bystander,30) %>%
  kbl(caption="Latent vs bystander markers DEG table") %>%
  kable_paper("hover", full_width = FALSE)

Latent vs bystander markers DEG table

	p_val	avg_log2FC	pct.1	pct.2	p_val_adj
pol	0	0.7774630	0.987	0.327	0
env	0	0.3710173	0.599	0.121	0
tat	0	0.2716689	0.434	0.060	0
gag	0	0.3160732	0.484	0.090	0
nef	0	0.3526073	0.582	0.155	0
vif	0	0.1121733	0.198	0.018	0
TMSB10	0	-0.1369325	1.000	1.000	0
RNASEK	0	-0.1501834	0.998	1.000	0
envelope	0	0.0733792	0.133	0.008	0
ACTB	0	-0.3452073	1.000	0.999	0
COX8A	0	-0.2021623	0.991	0.997	0
RPS28	0	-0.1072852	1.000	1.000	0
CD68	0	-0.1384003	0.998	0.999	0
RPL36	0	-0.1339130	0.998	1.000	0
ATP5MPL	0	-0.1892903	0.983	0.998	0
ATP5MF	0	-0.2094586	0.982	0.994	0
ATP5F1E	0	-0.1040254	1.000	1.000	0
ATP6V0C	0	-0.1511387	0.996	0.999	0
CTSD	0	-0.1081730	1.000	1.000	0
FCER1G	0	-0.1522394	1.000	1.000	0
PFN1	0	-0.1902416	0.996	0.999	0
CLIC1	0	-0.2141131	0.972	0.994	0
COX6C	0	-0.1775594	0.987	0.998	0
GAPDH	0	-0.1800718	0.996	0.998	0
ACTG1	0	-0.2913335	0.982	0.989	0
CYBA	0	-0.1141666	1.000	1.000	0
RPL35	0	-0.1303072	1.000	1.000	0
NDUFS5	0	-0.1578098	0.993	0.997	0
ATOX1	0	-0.1637638	0.993	0.999	0
RPS29	0	-0.1649383	1.000	1.000	0

write.csv(latent_v_bystander, "latent_vs_bystander_markers.csv")

DE2: latent vs productive

latent_v_productive <- FindMarkers(combined, ident.1 = "latent", ident.2 = "productive",
  logfc.threshold=0, only.pos=FALSE, group.by=NULL, subset.ident=NULL,
  assay=NULL, reduction=NULL, features=NULL,test.use="wilcox",
  min.pct=0.1, min.diff.pct=-Inf,verbose=TRUE,max.cells.per.ident=Inf,
  random.seed=1,latent.vars.=NULL,min.cells.feature=3, pseudocount.use=1)

dim(latent_v_productive)

## [1] 11129     5

head(latent_v_productive,30) %>%
  kbl(caption="Latent vs productive markers DEG table") %>%
  kable_paper("hover", full_width = FALSE)

Latent vs productive markers DEG table

	p_val	avg_log2FC	pct.1	pct.2	p_val_adj
gag	0	-1.3677864	0.484	0.956	0
nef	0	-1.0716866	0.582	0.953	0
env	0	-1.2486529	0.599	0.946	0
pol	0	-1.0877397	0.987	0.971	0
tat	0	-1.0329311	0.434	0.922	0
vif	0	-0.6945865	0.198	0.800	0
envelope	0	-0.4987595	0.133	0.741	0
vpu	0	-0.3640407	0.061	0.647	0
ARPC1B	0	-0.1922572	1.000	1.000	0
FTL	0	-0.0798613	1.000	1.000	0
eGFP	0	-0.3023021	0.065	0.603	0
RPLP1	0	-0.1182027	1.000	1.000	0
CTSS	0	0.1944087	1.000	1.000	0
BLVRB	0	-0.3279020	0.980	0.997	0
MIF	0	-0.3357282	0.950	0.991	0
GPNMB	0	0.1367290	1.000	1.000	0
CD44	0	0.3068821	0.993	0.942	0
ZFP36L1	0	0.3656252	0.989	0.932	0
RPS9	0	0.1429649	1.000	1.000	0
PDXK	0	-0.2313934	0.965	0.993	0
PTGDS	0	0.6787985	0.945	0.882	0
JUNB	0	0.3565184	0.943	0.879	0
MS4A6A	0	0.3248831	0.543	0.176	0
TTC3	0	0.2519835	0.985	0.970	0
vpr	0	-0.2167249	0.055	0.498	0
CFD	0	0.4215147	0.939	0.771	0
CSTB	0	-0.1404606	1.000	1.000	0
YBX1	0	-0.1381850	1.000	1.000	0
EEF1A1	0	0.1190320	1.000	1.000	0
RPS27L	0	-0.1526173	0.996	0.999	0

write.csv(latent_v_productive, "latent_vs_productive_markers.csv")

DE3: Productive vs all

productive_v_all <- FindMarkers(combined, ident.1 = "productive",
  logfc.threshold=0, only.pos=FALSE, group.by=NULL, subset.ident=NULL,
  assay=NULL, reduction=NULL, features=NULL,test.use="wilcox",
  min.pct=0.1, min.diff.pct=-Inf,verbose=TRUE,max.cells.per.ident=Inf,
  random.seed=1,latent.vars.=NULL,min.cells.feature=3, pseudocount.use=1)

dim(productive_v_all)

## [1] 10703     5

head(productive_v_all,20) %>%
  kbl(caption="productive vs all others markers DEG table") %>%
  kable_paper("hover", full_width = FALSE)

productive vs all others markers DEG table

	p_val	avg_log2FC	pct.1	pct.2	p_val_adj
FTL	0	0.0807056	1.000	1.000	0
RPS9	0	-0.1934307	1.000	1.000	0
CST3	0	-0.2091869	1.000	1.000	0
env	0	1.5865057	0.946	0.156	0
gag	0	1.6544926	0.956	0.124	0
pol	0	1.7865421	0.971	0.366	0
vif	0	0.7968048	0.800	0.033	0
vpr	0	0.2432598	0.498	0.009	0
tat	0	1.2824556	0.922	0.086	0
vpu	0	0.3925134	0.647	0.016	0
envelope	0	0.5656868	0.741	0.019	0
nef	0	1.4018344	0.953	0.170	0
eGFP	0	0.3379257	0.603	0.010	0
GPNMB	0	-0.1293370	1.000	1.000	0
RPL10A	0	-0.2556359	0.993	0.995	0
PTMA	0	-0.1928794	0.999	0.999	0
RPL10	0	-0.1210226	0.999	1.000	0
VIM	0	-0.1547093	1.000	1.000	0
ARPC1B	0	0.1271721	1.000	0.999	0
RPL34	0	-0.1305494	0.999	1.000	0

write.csv(productive_v_all, "productive_vs_all_markers.csv")

DE4: Productive vs mock

productive_v_mock <-FindMarkers(combined, ident.1 = "productive", ident.2= "mock-infected",
  logfc.threshold=0, only.pos=FALSE, group.by=NULL, subset.ident=NULL,
  assay=NULL, reduction=NULL, features=NULL,test.use="wilcox",
  min.pct=0.1, min.diff.pct=-Inf,verbose=TRUE,max.cells.per.ident=Inf,
  random.seed=1,latent.vars.=NULL,min.cells.feature=3, pseudocount.use=1)

dim(productive_v_mock)

## [1] 10713     5

head(productive_v_mock,20) %>%
  kbl(caption="productive vs mock markers DEG table") %>%
  kable_paper("hover", full_width = FALSE)

productive vs mock markers DEG table

	p_val	avg_log2FC	pct.1	pct.2	p_val_adj
env	0	1.6268460	0.946	0.101	0
gag	0	1.6841584	0.956	0.092	0
pol	0	1.8857008	0.971	0.269	0
vif	0	0.8061832	0.800	0.019	0
tat	0	1.3144762	0.922	0.039	0
vpu	0	0.3938672	0.647	0.015	0
envelope	0	0.5716983	0.741	0.009	0
nef	0	1.4627098	0.953	0.082	0
eGFP	0	0.3418937	0.603	0.004	0
RPS9	0	-0.1816366	1.000	0.999	0
vpr	0	0.2467055	0.498	0.003	0
FTL	0	0.0830717	1.000	1.000	0
PTMA	0	-0.2128902	0.999	0.998	0
CST3	0	-0.1912393	1.000	0.999	0
GPNMB	0	-0.1304474	1.000	1.000	0
VIM	0	-0.1636923	1.000	0.999	0
TMSB4X	0	-0.1142752	1.000	0.999	0
RPL10A	0	-0.2426967	0.993	0.992	0
RPS19	0	0.0925151	1.000	0.999	0
ARPC1B	0	0.1322292	1.000	0.998	0

write.csv(productive_v_mock, "productive_vs_mock_markers.csv")

DE5: latent vs mock

latent_v_mock <- FindMarkers(combined, ident.1 = "latent", ident.2= "mock-infected",
  logfc.threshold=0, only.pos=FALSE, group.by=NULL, subset.ident=NULL,
  assay=NULL, reduction=NULL, features=NULL,test.use="wilcox",
  min.pct=0.1, min.diff.pct=-Inf,verbose=TRUE,max.cells.per.ident=Inf,
  random.seed=1,latent.vars.=NULL,min.cells.feature=3, pseudocount.use=1)

dim(latent_v_mock)

## [1] 10951     5

head(latent_v_mock,20) %>%
  kbl(caption="latent vs mock markers DEG table") %>%
  kable_paper("hover", full_width = FALSE)

latent vs mock markers DEG table

	p_val	avg_log2FC	pct.1	pct.2	p_val_adj
pol	0	0.7979611	0.987	0.269	0
nef	0	0.3910233	0.582	0.082	0
env	0	0.3781931	0.599	0.101	0
tat	0	0.2815451	0.434	0.039	0
gag	0	0.3163720	0.484	0.092	0
ACTB	0	-0.3457681	1.000	0.998	0
IFI6	0	0.2834410	0.993	0.982	0
FCER1G	0	-0.1686309	1.000	1.000	0
vif	0	0.1115967	0.198	0.019	0
TMSB10	0	-0.1103753	1.000	1.000	0
ACTG1	0	-0.3171860	0.982	0.986	0
NPC2	0	-0.1153591	1.000	0.999	0
COX8A	0	-0.1749063	0.991	0.994	0
PFN1	0	-0.1820098	0.996	0.998	0
RNASEK	0	-0.1139961	0.998	0.999	0
CD52	0	-0.1679406	0.969	0.996	0
envelope	0	0.0729387	0.133	0.009	0
CTSD	0	-0.0999097	1.000	1.000	0
ATP5MPL	0	-0.1658110	0.983	0.996	0
PTMA	0	-0.1199486	1.000	0.998	0

write.csv(latent_v_mock, "latent_v_mock_markers.csv")

Stackplots

meta<-combined@meta.data
productive <-meta %>% filter(HIV == "productive") %>% group_by(seurat_clusters) %>% summarise(n())
latent<-meta %>% filter(HIV == "latent") %>% group_by(seurat_clusters) %>% summarise(n())
HIV_by_cluster <-meta %>% group_by(HIV, seurat_clusters) %>% summarise(n())

## `summarise()` has grouped output by 'HIV'. You can override using the
## `.groups` argument.

write.csv(HIV_by_cluster, "cluster_composition.csv")

table(combined$HIV, combined$seurat_clusters)

##                
##                    0    1    2    3    4    5
##   bystander     1033 1319  751  617  242   58
##   latent         224  114   77   46   78    2
##   mock-infected  541  502  333  295  161   32
##   productive     697  461  482  294   33    2

Session information

For reproducibility

save.image("processing.Rdata")

sessionInfo()

## R version 4.3.1 (2023-06-16)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 22.04.3 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] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] stringr_1.5.0        harmony_0.1.1        Rcpp_1.0.11         
##  [4] tidyr_1.3.0          SeuratWrappers_0.3.1 SeuratObject_4.1.3  
##  [7] Seurat_4.3.0.1       pkgload_1.3.2.1      GGally_2.1.2        
## [10] ggplot2_3.4.3        reshape2_1.4.4       beeswarm_0.4.0      
## [13] gtools_3.9.4         tibble_3.2.1         dplyr_1.1.3         
## [16] echarts4r_0.4.5      kableExtra_1.3.4     gplots_3.1.3        
## [19] mitch_1.12.0        
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3     rstudioapi_0.15.0      jsonlite_1.8.7        
##   [4] magrittr_2.0.3         spatstat.utils_3.0-3   farver_2.1.1          
##   [7] rmarkdown_2.25         vctrs_0.6.3            ROCR_1.0-11           
##  [10] spatstat.explore_3.2-3 webshot_0.5.5          htmltools_0.5.6       
##  [13] sass_0.4.7             sctransform_0.3.5      parallelly_1.36.0     
##  [16] KernSmooth_2.23-22     bslib_0.5.1            htmlwidgets_1.6.2     
##  [19] ica_1.0-3              plyr_1.8.8             plotly_4.10.2         
##  [22] zoo_1.8-12             cachem_1.0.8           igraph_1.5.1          
##  [25] mime_0.12              lifecycle_1.0.3        pkgconfig_2.0.3       
##  [28] rsvd_1.0.5             Matrix_1.6-1.1         R6_2.5.1              
##  [31] fastmap_1.1.1          fitdistrplus_1.1-11    future_1.33.0         
##  [34] shiny_1.7.5            digest_0.6.33          colorspace_2.1-0      
##  [37] reshape_0.8.9          patchwork_1.1.3        tensor_1.5            
##  [40] irlba_2.3.5.1          labeling_0.4.3         progressr_0.14.0      
##  [43] fansi_1.0.4            spatstat.sparse_3.0-2  httr_1.4.7            
##  [46] polyclip_1.10-4        abind_1.4-5            compiler_4.3.1        
##  [49] remotes_2.4.2.1        withr_2.5.0            highr_0.10            
##  [52] R.utils_2.12.2         MASS_7.3-60            caTools_1.18.2        
##  [55] tools_4.3.1            lmtest_0.9-40          httpuv_1.6.11         
##  [58] future.apply_1.11.0    goftest_1.2-3          R.oo_1.25.0           
##  [61] glue_1.6.2             nlme_3.1-163           promises_1.2.1        
##  [64] grid_4.3.1             Rtsne_0.16             cluster_2.1.4         
##  [67] generics_0.1.3         gtable_0.3.4           spatstat.data_3.0-1   
##  [70] R.methodsS3_1.8.2      data.table_1.14.8      sp_2.0-0              
##  [73] xml2_1.3.5             utf8_1.2.3             spatstat.geom_3.2-5   
##  [76] RcppAnnoy_0.0.21       ggrepel_0.9.3          RANN_2.6.1            
##  [79] pillar_1.9.0           later_1.3.1            splines_4.3.1         
##  [82] lattice_0.21-8         survival_3.5-7         deldir_1.0-9          
##  [85] tidyselect_1.2.0       miniUI_0.1.1.1         pbapply_1.7-2         
##  [88] knitr_1.44             gridExtra_2.3          svglite_2.1.1         
##  [91] scattermore_1.2        xfun_0.40              matrixStats_1.0.0     
##  [94] stringi_1.7.12         lazyeval_0.2.2         yaml_2.3.7            
##  [97] evaluate_0.21          codetools_0.2-19       BiocManager_1.30.22   
## [100] cli_3.6.1              uwot_0.1.16            xtable_1.8-4          
## [103] reticulate_1.32.0      systemfonts_1.0.4      munsell_0.5.0         
## [106] jquerylib_0.1.4        globals_0.16.2         spatstat.random_3.1-6 
## [109] png_0.1-8              parallel_4.3.1         ellipsis_0.3.2        
## [112] bitops_1.0-7           listenv_0.9.0          viridisLite_0.4.2     
## [115] scales_1.2.1           ggridges_0.5.4         leiden_0.4.3          
## [118] purrr_1.0.2            rlang_1.1.1            cowplot_1.1.1         
## [121] rvest_1.0.3