細胞通訊分析流程與方法

細胞通訊分析流程與方法：細胞通訊分析的基礎是基因表達數(shù)據(jù)和配體-受體數(shù)據(jù)庫信息，例如轉(zhuǎn)錄組數(shù)據(jù)表明A、B細胞分別表達了基因α和β，通過數(shù)據(jù)庫查詢α和β是配體-受體關系，則認為A-B通過α-β途徑進行了通訊。

輸入數(shù)據(jù)格式要求：

（1）細胞注釋文件（meta.txt，兩列，分別為細胞id和細胞類型）

（2）基因表達矩陣（counts.txt）

#下載軟件（按順序安裝）

conda?create?-n?cellphonedb2?python=3.7

conda?install?Cython

conda?install?h5py

conda?install?scikit-learn

conda?install?r

conda?install?rpy2

pip?install?https://pypi.tuna.tsinghua.edu.cn/simple?cellphonedb

pip?install?-i?https://pypi.tuna.tsinghua.edu.cn/simple?markupsafe==2.0.1

#R包

install.packages(‘ggplot2’)

install.packages(‘heatmap’)

#data1—> R

#使用RCTD注釋的空間數(shù)據(jù)進行后續(xù)分析

results?<-?myRCTD@results

results_df?<-?results$results_df

barcodes?=?rownames(results_df[results_df$spot_class?!=?“reject”?&?puck@nUMI?>=?1,])

my_table?=?puck@coords[barcodes,]

my_table$class?=?results_df[barcodes,]$first_type

meta?=?my_table?%>%?select(class) %>%?rownames_to_column(var?=?‘cell’) %>%?dplyr::rename(cell_type?=?class)

write.table(meta?,file?=?‘/xxx/Test/meta.txt’, sep?=?‘\t’, quote?=?F, row.names?=?T)

#data2—> R

expr?<-?Read10X(‘/xxx/BSTViewer_project/subdata/L13_heAuto/’,cell.column?=?1)

object?<-?CreateSeuratObject(counts?=?expr,assay?=?“Spatial”)

counts?=?object[[‘Spatial’]]@counts?%>%?as.data.frame() %>%?rownames_to_column(var?=?‘Gene’)

write.table(counts?,file?=?‘/xxx/Test/counts.txt’, sep?=?‘\t’, quote?=?F, row.names?=?T)

#鼠源基因轉(zhuǎn)化

#CellPhoneDB 只能識別人類基因，如果是鼠源基因需要進行同源基因的轉(zhuǎn)化

#下載人鼠同源基因?qū)P系

http://asia.ensembl.org/biomart/martview/b9f8cc0248e4714ba8e0484f0cbe4f02

#細胞通訊分析流程與方法參考步驟

https://www.jianshu.com/p/922a7f2c21fc

#下載的對應關系表應該有四列

Gene_stable_ID??|???Gene_name???|???Mouse_gene_stable_ID????|???Mouse_gene_name

#將對應關系(human_mouse_gene.txt)替換到 data2 矩陣—>Shell

awk?-F?‘\t’?‘BEGIN{OFS == ‘\t’}ARGIND==1{a[$4]=$1}ARGIND==2{if(FNR ==1){print}else{$1 = a[$1];print}}’?human_mouse_gene.txt?counts.txt?|?sed?‘/^ /d’?|?sed?‘1s/^/Gene\t/’?|?sed?‘s/\s/\t/g’?>?counts_new.txt

#從遠程存儲庫中列出cellphonedb可用數(shù)據(jù)庫版本:

cellphonedb?database?list_remote

#查看本地數(shù)據(jù)庫

cellphonedb?database?list_local

#下載遠程數(shù)據(jù)庫

cellphonedb?database?download?–version?<version_spec|latest>

#運行cellphonedb —> Shell

cellphonedb?method?statistical_analysis?–output-path?test_output?meta.txt?counts_new.txt

#輸出文件

deconvoluted.txt????|???pvalues.txt?|???significant_means.txt

test.count_network.txt??|???means.txt???|???test.interaction_count.txt

#氣泡圖—> Shell

cellphonedb?plot?dot_plot?–pvalues-path?test_output/pvalues.txt?–means-path?test_output/means.txt?–output-path?test_output??–output-name?test.dotplot.pdf

#僅繪制自己關注的細胞和互作—> Shell

cellphonedb?plot?dot_plot?–pvalues-path?test_output/pvalues.txt?–means-path?test_output/means.txt?–output-path?test_output??–output-name?test.dotplot2.pdf?–rows?test_output/row.txt?–columns?test_output/col.txt

注：橫坐標是細胞類型互作，縱坐標是蛋白互作，點越大表示 p 值越小，顏色代表平均表達量

#熱圖—> Shell

cellphonedb?plot?heatmap_plot?–pvalues-path?test_output/pvalues.txt?–output-path?test_output?–pvalue?0.05?–count-name?test.heatmap_count.pdf?–log-name?test.heatmap_log_count.pdf?–count-network-name?test.count_network.txt?–interaction-count-name?test.interaction_count.txt?meta.txt

注：橫縱坐標表示細胞類型，顏色深藍色到紫紅色代表互作數(shù)由低到高

#細胞通訊分析流程與方法：細胞類型網(wǎng)絡圖—> R

#這里借用另一款細胞通訊分析軟件——CellChat(https://github.com/sqjin/CellChat)的繪圖方法

#加載R包

library(igraph)

library(tidyr)

library(stats)

library(reshape2)

library(Matrix)

library(grDevices)

library(ggplot2)

#使用CellChat作者提供的繪制互作網(wǎng)絡圖的function

scPalette?<-?function(n) {

colorSpace?<-?c(‘#E41A1C’,’#377EB8′,’#4DAF4A’,’#984EA3′,’#F29403′,’#F781BF’,’#BC9DCC’,’#A65628′,’#54B0E4′,’#222F75′,’#1B9E77′,’#B2DF8A’,

‘#E3BE00′,’#FB9A99′,’#E7298A’,’#910241′,’#00CDD1′,’#A6CEE3′,’#CE1261′,’#5E4FA2′,’#8CA77B’,’#00441B’,’#DEDC00′,’#B3DE69′,’#8DD3C7′,’#999999′)

if?(n?<=?length(colorSpace)) {

colors?<-?colorSpace[1:n]

} else?{

colors?<-?grDevices::colorRampPalette(colorSpace)(n)

}

return(colors)

}

netVisual_circle?<-function(net, color.use?=?NULL,title.name?=?NULL, sources.use?=?NULL, targets.use?=?NULL, idents.use?=?NULL, remove.isolate?=?FALSE, top?=?1,

weight.scale?=?FALSE, vertex.weight?=?20, vertex.weight.max?=?NULL, vertex.size.max?=?NULL, vertex.label.cex=1,vertex.label.color=?“black”,

edge.weight.max?=?NULL, edge.width.max=8, alpha.edge?=?0.6, label.edge?=?FALSE,edge.label.color=’black’,edge.label.cex=0.8,

edge.curved=0.2,shape=’circle’,layout=in_circle(), margin=0.2, vertex.size?=?NULL,

arrow.width=1,arrow.size?=?0.2){

if?(!is.null(vertex.size)) {

warning(“‘vertex.size’ is deprecated. Use `vertex.weight`”)

}

if?(is.null(vertex.size.max)) {

if?(length(unique(vertex.weight)) ==?1) {

vertex.size.max?<-?5

} else?{

vertex.size.max?<-?15

}

}

options(warn?=?-1)

thresh?<-?stats::quantile(net, probs?=?1-top)

net[net?<?thresh] <-?0

if?((!is.null(sources.use)) |?(!is.null(targets.use)) |?(!is.null(idents.use)) ) {

if?(is.null(rownames(net))) {

stop(“The input weighted matrix should have rownames!”)

}

cells.level?<-?rownames(net)

df.net?<-?reshape2::melt(net, value.name?=?“value”)

colnames(df.net)[1:2] <-?c(“source”,”target”)

# keep the interactions associated with sources and targets of interest

if?(!is.null(sources.use)){

if?(is.numeric(sources.use)) {

sources.use?<-?cells.level[sources.use]

}

df.net?<-?subset(df.net, source?%in%?sources.use)

}

if?(!is.null(targets.use)){

if?(is.numeric(targets.use)) {

targets.use?<-?cells.level[targets.use]

}

df.net?<-?subset(df.net, target?%in%?targets.use)

}

if?(!is.null(idents.use)) {

if?(is.numeric(idents.use)) {

idents.use?<-?cells.level[idents.use]

}

df.net?<-?filter(df.net, (source?%in%?idents.use) |?(target?%in%?idents.use))

}

df.net$source?<-?factor(df.net$source, levels?=?cells.level)

df.net$target?<-?factor(df.net$target, levels?=?cells.level)

df.net$value[is.na(df.net$value)] <-?0

net?<-?tapply(df.net[[“value”]], list(df.net[[“source”]], df.net[[“target”]]), sum)

}

net[is.na(net)] <-?0

if?(remove.isolate) {

idx1?<-?which(Matrix::rowSums(net) ==?0)

idx2?<-?which(Matrix::colSums(net) ==?0)

idx?<-?intersect(idx1, idx2)

net?<-?net[-idx, ]

net?<-?net[, -idx]

}

g?<-?graph_from_adjacency_matrix(net, mode?=?“directed”, weighted?=?T)

edge.start?<-?igraph::ends(g, es=igraph::E(g), names=FALSE)

coords<-layout_(g,layout)

if(nrow(coords)!=1){

coords_scale=scale(coords)

}else{

coords_scale<-coords

}

if?(is.null(color.use)) {

color.use?=?scPalette(length(igraph::V(g)))

}

if?(is.null(vertex.weight.max)) {

vertex.weight.max?<-?max(vertex.weight)

}

vertex.weight?<-?vertex.weight/vertex.weight.max*vertex.size.max+5

loop.angle<-ifelse(coords_scale[igraph::V(g),1]>0,-atan(coords_scale[igraph::V(g),2]/coords_scale[igraph::V(g),1]),pi-atan(coords_scale[igraph::V(g),2]/coords_scale[igraph::V(g),1]))

igraph::V(g)$size<-vertex.weight

igraph::V(g)$color<-color.use[igraph::V(g)]

igraph::V(g)$frame.color?<-?color.use[igraph::V(g)]

igraph::V(g)$label.color?<-?vertex.label.color

igraph::V(g)$label.cex<-vertex.label.cex

if(label.edge){

igraph::E(g)$label<-igraph::E(g)$weight

igraph::E(g)$label?<-?round(igraph::E(g)$label, digits?=?1)

}

if?(is.null(edge.weight.max)) {

edge.weight.max?<-?max(igraph::E(g)$weight)

}

if?(weight.scale?==?TRUE) {

#E(g)$width<-0.3+edge.width.max/(max(E(g)$weight)-min(E(g)$weight))*(E(g)$weight-min(E(g)$weight))

igraph::E(g)$width<-?0.3+igraph::E(g)$weight/edge.weight.max*edge.width.max

}else{

igraph::E(g)$width<-0.3+edge.width.max*igraph::E(g)$weight

}

igraph::E(g)$arrow.width<-arrow.width

igraph::E(g)$arrow.size<-arrow.size

igraph::E(g)$label.color<-edge.label.color

igraph::E(g)$label.cex<-edge.label.cex

igraph::E(g)$color<-?grDevices::adjustcolor(igraph::V(g)$color[edge.start[,1]],alpha.edge)

if(sum(edge.start[,2]==edge.start[,1])!=0){

igraph::E(g)$loop.angle[which(edge.start[,2]==edge.start[,1])]<-loop.angle[edge.start[which(edge.start[,2]==edge.start[,1]),1]]

}

radian.rescale?<-?function(x, start=0, direction=1) {

c.rotate?<-?function(x) (x?+?start) %%?(2?*?pi) *?direction

c.rotate(scales::rescale(x, c(0, 2?*?pi), range(x)))

}

label.locs?<-?radian.rescale(x=1:length(igraph::V(g)), direction=-1, start=0)

label.dist?<-?vertex.weight/max(vertex.weight)+2

plot(g,edge.curved=edge.curved,vertex.shape=shape,layout=coords_scale,margin=margin, vertex.label.dist=label.dist,

vertex.label.degree=label.locs, vertex.label.family=”Helvetica”, edge.label.family=”Helvetica”) # “sans”

if?(!is.null(title.name)) {

text(0,1.5,title.name, cex?=?1.1)

}

# https://www.andrewheiss.com/blog/2016/12/08/save-base-graphics-as-pseudo-objects-in-r/

# grid.echo()

# gg <- ?grid.grab()

gg?<-?recordPlot()

return(gg)

}

#整體互作網(wǎng)絡圖繪制

#讀取cellphonedb的分析數(shù)據(jù)

count_net?<-?read.delim(“/xxx/test_output/test.count_network.txt”, check.names?=?FALSE)

#數(shù)據(jù)格式處理

count_inter?<-?count_net

count_inter$count?<-?count_inter$count/100

count_inter?<-?spread(count_inter, TARGET, count)

rownames(count_inter) <-?count_inter$SOURCE

count_inter?<-?count_inter[, -1]

count_inter?<-?as.matrix(count_inter)

#繪圖

netVisual_circle(count_inter,weight.scale?=?T)

#細胞通訊分析流程與方法：每種細胞與其他細胞的互作

par(mfrow?=?c(1,3), xpd=TRUE)

for?(i?in?1:nrow(count_inter)) {

mat2?<-?matrix(0, nrow?=?nrow(count_inter), ncol?=?ncol(count_inter), dimnames?=?dimnames(count_inter))

mat2[i, ] <-?count_inter[i, ]

netVisual_circle(mat2,

weight.scale?=?T,

edge.weight.max?=?max(count_inter),

title.name?=?rownames(count_inter)[i],

arrow.size=0.2)

}