人类胸腺发育的细胞图谱揭示了T细胞组库的形成

分享是一种态度

文章信息

文献标题：A cell atlas of human thymic development defines T cell repertoire formation 发表时间：2020.02.21 发表杂志：Science（IF=41.845） 原文链接：http://dx.doi.org/10.1126/science.aay3224

摘要

INTRODUCTION：胸腺（thymus）是 T 细胞发育和 T 细胞受体（T cell receptor, TCR）组库形成的重要器官，塑造了机体的适应性免疫。T 细胞的胸腺内发育有空间协调性，受胸腺微环境（thymic microenvironment）多种细胞类型的精细调节。尽管胸腺在多种动物模型中被广泛研究，目前尚缺少一份完整的人类胸腺图谱帮助我们理解人体免疫系统。

RATIONALE：作者将处于发育中、儿童期及成人期的胸腺标本解离后，利用单细胞转录组测序构建了跨越人类生命周期的完整胸腺图谱。作者总共采集了 15 个胎儿胸腺（从 7 - 17 孕周）和 9 个出生后至成人期的胸腺样本，同时采用了不同的细胞分选策略保证细胞亚群的覆盖率。根据 scRNA-seq 得到的 marker 基因，作者利用单分子荧光原位杂交（single-molecule fluorescence in situ hybridization, smFISH）对不同状态的细胞作了空间定位。另外，作者系统地比较了人和小鼠胸腺单细胞转录组数据的异同。最后，作者还获取了单细胞 TCR 组库数据以探索人类 TCR 重排和选择的时空偏倚。

RESULTS：作者从人类胸腺中鉴定出了超过 50 种状态的细胞。人胸腺中的细胞状态、丰度及基因表达谱随着胎儿及出生后发育呈动态变化。作者鉴定出了新的胸腺成纤维细胞和上皮细胞，并对其进行了定位。作者通过计算预测了人 T 细胞从胎儿肝脏造血干祖细胞到多种成熟 T 细胞类型的发育轨迹，进而构建了导致 T 细胞命运决定的转录因子网络框架。在非经典的胸腺 T 细胞中，作者注意到了一种 CD8αα+ T，以表达 GNG4 为特征，定位于胸腺髓质周围，它们高表达 XCL1，并且与 XCR+ 树突状细胞共定位。与小鼠数据相比，作者发现人体中这类非经典胸腺 T 细胞有着不同的基因表达谱。最后，作者发现人 VDJ 在重排和多轮选择下呈现出强烈偏倚，包括 CD8+ T 细胞中的 TCRα V-J 偏倚等。

CONCLUSION：该文章描绘了跨生命周期、跨物种的胸腺单细胞转录组图谱，提供了胸腺微环境下 T 细胞发育的高分辨率景观。通过系统地跨物种比较，作者揭示了人类特异性的细胞状态和基因表达特征。这份详细的细胞网络将有助于未来构建体外类器官模型，用以还原人类的在体胸腺组织。

测序数据

• 组织样本：
  1. 人胸腺样本：胚胎（7 - 17 PCW）、出生后（3 m - 15 y）、成人（25 y、35 y）
  2. 小鼠胸腺样本：C57BL/6J （4、8、24 w）
• 分选策略：
  1. DAPI 阴选富集活细胞
  2. CD45 分别富集免疫细胞和非免疫细胞
  3. CD3 分别富集 T 细胞和非 T 细胞
  4. EPCAM 富集上皮细胞
• 建库方法：
  1. 10x Genomics Single Cell 3' v2 and 5’ Reagent Kit：人类样本
  2. 10x Genomics Single Cell 3' v3 Reagent Kit：小鼠出生后样本
  3. 10x Genomics Single Cell 5' VDJ：部分人类样本
  4. Smart-seq2：部分人类样本
• 公共数据：
  1. scRNA-seq：人胚胎肝脏（E-MTAB-7407）、人骨髓造血干细胞（GSE117498）、小鼠胚胎胸腺（GSE107910）、小鼠胸腺基质细胞（GSE103967）
  2. Microarray：人出生后分选的胸腺细胞（E-MEXP-337）

数据分析

由于本文大量篇幅用在阐述 10x 数据，Smart-seq2 数据主要用于验证，这里我重点总结 10x 数据的分析方法。

原始数据处理

基因表达数据：Cell Ranger Single-Cell Software Suite (version 2.0.2 for 3’ v2 chemistry, version 3.0.0 for 3’ v3 chemistry and version 2.1.0 for 5’ chemistry)

VDJ 数据：Cell Ranger Single-Cell Software Suite (version 2.1.0)

下游分析

基本分析流程：

python3 下使用 scanpy 软件（version 1.3.4）

1. 归一化：scanpy.api.pp.normalize_per_cell
2. log 转换：scanpy.api.pp.log1p
3. 高可变基因（Highly variable genes, HVGs）：scanpy.api.pp.filter_gene_dispersion
4. 细胞周期：将预先定义的细胞周期相关基因从 HVGs 中移除
5. Scaling data：scanpy.api.pp.scale
6. PCA：scanpy.api.pp.pca 基于高可变基因
7. 数据整合：scanpy.api.pp.bbknn
8. 聚类：scanpy.api.tl.louvain

质控：

1. 保留 UMI > 2000、500 < Gene < 7000 的细胞
2. 鉴定 doublets：scrublet（参数：sim_doublet_ratio=2、n_neighbors=30、 expected_doublet_rate= 0.1），scrublet 评分 > 0.7 为 doublet
3. 过度预聚类（sc.tl.louvain 分辨率 20），平均 scrublet 评分 > 0.6 的聚类为 doublets。剩余的聚类中，若符合以下标准则同样推断为 doublets： (1) 根据先验知识，同时表达两种或以上细胞类型 marker 基因（例如同时表达 CD3 和 CD19） (2) UMI 数过高 (3) 缺少定义该聚类的 marker 基因
4. 移除甲状腺来源的污染细胞聚类：Epi_PAX8（PAX8, HHEX, TG, NKX2.1）和 Epi_GCM2（PTH, GCM2, GATA3, CHGA）

批次效应矫正：

考虑到批次效应可能来自技术差异（例如 10x 试剂版本）和生物学差异（例如发育阶段和组织供体），作者采用了迭代整合方式：

1. 利用 scanpy.api.pp.bbknn 整合相似样本（例如所有胎儿和儿童来源样本）
2. 利用整合后的数据获取粗略的细胞亚群注释
3. 利用批次信息（如 10x 试剂版本、供体等）或细胞亚群注释作为分类变量，拟合 L2-正则化的线性模型（L2-regularised linear model），将批次信息造成的差异回归掉（regress out），保留仅包含生物学信息的残差
4. 再次利用 scanpy.api.pp.bbknn 整合数据，用于下游的细化聚类、可视化和轨迹推断

聚类和细胞亚群注释：

1. 根据已知 marker 基因和差异表达基因（differentially expressed genes, DEGs）进行手动注释，首先将分群清晰、具有唯一定义的细胞聚类注释出来
2. 利用上述注释的聚类训练 logistic 回归模型（logistic regression model），用以预测因转录相似性而被软件聚到一起的混合细胞亚群
3. 对粗略注释的细胞分群（如上皮细胞、单阳性 T 细胞等），将其提取出来重复 HVGs 选择、降维和聚类流程，以获取更精细的细胞亚群（如髓质胸腺上皮细胞 mTEC 或调节性 T 细胞等）

估计细胞类型的比例：

1. 作者首先宽泛地定义细胞类型（如淋巴细胞、髓细胞等），并计算每种类型在选定的对比组之间的比例
2. 如果一次比较中的所有细胞类型均来自同一个分选门，则将比例简单定义为：特定类型细胞数 / 细胞总数
3. 如果细胞来自不同的分选门，则为每个分选门计算一个归一化因子：给定分选门的细胞数 / 所有分选门的总细胞数。对每个分选门的细胞数均乘以相应的归一化因子，再用归一化的细胞数来计算细胞比例
4. 用 t 检验评估细胞比例变化的显著性

轨迹推断：

作者利用前述批次矫正方法获得邻域图（neighbourhood graph），选择符合先验知识（包括已知 T 细胞分化、TCR 重排的顺序）的流形（manifold），然后利用 scanpy.api.tl.dpt 计算扩散拟时间（diffusion pseudotime），并计算拟时间中差异表达的基因

转录因子网络：

1. 人类转录因子信息下载自数据库 AnimalTFDB3
2. 基因表达估算：在 3D UMAP 空间对最近 30 邻域（30-nearest neighbors）取平均
3. 在 3D UMAP 的每个体素中随机抽样细胞，以去除冗余信息。在用于估算表达量的最近 30 邻域内计算细胞类型频率，对每个细胞类型计算注释评分
4. 选取高可变的转录因子计算相关性矩阵，建立图空间，并利用力导引图（force-directed graph, FDG）实现可视化

• VDJ 序列分析：
  1. 质控：选取包含全长的重排序列且 UMI 数 > 2 的 TCR 链
  2. 鉴定生产性重排（productive）的 TCR：能找到覆盖互补决定区（complementarity determining regions, CDR）CDR3 的开放读码框（Open Reading Frame, ORF）
  3. 比较 V、D、J 基因的使用率：对每个细胞类型的每个 VDJ 基因计数，利用该细胞类型的总 VDJ 计数进行归一化，并转换为 z-score，利用 t 检验在不同细胞类型之间比较 z-score。对 CD4+ T 和 CD8+ T 细胞的比较，还使用了 Cochran–Mantel–Haenszel 检验
• 细胞间相互作用：CellPhoneDB(www.CellPhoneDB.org)
  1. 每种细胞类型抽样 1000 个细胞
  2. 从 HUGO Gene Nomenclature Committee 数据库获取趋化因子信息
  3. 通过在细胞-细胞对内将配-受体对的平均表达水平相乘，计算相互作用评分，并使用最大值归一化
  4. 根据特异性的显著程度选取相互作用对（interaction pairs）进行可视化
• 整合公共数据： 同样使用了 BBKNN 方法
• 跨物种比较：
  1. 使用前述流程对小鼠数据集做预处理，包括归一化、批次整合、降维聚类等。利用人类数据集训练 logistic 模型对小鼠数据进行辅助注释
  2. 数据整合：结合了两种方式： (1) 基于互为训练集的 logistic 模型做双向预测，两套预测概率相乘作为最终的相似度评分 (2) 线性回归（linear regression）辅助的 BBKNN 方法去除物种特异性变异，保留生物学信息的结构

主要结果

人类生命过程中胸腺内的细胞组成

• 质控后保留了 255,901 个单细胞，包含了从胎儿期、幼年期、青春期到成人期的完整胸腺图谱
• 鉴定出超过 40 种细胞亚型：
  • 免疫细胞：包括分化的 T 细胞（DN、DP、CD4 / CD8 SP、Treg 、CD8αα+ T 和 γδ T）、B 细胞、NK 细胞、固有免疫细胞（ILCs）、巨噬细胞、单核细胞和树突状细胞（dendritic cells, DCs）等
  • 非免疫细胞：胸腺上皮细胞（thymic epithelial cells, TECs）、成纤维细胞、血管平滑肌细胞（vascular smooth muscle cells, VSMCs）、内皮细胞和淋巴管内皮细胞等
• DCs 进一步被分为经典的髓样 DCs（DC1、DC2）和浆细胞样 DC（plasmacytoid DC, pDC）
• 成纤维细胞细被细分为 Fb1（COLEC11, C7, GDF10）、Fb2（PI16, FN1, FBN1）和 Fb_cycling。Fb1 表达固有免疫相关的重要基因 COLEC11 以及调节上皮细胞发育的维甲酸反应酶 ALDH1A2。而 Fb2 表达细胞外基质基因和信号素（semaphorins），调控血管发育。smFISH 发现 Fb1 主要定位在胸腺小叶周围，而 Fb2 主要在小叶间区，并且和大血管相联系。
• 作者将 TECs 根据髓质（medullary）和皮质（cortical）分为 mTEC 和 cTEC，并利用小鼠数据集辅助注释，鉴定出 cTEC（PSMB11）、mTEC(I)（KRT14）、mTEC(II)（AIRE）、mTEC(III)（KRT1）以及罕见的 mTEC(IV)（DCLK1, POU2F3）。作者发现 cTECs 在早期发育阶段丰度较高，而在胎儿晚期和儿童期则出现一种中间态亚型 mcTEC（DLK2）。
• 作者还注意到两群 EpCAM+ 细胞：肌样 TEC（MYOD1, MYOG）和神经内分泌样 TEC（NEUROD1, NEUROG1, CHGA），定位于胸腺髓质。自身免疫性重症肌无力相关的基因 CHRNA1 在这两群细胞以及 mTEC(II) 中高表达。
• 最后，作者分析了已知的先天性 T 细胞免疫缺陷（congenital T cell immunodeficiencies）致病基因表达谱，揭示了这些基因可能发挥作用的位置和发育阶段

胸腺基质和 T 细胞的协调发育

• 早期胎儿胸腺（7 - 8 PCW）中，淋巴细胞主要包括 NK 细胞、γδ T 细胞和 ILC3；分化中的 αβ T 很少，主要处于 DN 阶段。随后 DN 逐渐分化出 DP 和 SP 阶段，并在 12 PCW 达到均衡；相反，固有淋巴细胞比例则逐渐减少
• 逐渐退化的成人胸腺中出现终末分化的 T 细胞，可能是来自归巢 T 细胞或外周血 T 细胞。退化胸腺中，细胞毒性 CD4+ T 细胞（CD4+ CTLs）、记忆 T 细胞和记忆 B 细胞比例均增加
• T 细胞发育的趋势和胸腺基质细胞相对应。随着 T 细胞成熟的启动，TEC 从 cTEC 主导逐渐向 cTEC 和 mTEC 平衡状态转变，提示了胸腺上皮细胞和成熟 T 细胞相互作用、同步分化的特征
• 成纤维细胞中，Fb1 在早期发育阶段占主导，而后 Fb1 和 Fb2 逐渐达到均衡
• 其他免疫细胞中，巨噬细胞在早期发育中较丰富，而 DCs 随着发育进程逐渐增加。其中，DC1 在 12 PCW 后占主导，而 pDC 在出生后迅速增加
• 作者最后预测了胸腺基质细胞和 T 细胞相互作用的配-受体对：淋巴毒性信号（LTB:LTBR）来自多种免疫细胞，并由大部分基质细胞接受；而 RANKL-RANK 信号（TNFRSF11:TNFRSF11A）局限于 ILC3 和 mTEC(II) / 淋巴管内皮细胞；FGF 信号（FGF7:FGFR2）从成纤维细胞到 TEC，且 FGFR2 在成人胸腺表达降低；NOTCH1 是早期胸腺祖细胞（early thymic progenitors, ETPs）的主要受体，而多种 Notch 配体中，JAG2 和 DLL4 主要由 cTECs 和内皮细胞表达，其他 TEC 则广泛表达 JAG1

传统 T 细胞分化轨迹

• 作者整合了胎儿肝脏造血干祖细胞和胸腺 T 细胞的数据，发现 ETP 在 UMAP 图上位于造血干细胞 / 多能干祖细胞（hematopoietic stem cells / multipotent progenitors, HSCs/MPPs）和 pre / pro-B 细胞之间。而进一步整合骨髓单细胞数据后，作者发现 ETP 位于骨髓多重淋巴祖细胞（multi-lymphoid progenitor, MLP）和胎儿肝脏早期淋巴祖细胞（early lymphoid progenitor）之间

• 作者进一步推断了 T 细胞后续分化的轨迹：始于 CD4-CD8- DN，逐渐变为 CD4+CD8+ DP，经过 CCR9 high 的 αβ(entry) 阶段，分化为成熟的 CD4+ 或 CD8+ SP。同时，作者注意到 γδ T 在 DN-DP 交界处另行分化。
• 对 DN 和 DP，作者根据其细胞周期增殖（CDK1）将其细分为增殖态（proliferating）和静息态（quiescent）。VDJ 重排基因（RAG1, RAG2）从增殖晚期开始增加，于静息态达到高峰，提示了 T 细胞分化中增殖和重排的关系
• 通过整合 TCR 数据，作者发现，在增殖晚期 DN 中检测到重排的 TCRβ 链，与重排事件和 pre-TCRα（PTCRA）表达的增加相符合。TCRβ 链非生产性 / 生产性重排的比率在 DN 中相对较高，而在进入 DP 阶段逐渐下降到基线水平，提示了 β-selection 的影响。在 DN(Q) 阶段，TCRβ 链非生产性重排达到高峰，可能提示第一轮等位基因重排失败的细胞进行了其他等位基因重排。在 DP 阶段，自增殖期开始检测到重排的 TCRα 链，且与 TCRβ 链不同，其在 DP(Q) 阶段非但没有富集，反而直接被丢弃了
• 通过推断拟时间差异表达基因，作者鉴定了新的 T 细胞发育阶段及 marker 基因，例如early DN（ST18）、DP （AQP3）和 DP-to-SP 过渡态（TOX2），并且构建了 T 细胞谱系定向相关的转录因子调控网络

• 最后，作者结合基于蛋白 marker 分选的 microarray 公共数据，发现 DN(P)、DN(Q) 和 DP(P) 分别对应 CD34+CD1A+ DN、CD4+ ISP 和 CD3- DP，而 DP(Q) 和 αβ(entry) 对应 CD3+ DP。其中 DN(Q) 和 CD4+ ISP 均表达 pre-β-selection 特征

Treg 和非经典 T 细胞的发育

• 作者鉴定出了多种非经典 T 细胞，例如 Treg（FOXP3）、CD8αα+ T、NKT 样细胞（EOMES, NKG7, IFNG, TBX21）以及 TH17 样细胞（CD4, CD40LG, RORC, CCR6）。这些细胞的发育需要激动剂选择（agonist selection），其非生产性 TCR 链比例很低，提示在胸腺内的时间比经典 T 细胞更久
• 非经典 T 细胞在胸腺的富集程度高于肝脏，并且均富集于成熟后（~10 PCW）的胸腺，提示其发育具有胸腺依赖性
• 在 UMAP 图中，Treg 和 αβ T 有分化轨迹联系。作者将连接处的 Treg 定义为 Treg(diff)，其表达低水平 FOXP3 和 CTLA4，而高表达 IKZF4、GNG8 和 PTGIR，这些基因均与自身免疫和 Treg 分化相关。此外，作者注意到另一群由非编码 RNA（MIR155HG）定义的不同于终末分化 Treg 及 Treg(diff) 的细胞，注释为 T(agonist) ，其表达 IL2RA 但低表达 FOXP3，与小鼠胸腺中 CD25+FOXP3- Treg 祖细胞相似。Treg(diff) 和 T(agonist) 可能代表了人类胸腺中的两种 Treg 祖细胞
• 作者进一步将 CD8αα+ T 细胞分为 3 群：CD8αα+ T(I)（GNG4, CREB3L3, CD72）、CD8αα+ T(II)（ZNF683, MME）和 CD8αα+ NKT 样细胞（EOMES）。CD8αα+ T(I) 和 CD8αα+ T(II) 在早期高表达 PDCD1，而在终末分化阶段逐渐下调。CD8αα+ T(I) 和晚期 DP 及 αβ(entry) 之间有明显的谱系区分，CD8αα+ T(II) 则混合了 αβ 和 γδ T 的特征，而 NKT 样细胞则富集 γδ TCR 链
• 利用 smFISH，作者发现 GNG4+ CD8αα+ T(I) 细胞位于胸腺髓质。CD8αα+ T(I) 和 Treg 都表达表面蛋白 marker CD137（TNFRSF9），于是作者通过流式分选和 Smart-seq2 建库测序，验证了 CD3+CD137+CD4- 可以作为 CD8αα+ T(I) 的分选策略。跨物种比较发现，人 CD8αα+ T(I) 与小鼠上皮内淋巴细胞前体 A 型（intraepithelial lymphocyte precursor type A, IELpA）最为相似，均表达 HIVEP3、NR4A3、PDCD1 和 TNFRSF9，但同时也存在物种间差异基因（人类 GNG4 和 XCL1，小鼠 ZEB2 和 CLDN10），并且人类 CD8αα+ T(I) 最终成为 CD8AhighCD8Blow 表型，而小鼠 IELpA 则表现为 CD8AlowCD8BlowCD4low 特征
• NKT 样细胞和 TH17 样细胞均表达 KLRB1 和 ZBTB16，表现出类似固有淋巴细胞的特征
• 最后，作者分析了 TRDV1 和 TRDV2 两个常见的 TCRδ V 基因，发现 TRDV2 主要用于 DN 阶段，而 TRDV2 主要用于 DP(Q) 和 αβT entry 阶段。进一步地，作者推断 CD8αα+ T(II) 主要来自晚期 DP，而 NKT 样和 TH17 样细胞主要来自早期分化的胸腺细胞

DCs 在胸腺细胞选择中的募集和活化

• DC 亚群及 marker 基因：前文所述的 DC1（XCR1, CLEC9A）、DC2（SIRPA, CLEC10A）和 pDC (IL3RA, CLEC4C），以及新鉴定的活化 DC（activated DCs, aDCs; LAMP3, CCR7）。aDC 高表达多种趋化因子、共刺激分子和转录因子（AIRE, FOXD4），它们可能代表了过去在人胸腺和扁桃体描述过的 AIRE+CCR7+ DCs
• aDC 进一步细分为 aDC1、aDC2 和 aDC3。作者发现 aDC1 与 aDC2 分别同 DC1 及 DC2 有相似表达谱，提示它们之间的激活关系（aDC1 - DC1 和 aDC2 - DC2），并且两群 DC 分别表达不同的趋化因子。aDC3 细胞中，主要组织相容性复合物（major histocompatibility complex, MHC）和共刺激分子的表达均相对降低，提示其可能为 DC 的活化后状态
• 作者重点关注了趋化因子介导的细胞间相互作用，发现 aDC 同时表达 CCR7 和 CCL19，提示其具有募集 T 细胞进入胸腺髓质的能力。aDC 还表达 CCL17 和 CCL22，相应的受体 CCR4 则由 CD4+ T 和某些特定的 Treg 表达。aDC 还可能通过 CXCL9/10:CXCR3 募集其他 DC 和成熟 Treg，并在 Treg 成熟过程中提供共刺激信号。作者还注意到 CD8αα+ T(I) 表达 XCL1，可能和 DC1（XCR1）的募集与活化有关
• 最后，作者利用 smFISH 验证了上述发现

人 TCR 组库形成和选择的偏倚

• TCRβ 的 VDJ 基因重排偏倚从 DN 到成熟 T 细胞的过程中持续存在，该偏倚可能和基因座的襻环结构有关。D2 基因与 J2 基因之间具有偏好性，而 D1 基因和 J1、J2 重排的频率差异不大。V - D 或 V - J 之间则没有明显相关性。在 β-selection 后，某些 TRBV 基因被丢弃或富集，提示在不同 Vβ 基因对 MHC 肽刺激的反应能力具有种系编码的差异
• TCRα 基因座中，发育时间和 V - J 配对之间有清晰的相关性：近端先于远端重排，V - J 基因严格配对。而在成熟 T 细胞中，远端 V - J 配对被丢弃，可能反映了阳性选择（positive selection）阶段的偏倚
• 作者比较了不同细胞亚群之间的 TCR 组库差异，发现 CD8+ T 细胞和其他细胞有明显差别，其 TRAV - TRAJ 组库偏好远端 V - J 对，提示这种差异可能与 CD8+ T 谱系的缓慢或低效定向有关。和胚胎期相比，出生后胸腺中的 CD8αα+ T(I) 则偏好近端的 V - J 对

总结

本文构建了人类经典和非经典 T 细胞的分化轨迹以及 TCR 组库信息，发现了 TCR 组库在成熟 T 细胞中的偏倚，可能提示了人体如何应对不同的抗原挑战。本文对胸腺微环境的分析揭示了构成胸腺的复杂细胞类型，以及基质细胞和固有免疫细胞如何相互作用以支持 T 细胞分化。这些细胞间交流网络可以帮助体外构建产生 T 细胞的培养系统，并影响未来 T 细胞治疗的工程学策略。

我的评价

本文的生物学意义已经说得很清楚了，我来从生信分析的角度谈一些感想吧。

面对这样大量的单细胞样本，无论是数据读写还是内存占用，光靠基于 R 语言的分析平台已经远远不够。人类发育细胞图谱（HDCA）项目里，过去两篇母-胎界面（70,000 细胞）和胎儿肝脏造血（大约 130,000 细胞）的 Nature 基本还是以 Seurat 包为主，本文可能是该项目第一篇纯 python 分析的大文章吧。本文的代码已经以 jupyter notebook 的形式上传 github，包含大量经典套路和高级分析，其中 logistic 模型辅助细胞注释的部分感觉是 Teichmann 团队后来发表的 SCCAF(https://github.com/SCCAF) 方法的前身。因此，本文理论上说可以作为学习 python 环境下分析单细胞数据的案例。

于是乎，我们又来到了“可复现的生物信息学分析”这个话题。前两年，一位中科院计算所的老师曾告诉我：”即使公开了代码，生信分析的文章仍有 70 - 80 % 都不能完全重复出来。“然而，早在 2017 年，《Reproducible bioinformatics project: a community for reproducible bioinformatics analysis pipelines》(https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-018-2296-x)一文已经指出以 docker 为代表的容器化技术可能为我们带来更灵活可控的生信分析环境。理论上说，使用给定镜像的容器在任何设备上做分析，只要喂给计算机相同的数据，理应得出相同的结果。除此之外，各大杂志也在利用 github、codeocean 等代码托管或云计算平台推进生信分析规范化的进程，上文提到的 SCCAF作者同样提供了 docker 镜像用于重复该软件的结果。综上，未来的生信分析文章一定会面临更加严苛的检验。

话说我真的不是给 docker 打广告啊。醒醒，copy一下教程、随便跑跑代码出一堆图、水得一手文章的日子快要到头啦！

如果你对单细胞转录组研究感兴趣，但又不知道如何入门，也许你可以关注一下下面的课程