How to treat zeros and NA values from normalized Anndata object?

[estebanelias]

R version 4.4.0 (2024-04-24) -- "Puppy Cup"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: aarch64-apple-darwin20

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

Natural language support but running in an English locale

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

[Workspace loaded from ~/.RData]

Loading required package: SeuratObject
Loading required package: sp

Attaching package: ‘SeuratObject’

The following objects are masked from ‘package:base’:

intersect, t

library(anndata)

Attaching package: ‘anndata’

The following object is masked from ‘package:SeuratObject’:

Layers

library(patchwork)
library(CellChat)
Loading required package: dplyr

Attaching package: ‘dplyr’

The following objects are masked from ‘package:stats’:

filter, lag

The following objects are masked from ‘package:base’:

intersect, setdiff, setequal, union

Loading required package: igraph

Attaching package: ‘igraph’

The following objects are masked from ‘package:dplyr’:

as_data_frame, groups, union

The following objects are masked from ‘package:stats’:

decompose, spectrum

The following object is masked from ‘package:base’:

union

Loading required package: ggplot2

read the data .h5ad

ad <- read_h5ad("/Users/estebanelias/Library/CloudStorage/OneDrive-UniversityofCalgary/Ucalgary/Single cell/cell_communication/adata_IRI_D1_annotated.h5ad")
Would you like to create a default Python environment for the reticulate package? (Yes/no/cancel) gc()
Error in utils::askYesNo(sprintf("Would you like to create a default Python environment for the %s package?", :
Unrecognized response “gc()”
gc()
used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
Ncells 5236976 279.7 9082685 485.1 NA 8113588 433.4
Vcells 9546976 72.9 1522930256 11619.1 16384 1591894266 12145.2
library(anndata)
library(patchwork)
library(CellChat)

read the data .h5ad

ad <- read_h5ad("/Users/estebanelias/Library/CloudStorage/OneDrive-UniversityofCalgary/Ucalgary/Single cell/cell_communication/adata_IRI_D1_annotated.h5ad")
Would you like to create a default Python environment for the reticulate package? (Yes/no/cancel) no

access count data matrix

counts <- t(as.matrix(ad$X))
#Access metadata
meta <- ad$obs
meta$labels <- meta[["cell_id"]]
cellchat <- createCellChat(object = counts, meta = meta, group.by = "labels")
[1] "Create a CellChat object from a data matrix"
Set cell identities for the new CellChat object
The cell groups used for CellChat analysis are B cell DC DTC EC Granulocyte Gzmb+ cell LOH Macrophage Monocyte NK cell PCT PST Podocyte Stroma T cell cDC1 cDC2 pDC
#Set the ligand receptor interaction database
CellChatDB <- CellChatDB.mouse
showDatabaseCategory(CellChatDB)
#Show the structure of the database
dplyr::glimpse(CellChatDB$interaction)
Rows: 2,019
Columns: 11
$ interaction_name "TGFB1_TGFBR1_TGFBR2", "TGFB2_TGFBR1_TGFBR2", "TGFB3_TGFBR1_TGFBR2", "TGFB1_ACVR1B_TGFBR2", "TGFB1_ACV…
$ pathway_name "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "TGFb", "BMP",…
$ ligand "Tgfb1", "Tgfb2", "Tgfb3", "Tgfb1", "Tgfb1", "Tgfb2", "Tgfb2", "Tgfb3", "Tgfb3", "Tgfb1", "Tgfb2", "Tg…
$ receptor "TGFbR1_R2", "TGFbR1_R2", "TGFbR1_R2", "ACVR1B_TGFbR2", "ACVR1C_TGFbR2", "ACVR1B_TGFbR2", "ACVR1C_TGFb…
$ agonist "TGFb agonist", "TGFb agonist", "TGFb agonist", "TGFb agonist", "TGFb agonist", "TGFb agonist", "TGFb …
$ antagonist "TGFb antagonist", "TGFb antagonist", "TGFb antagonist", "TGFb antagonist", "TGFb antagonist", "TGFb a…
$ co_A_receptor "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", "", ""…
$ co_I_receptor "TGFb inhibition receptor", "TGFb inhibition receptor", "TGFb inhibition receptor", "TGFb inhibition r…
$ evidence "KEGG: mmu04350", "KEGG: mmu04350", "KEGG: mmu04350", "PMID: 27449815", "PMID: 27449815", "PMID: 27449…
$ annotation "Secreted Signaling", "Secreted Signaling", "Secreted Signaling", "Secreted Signaling", "Secreted Sign…
$ interaction_name_2 "Tgfb1 - (Tgfbr1+Tgfbr2)", "Tgfb2 - (Tgfbr1+Tgfbr2)", "Tgfb3 - (Tgfbr1+Tgfbr2)", "Tgfb1 - (Acvr1b+T…

Usar un subconjunto de la base de datos para el análisis

CellChatDB.use <- subsetDB(CellChatDB, search = "Secreted Signaling")
cellchat@DB <- CellChatDB.use
cellchat <- subsetData(cellchat)
future::plan("multisession", workers = 4) # Paralelizar la computación
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)

Inferencia y análisis de comunicación célula a célula

cellchat <- computeCommunProb(cellchat)
triMean is used for calculating the average gene expression per cell group.
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
Warning in log(x) : NaNs produced
[1] ">>> Run CellChat on sc/snRNA-seq data <<< [2024-06-20 12:48:31.880874]"
|====================== | 18%Error in if (sum(P1_Pspatial) == 0) { :
missing value where TRUE/FALSE needed

Thanks,

[sqjin]

@estebanelias You should use the normalized data instead of count data. Please check the tutorial:

access count data matrix

counts <- t(as.matrix(ad$X))

normalize the count data if the normalized data is not available in the .h5ad file

library.size <- Matrix::colSums(counts)
data.input <- as(log1p(Matrix::t(Matrix::t(counts)/library.size) * 10000), "dgCMatrix")