推荐理由
这篇文章涉及单细胞的细胞类型鉴定、拟时分析、细胞交流等主流的单细胞分析方式，提供数据及代码，因而学习和参考价值较高。
卡路里限制重编程老化小鼠的单细胞转录景观
老化导致的全身组织功能减退或许可以通过卡路里限制（caloric restriction ，CR)推迟,但但是CR造成的细胞、衰老信号的变化以及减弱衰老的机制仍不清楚。在此我们对衰老以及CR小鼠建立不同组织的单细胞以及单细胞核转录组图谱。我们发现CR减弱了衰老相关的变化，包括细胞组成类型、基因表达、核心转录调控网络。在衰老过程中免疫细胞增多，CR有利于逆转这种衰老造成的免疫微环境。计算预测揭示了在衰老过程中观察到的异常细胞间通讯模式，包括过多的促炎配体与受体的相互作用，这些现象可以被CR逆转。我们的工作提供了衰老和CR动物的多组织的单细胞转录图谱，增强了我们对CR这种强有力的干预保护措施的理解，揭示了代谢干预可以通过免疫系统改变衰老的进程。
1.CR对寿命以及衰老特征的影响
实验方法：最初小鼠都采用随意饮食（ad libitum feeding，AL），CR开始于中年（取部分18个月大的AL小鼠）限制饮食为正常的70%。进过9个月的饮食干预，27个月大的小鼠（相当于人类70岁）可分为O-AL、O-CR。这两组小鼠都使用5个月大的年轻AL小鼠（Y-AL）做对照。
衰老特征以及寿命延长：
O-CR小鼠的血糖相比于O-AL小鼠没有变化，但是体重减轻了约1/3(32.6% ± 3.9%）
-
O-CR的中位生存时间以及最长寿命都有所增加
-
CR影响组织中年龄相关的特征变化。通过肝组织染色（Oil
Red O staining）发现Y-AL和O-CR相比于O-AL有更少的脂肪聚集。
此外在O-AL小鼠的棕色脂肪组织 （brown adipose tissue ，BAT）中找到较少的多叶小脂肪细胞和较多含有较大脂滴的米色脂肪细胞，CR可以逆转这一现象。并且CR可以减少白色脂肪组织（white adipose tissue，WAT）中脂肪细胞的体积。
-
老化相关的β半乳糖苷酶（SA-b-Gal）染色显示老化细胞在WAT、BAT、肝组织和肾组织累积，这种现象可以被CR逆转。
2.构建小鼠单细胞图谱
使用液滴微流控的测序方法，质控后获得166,111个细胞用于后续分析。
最终得到109个细胞群，隶属于7种组织42类细胞。
3.卡路里限制重塑衰老的细胞的微环境
通过分析rescue（挽救）实验发现一些组织的细胞类型的变化是可逆的（通过卡路里限制逆转衰老）
在多种组织的免疫细胞（Ptprc+ cells）（除了骨髓）都是随年龄增长儿增多，而CR可以阻止这种趋势。
这种现象在中性粒细胞以及浆细胞中尤为明显。
上文rescue实验结果的免疫印迹也验证了这一结论。
并且定居外周组织的中性粒细胞低表达骨髓的中性粒细胞前体标记(Mki67+/Lcn2+) ，但是高表达趋化因子（Cxcl2, Ccl2- (MCP-1), and Cxcr4）。拟时间分析证实中性粒从骨髓迁移到外周组织，CR可以逆转这一过程。
CR对巨噬细胞的状态也有调节作用
4.CR逆转了多个组织中衰老相关基因的表达
O-AL 和 Y-AL之间的差异基因定义为“衰老基因”，O-CR 和 O-AL之间的差异基因定义为“CR基因”，衰老基因和CR基因的重合部分（即表达可以被CR挽救的基因）称为“挽救基因”。
玫瑰图显示主动脉、褐色脂肪组织以及白色脂肪组织受这些差异基因影响最大。
但是也存在一些下调衰老基因被CR进一步下调，或者一些上调的衰老基因被CR进一步上调。但是这些基因仅占所有衰老基因的1.8%，他们产生的效应或许不足为虑。
我们对“衰老基因”、“CR基因”以及“挽救基因”进行了GO分析，下调的基因主要集中在血管发育、细胞再生、细胞对生长因子应答，细胞外基质通路上。并且这种下调可以被CR逆转。衰老同时也加强了炎症应答，这种状况同样可以被CR逆转。
46 个下调的衰老基因 在至少3中组织中可以被CR挽救, 并且34个上调的衰老基因在至少3中组织中可以被CR抑制。其中
Ybx1 (freq = 5), Cebpb (freq = 5), Klf4 (freq = 3), and Atf3 (freq = 3)是重要的转录因子，S100a9 (freq = 5), Igkc (freq = 5), S100a8 (freq = 4), Cxcl2
(freq = 3), Il1b (freq = 3)是重要的炎症因子。
5.CR逆转年龄相关的细胞类型特异的基因表达
总的来说主动脉定植的增值细胞、中性粒细胞、自然杀伤细胞在所有组织中老化都是最严重的。
-
主动脉定植的增值细胞，巨噬细胞M2、和表皮细胞被CR挽救的程度最大。
不同组织的不同类型细胞中三类差异基因的表达情况 -
不同组织中的不同类型的细胞受老化影响，以及CR后得到挽救的程度也不同。
不同组织细胞中差异基因的情况
S100a9, S100a8 和 Igkc在多达40种不同类型衰老细胞中表达上调，并且可以被CR挽救。Ybx1则是在衰老细胞中表达上调，我们通过在人的白色脂肪细胞、脂生成干细胞（adipose-derived stem cells ，ASCs）中敲除Ybx1。该基因的敲除使这类干细胞逐渐丧失了增殖能力，与干细胞的体内耗竭一致。
随后对敲除了Ybx1的ASC进行了RNA测序GO分析，发现下调的基因与细胞分化以及DNA复制有关——这与我们免疫荧光观察ASC增殖能力的实验结果一致。
敲除Ybx1的ASC
RNA-Seq
Yxb1可能是“衰老-CR挽救”之间重要的分子转换器。
另一个问题是：受衰老和CR调控的经典的衰老标记和营养感受通路是否也存在单细胞水平上细胞类型的差别的调控模式？
6.CR重建转录调控网络
使用SCENIC预测（衰老&CR）核心转录调控因子
在衰老时表达下调，并且可以被CR恢复的转录因子称为“rescue TFs"，在不同组织中我们鉴定出来不同的rescue TFs
其中亮氨酸拉链蛋白Cebpd被证实与炎症和脂代谢相关，Nfkb1 、Rel与验证相关。
7.通过计算预测基于CR清除细胞间交流模式异常
我们对每种组织内不同细胞类型间进行了“配体-受体”配对，也考虑了自分泌，对细胞间交流进行了鉴定。
其中我们发现衰老普遍增强了细胞间交流，CR可以逆转这种趋势。
然后具体分析了那些可以被逆转的途径（对受体-配体进行GO分析）发现这些可以被消除得到联系通路都是“细胞因子-细胞因子受体”，这与CR逆转炎症的作用一致。
8.在衰老和CR过程中细胞类型以及转录组在两性的变化
在衰老过程中细胞类型分布、差异基因在两性中都不存在显著差异。但是“挽救基因”在雄性大鼠中的比例要更高，这暗示饮食干预或许对雄性的作用更大。
但是这里都没有放统计学检验的P值？？？
免疫细胞，尤其是中性粒细胞在雌雄大鼠中增加显著。
9.衰老和CR过程中脑和骨骼肌细胞的单细胞核转录组测序
大脑和肌肉都是受衰老影响很大的器官，但是神经元与骨骼肌的分离都十分困难，在这里我们使用单细胞核测序对我们的上述研究做一补充（我们上文研究得到其中组织中不包括脑和肌）。
我们从三群小鼠（Y-AL，O-AL，O-CR）脑和骨骼肌中分别分离得到25,021 和 28,131个细胞核。,通过鉴定marker基因在脑组织中发现10种细胞、在骨骼肌中发现11中（C图）
与我们之前的结果类似，大脑和骨骼肌中的脂肪组织、表皮组织以及脑中的神经元都是随衰老而衰退，但是可以被CR挽救。然而肌卫星细胞呈现不可逆的衰退，不可被CR改善（图E）。
用到的软件及算法
image.png
image.png
image.png
methods
测序
仪器：Chromium Single-Cell instrument (10x Genomics)中文网址
建库：10x Genomics protocol using a Chromium Single-Cell 30 Gel Bead and Library V2 Kit
测序平台： NovaSeq 6000 Sequencing System (Illumina 20012866)
数据处理
利用bcl2fastq (version 2.20.0.422)转换格式（Illumina测序结果BCL files to FASTQ format）
利用10×Genomics的官方软件cell Ranger将reads进行比对以及计数version 2.2.0
测序结果
确定细胞类型
使用DoubletFinder 和 Seurat R 包完成质控、归一化、批次监控以及聚类。根据经典的细胞特异的marker的表达量确定细胞类型。使用 vignette of DropletUtils (version 1.6.1) package
处理10×Genomics droplet鉴定含有细胞的微滴以及空的背景。数据过滤标准
1.对于scRNS-seq舍弃少于500个基因或者线粒体基因的比例高于10%的细胞（肾细胞的线粒体比例低于50%即可），对于snRNA-seq舍弃少于300基因的细胞。
- Doublets(指的是一个液滴中出现一个以上细胞，对测序以及后续细胞类型鉴定造成困难）
Doublets用DoubletFinder package 鉴定，每个样本的均值方差标准化双峰系数（mean-variance-normalized bimodality coefficient ，BCMVN) 用于决定pK值，pN设定为0.25，之后用‘‘doubletFinder_v3’’功能求每个样本的doublet.3.使用‘‘NormalizeData’’ 功能对每个细胞进行log标准化，log(1 + counts per 10,000)。
4.使用‘‘FindVariableGenes’’ 功能寻找差异基因。
5.在样本间使用经典相关性分析（canonical correlation analysis ,CCA) 将差异基因投射到低维度以矫正批次效应，用‘RunMultiCCA’’ 和‘‘AlignSubspace’’ 函数。
利用欧氏距离作图，在低维度空间展示显著的主成分。细胞利用‘‘FindClusters’’功能进行聚类。
用‘‘RunTSNE’’功能进行2个维度上的t-SNE可视化展示。
8.每个类群的差异表达分析利用了‘‘FindAllMarkers’’功能的Wilcoxon秩检验。
- 细胞类型的聚类利用已知的marker基因的丰度。
细胞组成类型分析
排除了少于50个细胞的细胞类型（大动脉和中性粒细胞），在三种小鼠中（ (Y-AL, O-AL, and O-CR) 每种类型的细胞都被单独计数并且计算比例。O-AL and Y-AL细胞群之间细胞比例的变化以Log2FC展示衰老的影响(|Log2FC| > 0.5)。O-CR and O-AL 以Log2FC展示挽救对衰老的影响。
拟时间分析
拟时间分析使用筛选过的数据集
- 筛选条件
来源于骨髓的前体粒细胞、中性粒细胞，其他组织中的外周粒细胞。
该组织中粒细胞数不少于30.使用Monocle2进行基因排序，表达该基因的细胞不少于10个（cutoff of expression>=10），内部差异表达以及离散度的q-value cutoff < 0.1。轨迹利用二维展示（DDRTree降维算法）。
差异表达和差异基因网络分析
对样本进行筛选 (O-AL/Y-AL and O-CR/O-AL)各对照组间样本少于三个细胞则剔除（排除了 B cells, CD8+ T cells and plasmocytes in the aorta, B cells, NK cells and neutrophils in the kidney, plasmocytes in BAT, and NKTs in skin），不同样本(O-AL/Y-AL and O-AL/O-CR)以及不同类型细胞的差异分析使用 Seurat package的‘‘FindMarkers’’功能的Wilcoxon秩和检验。
- First, DEGs between the O-AL and Y-AL groups were identified
to generate an aging-related DEG dataset (aging DEGs) (|LogFC| > 0.5, adjusted p value < 0.05)- Secondly, DEGs between the O-CR and O-AL groups were identified to generate a CR-related DEG dataset (CR DEGs) (|LogFC| > 0.5, adjusted p value < 0.05).
- 基于前两个结果 ‘‘rescue DEGs’’ 被定义为在aging DEGs和CR DEGs中变化趋势相反的基因，DEGs (|LogFC| > 0.5, adjusted p value < 0.05) between the O-CR and Y-AL groups were identified to compare the rescue effect between the aging and CR DEGs.
接下来整合差异基因（DEG)的组织和细胞类型信息以及相关基因的LogFC值，构建 cell-DEG networks 。利用Cytoscape (version 3.7.1)可视化结果GO分析
使用Metascape , 结果利用 ggplot2 R package (p < 0.01)可视化 。10 GO terms 与 aging and/or CR 相关，热图 和 基因表达谱聚类图 使用 pheatmap R package。
转录因子-靶基因网络分析
转录因子结合结构域使用GENIE3 R/Bioconductor packages以及RcisTarget database of the SCENIC预测，工作流程以及参数设置参照(http://scenic.aertslab.org/)。
GENIE3用于从基因表达矩阵中推断调控网络，RcisTarget用于识别富集的转录因子结合结构域以及靶基因，高度可信的注释过的转录因子及靶标用Cytoscape (version 3.7.1)展示。
细胞间交流分析
基于单细胞数据使用CellPhoneDB software进行预测，针对某一对细胞量类型当受体-配体在超过10%的细胞中出现才纳入下一步分析，否则认为细胞间联系不存在。受体-配体对的平均表达量在不同的细胞间进行比较，p < 0.05 纳入下一步三种小鼠间细胞通讯的分析 (Y-AL, O-AL, and O-CR) 。
bulk RNA-seq分析
首先质控去接头， 使用HISAT2将reads比对到 UCSC human hg19 genome 。 The featureCounts 用于计数注释过的基因。使用 DESeq2 找差异基因 cutoff Benjamini-Hochberg adjusted p <0.05 ，|Log2FC| > 0.58. 使用 DESeq2 regularized logarithm (rlog)-标准化之后的数据计算样本之间的欧氏距离，来看样本之间的相关性。基因表达水平使用FPKM标准化。
统计学分析
数据使用 means ± SEM，GraphPad 6 Software双尾t-test 判断不同实验组和对照组之间的差别 p < 0.05 认为有显著差异。 *, **, *** and **** indicate p < 0.05, p < 0.01, p < 0.001 and p < 0.0001， NS indicates
not significant.
数据来源
GEO: GSE137869. https://github.com/MSZaki/Aging-and-CR-rat-atlas.
