哇！单细胞测序-配体受体互作分析原来可以这么简单又高大上

NGS系列文章包括NGS基础、转录组分析（Nature重磅综述|关于RNA-seq你想知道的全在这）、ChIP-seq分析（ChIP-seq基本分析流程）、单细胞测序分析(重磅综述：三万字长文读懂单细胞RNA测序分析的最佳实践教程（原理、代码和评述）)、DNA甲基化分析、重测序分析、GEO数据挖掘（典型医学设计实验GEO数据分析(step-by-step)-Limma差异分析、火山图、功能富集）等内容。

单细胞转录组测序发展至今，我们发现许多文章的最后一部分都会落到配受体结合，可是如何挑配受体，哪些基因可能是配受体，我一脸懵逼。。。

于是，我一不小心发现了celltalker()，大家可以尝试一下哦，嘻嘻，废话不多说。

对单细胞RNAseq数据可能进行的多种分析之一是评估细胞间的交流(cell-cellcommunication)。celltalker通过寻找细胞群内和细胞群之间已知的配体和受体对的表达来评估细胞之间的交流。我们采用的配受体数据库来自Ramilowski等人于2015年在Naturecommunication上发表的Adraftnetworkofligand-receptor-mediatedmulticellularsignallinginhuman描述的一组配体和受体。我们建议使用此数据集作为起点，并整理自己的已知配体和受体列表。另外Tormo2018年发表的Nature文章Single-cellreconstructionoftheearlymaternal-fetalinterfaceinhumans扩展了受体和配体对，也会应用于cellTalker的更新版中。

为了获得可靠的结果，我们要求每个组中都有多个重复样品，并且只对不同组间一致性表达的配体和受体感兴趣（而仅在单个重复中发现的互作可信度低）。我们通过评估每组中各个重复的表达矩阵并仅对满足一定阈值（这个阈值随意性也比较强）的相互作用进行提取。

配体和受体相互作用的差异至少在三种方面具有生物学意义：

在一组细胞中独特地存在；

各个cluster间配体或受体的互作差异；

参与组间配体和受体相互作用的细胞网络的差异。

我们提供了评估每种潜在生物学差异的方法，并提供了具体示例。

在这个vignette中，我们展示了cellTalker在评估健康捐献者外周血（N=2）和扁桃体（N=3）中鉴定配体/受体相互作用的基本应用。该数据可从我们最近发布的数据集GSE139324中获得(Cilloetal,Immunity2020)。

展示Celltalker应用于10XGenomics数据的的标准用法。具体分为下面几步：

使用标准的Seurat工作流程（）对数据进行聚类;

使用Celltalker建立样品组中稳定表达的配体和受体的列表;

确定配体/受体相互作用;

评估组之间特异表达的配体/受体对;

识别和可视化组特异的配体/受体对;

library(devtools)install_github("arc85/celltalker")library(celltalker)

使用Seurat进行标准的聚类分析和免疫谱系识别（假设已从GEO下载了rawmatrix）。（重磅综述：三万字长文读懂单细胞RNA测序分析的最佳实践教程（原理、代码和评述））

suppressMessages({library(Seurat)library(celltalker)})读取rawdataGRCh38根据需要调整为其它基因组版本(,"data_matrices",,"GRCh38",sep="/")setwd()()标准Seurat工作流程()()()

PCA分析

()

获得对各个主成分贡献比较大的基因(用了这么多年的PCA可视化竟然是错的！！！)

PC_1Positive:S100A6,IL32,S100A4,ANXA1,VIM,FTL,TRBC1,SRGN,S100A9,S100A8TYROBP,LYZ,CTSW,XIST,NEAT1,VCAN,S100A12,FCER1G,S100A11,FCN1PLAC8,ID2,CCL5,NKG7,CST3,CSTA,ZFP36,IL1B,MT2A,KLRB1Negative:RGS13,KIAA0101,NUSAP1,AURKB,MKI67,BIRC5,TYMS,TOP2A,TK1,CDKN3UBE2C,PTTG1,CDK1,STMN1,CCNB2,GTSE1,BIK,RRM2,TCL1A,SHCBP1CDCA3,CDC20,TPX2,LRMP,CCNA2,MND1,CCNB1,PBK,ZWINT,RMI2PC_2Positive:CST3,LYZ,FCN1,CSTA,S100A9,S100A8,TYROBP,LST1,FGL2,VCANS100A12,SERPINA1,MNDA,FCER1G,CLEC7A,MS4A6A,CD14,CFD,IL1B,TYMPLGALS1,,AIF1,CTSS,NAMPT,CFP,TNFSF13B,CSF3R,MPEG1,TMEM176BNegative:IL32,NPM1,CD69,TRBC1,ISG20,ITM2A,IGKC,IGHA1,HSP90AB1,DDIT4HIST1H4C,PSIP1,AQP3,MYC,PIM2,HMGN1,PASK,NUCB2,HSPA1B,HSPB1CD79A,SUSD3,KLRB1,SYNE2,CHI3L2,IGHG3,IGLC2,FKBP11,IGHG1,SH2D1APC_3Positive:IL32,NKG7,CTSW,TRBC1,GZMA,CST7,GNLY,MKI67,ANXA1,TOP2ACCL5,PRF1,BIRC5,S100A4,KLRB1,CCNA2,AURKB,CENPF,GTSE1,CDKN3KLRD1,UBE2C,CDK1,TYMS,TPX2,RRM2,ID2,S100A6,FGFBP2,CDC20Negative:HLA-DRA,HLA-DQA1,HLA-DQB1,CD79A,HLA-DRB1,MS4A1,CD74,HLA-DPA1,HLA-DPB1,CD79BHLA-DMA,HLA-DMB,BANK1,VPREB3,IGKC,HLA-DRB5,MEF2C,CD22,IRF8,CD19SMIM14,FCRLA,HLA-DOB,CD24,CD40,FCER2,BLK,HLA-DQA2,IGHD,CTSHPC_4Positive:TOP2A,UBE2C,MKI67,GTSE1,CENPF,AURKB,PLK1,CCNA2,CDK1,CDCA8HMMR,CDCA3,CDC20,TPX2,CDKN3,DLGAP5,CENPE,BIRC5,CCNB2,CENPAKIF2C,CKAP2L,PBK,NUSAP1,KIFC1,AURKA,SPC25,NUF2,KIF23,ASPMNegative:NKG7,GNLY,CST7,GZMB,GZMA,PRF1,KLRD1,FGFBP2,CCL5,KLRF1HOPX,CTSW,GZMH,TRDC,FCGR3A,SPON2,CLIC3,MATK,ADGRG1,S1PR5CCL4,CMC1,XCL2,PFN1,CD160,FCRL6,IL2RB,TRGC1,KLRC1,C12orf75PC_5Positive:ICA1,PDCD1,TBC1D4,ITM2A,ICOS,MAF,TOX2,IL32,TNFRSF4,PASKPKM,SMCO4,ACTG1,CORO1B,CTLA4,NPM1,TRBC1,PCAT29,TIGIT,,ANP32B,ENO1,GBP2,COTL1,GAPDH,SUSD3,PIM2,AQP3,SERPINA9Negative:NKG7,GNLY,KLRD1,FGFBP2,GZMB,GZMA,KLRF1,CCL5,PRF1,TRDCGZMH,CST7,CTSW,BANK1,MATK,PLK1,HMMR,HLA-DPB1,CENPA,CLIC3GTSE1,CENPE,CCL4,SPON2,PDLIM1,HLA-DPA1,CDCA8,DLGAP5,TPX2,IGHD

拐点法寻找top可用的主成分用于后续分析(具体选择方式见：（重磅综述：三万字长文读懂单细胞RNA测序分析的最佳实践教程（原理、代码和评述））)

ElbowPlot()

降维可视化

＃我们选择top15PCs用于后续分析(,reduction="pca",dims=1:15)

聚类

(,reduction="pca",dims=1:15)(,resolution=0.5)

:15524Numberofedges:543084RunningLouvainalgorithmMaximummodularityin10randomstarts:0.9185Numberofcommunities:17Elapsedtime:2seconds

画图看一看！ggplot2高效实用指南(可视化脚本、工具、套路、配色)

p1-DimPlot(,reduction="umap",="")p2-DimPlot(,reduction="umap",="")p3-DimPlot(,reduction="umap",="RNA_snn_",label=T)+NoLeg()cowplot::plot_grid(p1,p2,p3)

让我们看看部分基因的表达情况！

FeaturePlot(,reduction="umap",features=c("CD3D","CD8A","CD4","CD14","MS4A1","FCGR3A","IL3RA"))

命名细胞簇并移除红细胞(cellassign：用于肿瘤微环境分析的单细胞注释工具（9月Nature）)

去除红细胞()[@$RNA_snn_=="16"][,!colnames()%in%]

现在，我们已经在数据中识别并命名了cluster，我们将继续进行celltalker分析。随该软件包一起提供的有一个ramilowski_pairs，它是一个由配体、受体和推测的配体受体对组成的。

首先，根据通用型配体和受体从我们的数据集中选出对应的基因，然后进行差异基因分析，只保留在样品组之间差异表达的配体受体。

然后，我们将为每个重复样本创建单独的数据矩阵，存储为tibble格式，以便于使用tidyverse进行后续处理。

(生信宝典注：如果我们自己有受体配体对，也可以整理成这样一个三列的格式，导入进来，替换掉数据包原有的配体受体对信息，实现更加定制的分析。)

该数据集中有2,557个特异的配体/受体对在我们的数据集中识别配体和受体(unique(ramilowski_pairs$ligand))(unique(ramilowski_pairs$receptor))()[rownames()%in%ligs]()[rownames()%in%recs](,)使用FindAllMarkers区分组之间差异表达的配体和受体Idents()-""markers-FindAllMarkers(,assay="RNA",features=,=TRUE)

(markers$gene)length()

[1]61

过滤ramilowski配受对_pairs[(ramilowski_pairs$ligand)%in%,]_pairs[(ramilowski_pairs$receptor)%in%,](,)dim()

[1]2413

生成celltalker的输入数据

(,slot="counts")@$@$@$_matrices(=,clusters=,groups=,replicates=,=)Atibble:2x2groupsampleschrlist1pbmctibble[3×2]2tonsiltibble[3×2]

数据展开为单个样品展示

unnest(,cols="samples")

=10：每个cluster中至少有10个细胞表达了配体/受体Atibble:2[8×2]2tonsiltibble[8×2]

unnest([1,2],cols="")

:只对包含细胞数大于总细胞数5%的簇进行互作分析_interactions(=,clusters=,groups=,=0.05,=)

Warning:`cols``cols=c()`Warning:`cols``cols=c()`

现在我们有了配体/受体相互作用的列表，我们可以使用unique_interactions函数鉴定组特异的互作，并使用circos_plot函数可视化组之间的差异。

([1,2])[[1]]([1,2])[[1]][]circos_plot(interactions=,clusters=)

PBMC组特有的受体-配体互作

从以上图中我们可以看出研究人员用黄色表示配体基因，绿色表示受体基因，然后将可以相互配对的基因连在一起构成簇之间的互作关系。

Tonsil组特有的受体-配体互作

#([2,2])[[1]]([2,2])[[1]][]circos_plot(interactions=,clusters=)

Tonsil组特有的受体-配体互作

CIRCOS圈图绘制-circos安装

CIRCOS圈图绘制-最简单绘图和解释

CIRCOS圈图绘制-染色体信息展示和调整

CIRCOS增加热图、点图、线图和区块属性

作者：May（协和医院）