L–R pairs mediating the cell-type pairs interaction (Fig. 5g)

[SoniaHermoso]

Hello!

First, thank you by advanced for the package developed for users.
Reading your article in detail ("Estimation of cell lineages in tumors from spatial transcriptomics data"), I wonder how is it possible to plot the figure 5.g: L–R pairs mediating the CAF–M2 interaction in the current breast cancer tissue. The direction of an arrow source from ligand to receptor.
If you can share the code to know which proteins are interacting in our data.

thanks
Sonia H.

[beibeiru]

Hi, @SoniaHermoso ,

The results in Fig. 5g were calculated using both the current breast cancer Visium data and multiple scRNA-seq datasets.

Here are the steps:
In step 1, we can get the co-expressed L-R pairs in CAF–M2 colocalized spots. However, we can not determine the cell-type source of ligands and receptors. The scRNA-seq datasets would be integrated to address this issue.

In step 2, we examined which L–R pairs have significant co-expression between CAF and M2 cells across our collected scRNA-seq datasets.

In step 3, we selected the overlapped L–R pairs between steps 1 and 2, which are shown in Fig. 5g.

(1) We computed the Spearman correlation of ligand and receptor pairs across CAF–M2 colocalized spots.

# extract CAF–M2 colocalized spots
groupMat <- SpaCET_obj@results$CCI$interaction$groupMat
CAFM2_colocalized_spots <- names(groupMat)[unlist(groupMat[1,])%in%"Both"]

# extract expression matrix
st.matrix.data <- as.matrix(SpaCET_obj@input$counts)

# remove low expressed genes
top5p <- round(ncol(st.matrix.data)*0.05)
st.matrix.data <- st.matrix.data[rowSums(st.matrix.data>0)>top5p,] 

st.matrix.data <- t(t(st.matrix.data)*1e5/colSums(st.matrix.data))

# extract expression matrix in the CAF–M2 colocalized spots
expr <- st.matrix.data[, CAFM2_colocalized_spots]

# load ligand-receptor pairs
dataPath  <- file.path(system.file(package = "SpaCET"), "extdata")
L_R <- read.csv(paste0(dataPath,"/Ramilowski2015.txt"),as.is=T,sep="\t",row.names=1)
rownames(L_R) <- paste0(L_R[,"Ligand.ApprovedSymbol"],"->",L_R[,"Receptor.ApprovedSymbol"])

LR <- c(L_R[,"Ligand.ApprovedSymbol"], L_R[,"Receptor.ApprovedSymbol"])

# filter the common items
olp <- intersect(LR, rownames(expr))
L_R <- L_R[L_R[,"Ligand.ApprovedSymbol"]%in%olp&L_R[,"Receptor.ApprovedSymbol"]%in%olp,]
expr <- expr[olp,]

# deconvolution results
res_deconv <- SpaCET_obj@results$deconvolution$propMat

summary_df <- data.frame()

for(i in 1:nrow(L_R) )
{
	library(ppcor)
	ppcor_res <- pcor.test(
		expr[L_R[i,"Ligand.ApprovedSymbol"],],
		expr[L_R[i,"Receptor.ApprovedSymbol"],],
		res_deconv["Malignant",CAFM2_colocalized_spots],
		method="spearman")
	# Here, we use partial correlation by using malignant fraction as a covariate
	
	summary_df[i,"L_R"] <- rownames(L_R)[i]
	summary_df[i,"rho"] <- ppcor_res[1,"estimate"]
	summary_df[i,"p"] <- ppcor_res[1,"p.value"]
	summary_df[i,"padj"] <-0
}

summary_df[,"padj"] <- p.adjust(summary_df[,"p"], method = "BH")
summary_df <- summary_df[order(summary_df[,"rho"],decreasing=T),]

# show all significant co-expressed Ligand-Receptor pairs
stLR <- summary_df[!is.na(summary_df[,"padj"])&summary_df[,"rho"]>0.15, 1 ]

(2) We examined which L–R pairs have significant co-expression between CAF and M2 cells across our collected scRNA-seq datasets by using a similar strategy from CellPhoneDB. These scRNA-seq can be downloaded at https://doi.org/10.5281/zenodo.7466024.

# 5 out of 10 scRNA-seq datasets we collected have both CAF and M2 cell types
scs <- c(
	"E-MTAB-6149_NSCLC",
	"GSE103322_HNSCC",
	"GSE115978_Melanoma",
	"GSE127465_NSCLC",
	"GSE146771_CRC_SS2"
)

for(sc in scs)
{
	load(paste0(sc,"_TPM.rda"))
	load(paste0(sc,"_meta.rda"))
	
	sc.matrix.data <- round(log2(sc.matrix.data/10+1),1)
	
	sc.CAF <- rownames(sc.matrix.anno)[sc.matrix.anno[,"uniformCellType"]%in%"CAF"]
	sc.M2 <- rownames(sc.matrix.anno)[sc.matrix.anno[,"uniformCellTypeSub"]%in%"M2"]
	
	sc.matrix.data.sub <- sc.matrix.data[,c(sc.CAF,sc.M2)]
	sc.matrix.anno.sub <- sc.matrix.anno[c(sc.CAF,sc.M2),]
	
	counts <- sc.matrix.data.sub
	meta <- rbind(cbind(sc.CAF,"CAF"),cbind(sc.M2,"M2"))
	

	# load ligand-receptor pairs
	dataPath  <- file.path(system.file(package = "SpaCET"), "extdata")
	L_R <- read.csv(paste0(dataPath,"/Ramilowski2015.txt"),as.is=T,sep="\t",row.names=1)
	rownames(L_R) <- paste0(L_R[,"Ligand.ApprovedSymbol"],"->",L_R[,"Receptor.ApprovedSymbol"])

	LR2015 <- data.frame(
		L=L_R[,"Ligand.ApprovedSymbol"], 
		R= L_R[,"Receptor.ApprovedSymbol"], 
		stringsAsFactors=FALSE
	)
	
	Ls <- unique(LR2015[,"L"])
	Rs <- unique(LR2015[,"R"])
	
	# intersect(Ls,Rs)
	# [1] "DLK2"
	
	LR2015 <- LR2015[LR2015[,"R"]!="DLK2",]
		
	
	L_R <- LR2015
	L_R <- L_R[!duplicated(paste0(L_R[,1],"->",L_R[,2])),]
	rownames(L_R) <- paste0(L_R[,1],"->",L_R[,2])
	LR <- c(L_R[,1],L_R[,2])
	
	
	olp <- intersect(LR,rownames(counts))
	L_R <- L_R[L_R[,1]%in%olp&L_R[,2]%in%olp,]
	counts <- counts[olp,]

	cellTypes <- unique(meta[,2])
	cellTypeMat <- matrix(0,ncol=length(cellTypes),nrow=length(cellTypes))
	colnames(cellTypeMat) <- cellTypes
	rownames(cellTypeMat) <- cellTypes
	cellTypeMat.m <- reshape::melt(cellTypeMat)
	cellTypePairs <- paste0(cellTypeMat.m[,1],"->",cellTypeMat.m[,2])


	collaps <- function(x){
		counts_sub <- counts[,colnames(counts)%in%meta_s[meta_s[,2]==x,1]]
		rowMeans(counts_sub)
	}
	collaps2 <- function(x){
		counts_sub <- counts[,colnames(counts)%in%meta_s[meta_s[,2]==x,1]]
		temp <- rowMeans(counts_sub)
		expr_flag <- rowSums(counts_sub>0.5)
		temp[expr_flag<dim(counts_sub)[2]*0.1] <- 0
		temp
	}
	
	meta_s <- meta
	LR_cellType <- sapply(cellTypes, collaps2)
	
	
	inter <- function(x){
		temp <- t(LR_cellType[x[1],,drop=F])%*%LR_cellType[x[2],,drop=F]
		temp.m <- reshape::melt(temp)
		rownames(temp.m) <- paste0(temp.m[,1],"->",temp.m[,2])
		temp.m[,3]
	}
	L_R_cellType_cellType <- t(apply(L_R,1,inter))
	colnames(L_R_cellType_cellType) <- cellTypePairs
	
	L_R_cellType_cellType_p <- L_R_cellType_cellType
	L_R_cellType_cellType_p[L_R_cellType_cellType_p>-1] <- 0
	
	for(i in 1:1000)
	{
		meta_s <- meta
		meta_s[,2] <- sample(meta_s[,2])
		LR_cellType <- sapply(cellTypes, collaps2)
	
		L_R_cellType_cellType_s <- t(apply(L_R,1,inter))
		colnames(L_R_cellType_cellType_s) <- cellTypePairs
		
		L_R_cellType_cellType_p <- L_R_cellType_cellType_p+(L_R_cellType_cellType_s>=L_R_cellType_cellType)
	}
	
	L_R_cellType_cellType_p <- (L_R_cellType_cellType_p+1)/1001


	score <- as.matrix(L_R_cellType_cellType)
	pv <- as.matrix(L_R_cellType_cellType_p)
	
	c1 <- sapply(strsplit(colnames(pv),"->",fixed=T),function(x) return(paste0(x[1])))
	c2 <- sapply(strsplit(colnames(pv),"->",fixed=T),function(x) return(paste0(x[2])))

	score <- score[,c1!=c2]
	pv <- pv[,c1!=c2]
	
	score <- log2(score+1)
	pv <- -log10(pv)
	
	score.m <- reshape::melt(score)
	pv.m <- reshape::melt(pv)

	comb <- cbind(score.m,pv.m[,3])
	colnames(comb) <- c("L_R","Cell_Cell","Expr_Product","log10pv")
	
	
	write.csv(comb,paste0("sc_",sc,"_LR.csv"),quote=F)

} # sc 


summary_table <- data.frame()
for(sc in scs)
{
	comb <- read.csv(paste0("sc_",sc,"_LR.csv"),row.names=1,header=T,as.is=T)

	comb <- comb[comb[,4]>3,]
	
	rownames(comb) <- paste0(comb[,1],"_",comb[,2])
	
	summary_table[rownames(comb),"L_R"] <- comb[rownames(comb),"L_R"]
	summary_table[rownames(comb),"Cell_Cell"] <- comb[rownames(comb),"Cell_Cell"]
	summary_table[rownames(comb),sc] <- comb[rownames(comb),4]
}

# Compute how many scRNA-seq datasets each L-R pair is significant co-expression between CAF and M2 cells in
rowSums(!is.na(as.matrix(summary_table[,scs])))-> stat

summary_table <- cbind(summary_table,stat)

# we select the L-R pairs that are significant co-expression between CAF and M2 cells in at least two scRNA-seq datesets
scLR <- unique(summary_table[summary_table[,"stat"]>=2,1])

(3) We get the overlapped L-R pairs from step 1 and 2.

olpLR <- intersect(stLR, scLR)

# show cell-type sources of ligands and receptors
summary_table[summary_table[,1]%in%olpLR, ]

Best,
Beibei