目的

通過(guò)celltalker找出小鼠的單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)的受體配體互作信息蛛枚，并繪制類似這樣的圈圖

具體用法

看如下網(wǎng)站：

• http://www.reibang.com/p/177bbf7e9d67
• https://arc85.github.io/celltalker/articles/celltalker.html#celltalker

問題解決

與GSVA或者CellPhoneDB的情況相同

解決方案：像我之前構(gòu)建GSVA的小鼠gmt文件的方法相同http://www.reibang.com/p/2d659a4ec167伏穆，將 ramilowski_pairs 表示的人的受體配體的關(guān)系文件厢拭，轉(zhuǎn)換成小鼠的同源基因表示形式

實(shí)際操作

數(shù)據(jù)導(dǎo)入

#Robin 2020/03/04
#For celltalker analysis mouse scRNA-seq Data
library(Seurat)
library(ggplot2)
library(celltalker)
library(dplyr)

setwd('/share/nas1/Data/Users/luohb/20200303_Ligand_interaction')
#Load data
mouse<-readRDS("../20200221/scaled_merge_seurat_rna.rds")
ser.obj<-mouse

用Python轉(zhuǎn)換成小鼠的同源基因的受體配體關(guān)系表

在R中：

write.table(ramilowski_pairs, file='./ramilowski_pairs.xls', sep='\t', quote=F)

在Python中：

#!/usr/bin/python
#Author:Robin 20200220
#For transfer human to mouse

human_set = set()
h2m = dict()
m2m = dict()

# human gene ID to mouse Entrez ID
with open('./mart_export.txt') as f1:     
    for line in f1:
        line = line.strip('\n')
        lst = line.split('\t')    

        if lst[3]:
            if lst[3] in human_set:
                h2m[lst[3]].append(lst[0])
            else:
                h2m[lst[3]] = [ lst[0] ]
                human_set.add(lst[3])        


# mouse Entrez ID to mouse symbol ID
with open('./features.tsv') as f2:
    for line in f2:
        line = line.strip('\n')
        lst = line.split('\t')
        m2m[lst[0]] = lst[1]

# 遍歷ramilowski_pairs文件每行足绅，進(jìn)行替換
with open('./ramilowski_pairs.xls') as f3:
    with open('./mm.ramilowski_pairs.xls', 'w') as f4:
        for line in f3:
            line = line.strip('\n')
            lst = line.split('\t')
            
            ligand = lst[0]
            receptor = lst[1]
            mouse_pair = lst[2]
           
            mouse_entrez_ligand_lst = h2m.get(ligand, "")         
            mouse_entrez_receptor_lst = h2m.get(receptor, "") 
            
            for mouse_entrez_ligand in mouse_entrez_ligand_lst:
                for mouse_entrez_receptor in mouse_entrez_receptor_lst:
                    mouse_symbol_ligand = m2m[mouse_entrez_ligand]
                    mouse_symbol_receptor = m2m[mouse_entrez_receptor]
                    mouse_pair = mouse_symbol_ligand + "_" + mouse_symbol_receptor
                    
                    line = [mouse_symbol_ligand, mouse_symbol_receptor, mouse_pair]
                    line = '\t'.join(line)
                    line = line + '\n'
                    print(line)
                    f4.write(line)

回到R中：

# 讀入已經(jīng)構(gòu)建好的ramilowski_pairs文件
ramilowski_pairs<-read.table('./mm.ramilowski_pairs.xls', header=F, sep='\t', col.names=c('ligand','receptor','pair'))

head(ramilowski_pairs)
dim(ramilowski_pairs) #該數(shù)據(jù)集中有2,557個(gè)特異的配體/受體對(duì)
View(ramilowski_pairs)


#獲取小鼠與人的同源基因的對(duì)應(yīng)關(guān)系，然后用人的同源基因替換掉原本小鼠對(duì)應(yīng)的基因
mart_export<-read.csv('./mart_export.txt', sep='\t', header=T)
barcode_feature<-read.csv('./features.tsv', sep='\t', header=F, col.names=c('mus_entrezID', 'mus_symbol', "type"))

expr.mat <- GetAssayData(ser.obj,slot="counts")

#在我們的數(shù)據(jù)集中識(shí)別配體和受體
ligs <- as.character(unique(ramilowski_pairs$ligand))
recs <- as.character(unique(ramilowski_pairs$receptor))

gene_list<-rownames(ser.obj)
ligs.present <- gene_list[gene_list %in% ligs]
recs.present <- gene_list[gene_list %in% recs]
genes.to.use <- union(ligs.present, recs.present) # 取并集


# 使用之前已經(jīng)使用FindAllMarkers跑好了的所有cluster差異表達(dá)基因淳玩，用于區(qū)分組之間差異表達(dá)的配體和受體
Idents(ser.obj) <- ser.obj@meta.data$seurat_clusters
markers <- read.table('/share/nas1/Data/Users/yinl/Project/personality/20191227/markgene/all_sig_markgene.xls',
                                   header=T)
marker_gene<-markers$gene

ligs.recs.use <- unique(marker_gene)
length(ligs.recs.use) #產(chǎn)生7648個(gè)配體和受體以進(jìn)行評(píng)估


#過(guò)濾ramilowski配受對(duì)
interactions.forward1 <- ramilowski_pairs[as.character(ramilowski_pairs$ligand) %in% ligs.recs.use,]
interactions.forward2 <- ramilowski_pairs[as.character(ramilowski_pairs$receptor) %in% ligs.recs.use,]
interact.for <- rbind(interactions.forward1, interactions.forward2)
dim(interact.for)

構(gòu)建celltalker的輸入文件：

#生成celltalker的輸入數(shù)據(jù)
#Create data for celltalker
row_expr.mat <- GetAssayData(ser.obj,slot="counts")
expr.mat<-row_expr.mat
rownames(expr.mat)<-rownames(expr.mat)
rownames(expr.mat)<-gsub('\\.[0-9]*', '', rownames(expr.mat))  #去除版本號(hào)

defined.clusters <- ser.obj@meta.data$seurat_clusters
names(defined.clusters)<-colnames(ser.obj)

defined.groups <- ser.obj@meta.data$group
names(defined.groups)<-colnames(ser.obj)

defined.replicates <- ser.obj@meta.data$orig.ident
names(defined.replicates)<-colnames(ser.obj)

reshaped.matrices <- reshape_matrices(count.matrix=expr.mat,
                                      clusters=defined.clusters,
                                      groups=defined.groups,
                                      replicates=defined.replicates,
                                      ligands.and.receptors=interact.for)

unnest(reshaped.matrices,cols="samples")
names(pull(unnest(reshaped.matrices,cols="samples"))[[1]])
#在此初始步驟中旁瘫，我們要做的是將我們的整體表達(dá)矩陣中給每個(gè)樣本中分離出單獨(dú)的表達(dá)矩陣祖凫。
#接下來(lái)，使用create_lig_rec_tib函數(shù)為每個(gè)組創(chuàng)建一組一致表達(dá)的配體和受體酬凳。
# cells.reqd=10：每個(gè)cluster中至少有10個(gè)細(xì)胞表達(dá)了配體/受體
# freq.pos.reqd=0.5：至少有50％重復(fù)個(gè)體中表達(dá)的配體/受體
consistent.lig.recs <- create_lig_rec_tib(exp.tib=reshaped.matrices,
                                          clusters=defined.clusters,
                                          groups=defined.groups,
                                          replicates=defined.replicates,
                                          cells.reqd=300,
                                          freq.pos.reqd=0.5,
                                          ligands.and.receptors=interact.for)


unnest(consistent.lig.recs, cols="lig.rec.exp")
pull(unnest(consistent.lig.recs[1,2],cols="lig.rec.exp")[1,2])[[1]]
#從每個(gè)實(shí)驗(yàn)組的每個(gè)cluster中獲得了一致表達(dá)的配體和受體的列表惠况。

生成Celltalker對(duì)象，并繪制circos圈圖

#Determine putative ligand/receptor pairs
# freq.group.in.cluster: 只對(duì)包含細(xì)胞數(shù)大于總細(xì)胞數(shù)5%的簇進(jìn)行互作分析
put.int<-putative_interactions(ligand.receptor.tibble=consistent.lig.recs,
                               clusters=defined.clusters,
                               groups=defined.groups,
                               freq.group.in.cluster=0.05,
                               ligands.and.receptors=interact.for)

# Identifying and visualizing unique ligand/receptor pairs in a group
###Identify unique ligand/receptor interactions present in each sample
unique.ints <- unique_interactions(put.int,group1="Fat",group2="Health",interact.for)

#Get data to plot circos for Fat
fat.to.plot <- pull(unique.ints[1,2])[[1]]
for.circos.fat <- pull(put.int[1,2])[[1]][fat.to.plot]

par(mar=c(1, 1, 1, 1))
p_fat<-circos_plot(interactions=for.circos.fat, clusters=defined.clusters)

#Get data to plot circos for Health
health.to.plot <- pull(unique.ints[2,2])[[1]]
for.circos.health <- pull(put.int[2,2])[[1]][health.to.plot]

p2<-circos_plot(interactions=for.circos.health, clusters=defined.clusters)

#all
to.plot <- pull(unique.ints[2])[[1]]
for.circos <- pull(put.int[2])[[1]][to.plot]

p3<-circos_plot(interactions=for.circos, clusters=defined.clusters)