癌睾丸抗原（cancer-testis antigen，CTA）是一类特异性表达于睾丸和肿瘤中的抗原，已应用于肿瘤治疗^［1-2］。由于CTA具有局限表达于睾丸且在肿瘤中高表达特征，有研究者采用生物信息分析等方法鉴定癌睾丸基因（cancer-testis gene，CTG）^［3-4］。这些技术不能确定基因表达的蛋白质是否具有免疫原性，故这些基因连同之前的CTA合称为CTG。在2009年，研究者分析了已鉴定的CTG，将其信息整理成数据库CTdatabase^［3］。CTdatabase数据库中所有内容都经过了相关领域内专家的注释，包括基因名称、基因组位置等。此外，研究者根据CTG是否分布在X染色体上，将CTG分成2类：分布在X染色体上的CTG（CT-X）和分布在常染色体上的CTG（Non CT-X）^［5-6］。

肿瘤相关基因特异性表达的研究，使人们注意到精子和肿瘤这2种生理与病理上截然不同的组织，其发生过程可能存在相似之处^［5，7］。研究者认为，肿瘤中CTG的异常表达，可能反映了体细胞中本该沉默的代表了精子发生过程的基因表达程序（gene expression program，GEP）被再次激活，并且这些GEP可能参与肿瘤的发生和发展^［7］。然而，CTG在精子发生和肿瘤发展过程中的GEP特征，及其与肿瘤患者预后的关系等问题，尚未得到解决。

本研究采用生物信息学技术，研究CTG的分布特征和表达规律。通过多个数据库的联合分析，重新对CTG进行筛选；分析睾丸精子发生过程中的基因表达谱，鉴定出与CTG相关的GEP，并分析其在各种肿瘤类型中的活跃程度及其与患者预后的关联性。

GTEx（Genotpye-Tissue Expression）项目致力于提供给研究者用于研究遗传变异与基因表达调控的资源。本研究中使用的是移除个体水平身份信息的基因表达数据，即正常组织的转录组数据（https：//www.gtexportal.org/home/datasets）。其发布的时间是2019 年8月20日，版本号是V8。共有54种组织类型，包含睾丸组织。

Illumina Human Body Map（HBM）项目包含16种组织类型的转录组数据。为了与GTEx数据库的转录组数据相匹配，从HBM下载原始的测序短读测序片段（reads）数据（FASTQ格式存储），并用GTEx数据库的分析流程对HBM的数据进行比对和基因定量。HBM的FASTQ原始数据存储于EMBL-EBI数据库，访问号为E-MTAB-315。

通过TCGAbiolink^［8］检索TCGA数据库中33种癌症类型样本的转录组，并将用于gdc-client下载的索引保存至gdc_manifest.txt文件中。之后，用gdc-client下载33种癌症类型样本的转录组。本研究中使用的数据版本为v33.1，发布时间为2022年5月31日。

使用基因特异性度量（specific measure，SPM）衡量基因在睾丸组织中的特异表达情况^［4］。具体定义如下：基因在不同组织中的表达矩阵表示为X∈Rm×n，共有m个基因，n个不同的组织。基因i在不同组织中的表达量可以表示为向量xi∈Rn，如式1：

其中，i表示第i个基因，xij表示基因i在 j组织中的表达量，共有n个不同的组织类型。

基因特异性表达是指基因仅在某个组织中表达，而不在其他组织中表达。因此当基因i 特异地在组织 t中表达，那么基因i在不同组织中表达量的向量可以表示为式2：

根据式1和式2，可以定义基因i 在t 组织中的特异表达程度，即基因i在不同组织中的真实表达情况 xi和基因i特异表达于组织t时在不同组织中的理想表达情况xit 之间的cos相似度，SPMit。其定义如式3所示：

CTG的主要特征是在正常的睾丸组织中高表达，并且在肿瘤中也高表达。因此，可以分成两步筛选：先筛选睾丸特异表达基因（testis specific gene，TSG），再基于TSGs在肿瘤中的表达情况筛选出CTG。

使用正常转录组数据（GTEx数据库和HBM数据库）和肿瘤转录组数据（TCGA数据库）筛选CTG。首先，筛选TSG，具体步骤如下：①计算GTEx数据库中每个基因在不同组织类型中的平均表达量。②计算GTEx数据库中每个基因在睾丸组织中的特异性，即SPMm×1Testis,GTEx。③计算HBM数据库中每个基因在不同组织类型中的平均表达量。④计算HBM数据库中每个基因在睾丸组织中的特异性，即SPMm×1Testis,HBM。⑤筛选出满足SPMm×1Testis,GTEx>0.9和SPMm×1Testis,HBM>0.9，并且在睾丸组织中的平均表达量高于1 TPM的蛋白编码基因。其次，CTG还需要TSG满足在肿瘤中高表达（至少在1种肿瘤的1%样本中表达，表达的阈值为5 TPM）。

精子发生过程中的GEP是一组表达模式相似的基因构成的基因集^［9］。因此，为了鉴定精子发生过程中CTG所包含的GEP，对CTG在睾丸单细胞的表达谱进行聚类分析。聚类结果可以反映CTG在精子发生过程中的表达模式。

睾丸单细胞转录组数据来自2018—2019年发表的3篇睾丸单细胞转录组研究文献，其数据在GEO数据库中的标识号分别为GSE109037^［10］、GSE112013^［9］和GSE124263^［11］。这3套睾丸的单细胞转录组均是10X平台产生的数据。由于选用的参考基因组的版本不同，导致3套数据在基因名称等方面不兼容，不利于后续分析。因此，从原始的测序reads出发，重新处理这3套数据的原始数据，比对到统一的参考基因组上，有利于整合这3组数据。

将下载的原始测序reads转换成FASTQ格式，并将文件名命名为cellranger可以接受的格式，最后使用cellranger对基因的表达水平进行定量。下游分析使用Python语言的scanpy库进行^［12］。定量后的数据，经过质控、过滤、合并，使用scvi-tools进行去批次效应和过滤双细胞^［13-14］。随后，使用leiden图聚类算法对细胞进行聚类^［15］。leiden算法中影响聚类的参数主要是resolution。为了确定最优的聚类参数，在一系列resolution下使用leiden对细胞进行聚类，并用轮廓系数（silhouette score）评估聚类的效果。Silhouette score是一个经典的用于评估聚类结果的算法，数值越大说明聚类的结果越好。Silhouette score是每个样本单独计算的，由2个部分组成：a表示样本与同一类中其他点之间的平均距离，b为样本与次近类中其他点之间的平均距离。单个样本的silhouette score（s），可以表示为式4：

s的值域为 s∈-1,+1，-1表示该样本分类错误，+1表示该样本分类良好，s位于0附近表示该样本的分类不明确。为了衡量数据集（K个样本）的整体聚类情况，使用数据集中每个样本的silhouette score的平均值，记为S，表示为式5：

S 值越大，表示数据集中整体样本的分类情况越好。选择一系列resolution参数，根据每个参数的聚类结果计算出 S。最终选择 S 最大时对应的resolution作为聚类的参数，并以此参数确定细胞类群。随后，通过查阅文献获取之前研究^{［9-11，16-17］}中睾丸组织不同细胞类型的标志基因（表1），并根据这些基因在不同细胞类群中的表达情况对所有细胞类群进行注释。

为了鉴定精子发生过程中CTG包含的GEP，首先查看所有CTG在睾丸单细胞转录组中的表达情况。由于CTG主要表达于生殖系细胞（见结果部分），随后选择生殖系细胞的表达谱用于鉴定CTG包含的GEP，具体步骤如下：①提取CTG在生殖系细胞中的表达矩阵。②使用leiden图聚类算法对CTG进行聚类。③类似于细胞聚类，选定一系列resolution，并用silhouette score评估聚类结果，选择最终的resolution。④使用最优的聚类参数对CTG进行聚类，得到表达模式相似的CTG组成的基因集。⑤使用热图展示CTG的表达情况，观察不同基因集表达模式的异同。

这些基因集包含不同的基因表达模式，每个基因集表示一个GEP，并且这些GEP之间互不重叠。每个GEP表示为G，其包含的基因数为N。在GEP内部，以每个基因与该基因集中心的斯皮尔曼相关系数表示该基因的表达模式与其所属GEP的相似程度，表示为 ρg，则 ρg可以表示为式6：

式6中，spearmanr表示斯皮尔曼相关系数的计算函数，由SciPy库提供。xg表示GEP中的一个基因在所有细胞中的表达量组成的向量。xg¯ 表示基因集的中心，定义为式7：

为了定量表示GEP在每个细胞（scRNA-seq）或者样本（Bulk RNA-seq）中的活跃程度，使用DecoupleR包中的wmean函数计算每个GEP的活跃程度^［18］，表示为 aG。通过wmean先计算出GEP内基因表达量的权重平均值，再根据权重平均值的零分布标准化前面计算出来的权重平均值。标准化后的权重平均值定义为GEP的活跃程度。其中权重为前文中定义的基因表达模式与其所属GEP的相似程度，即 ρg。GEP内基因表达量的权重平均值定义为式8，yg 表示一个基因在一个细胞/样本中的表达量：

对一个细胞中所有基因，按照每个GEP对应的基因数重复抽样1 000次，并计算其权重平均值，得到该GEP权重平均值的零分布。零分布的均值和方差表示为 μG 和 σG2。由此，GEP在一个细胞/样本中的活跃程度可以表示为式9：

生存分析用于探索GEP的活跃程度与患者预后的关系，使用R语言中的survival包实现。首先，通过DecoupleR包的wmean函数计算TCGA中每例患者GEP的活跃程度；再将患者分为GEP活跃程度高组和GEP活跃程度低组，比较2组生存曲线；同时，基于患者的 GEP 活跃程度与患者的生存情况，进行Cox回归分析。

共筛选出1 271个TSG，进一步筛选出917个CTG。其中，鉴定到的CTG与WANG等^［4］鉴定到的CTG有705个重合，其中收录于CTdatabase数据库^［3］的有117个。新鉴定到212个CTG，其中收录于CTdatabase数据库的有16个（图1A）。

图1B比较了CTG和非CTG在睾丸组织中的表达情况，结果表明CTG整体上的表达水平高于非CTG的表达水平。图1C、D统计了CTG在染色体上的分布情况，其中CTG在X染色体上显著富集。图1C表明CTG在X染色体上的占比远大于其他染色体。通过Fisher's检验计算CTG在每条染色体上的富集情况（图1D），结果显示CTG在X染色体上的比值比（odds ratio，OR）高于其他染色体，且比值比显著大于1（P=0.000），说明CTG显著富集于X染色体上。

为了得到睾丸单细胞转录图谱，统一定量3套数据的原始测序数据，并进行质控和过滤。采用scvi-tools的整合算法将3组数据整合，共得到39 601个细胞。对整合后的睾丸组织单细胞转录组进行细胞聚类分析（图2）。结果显示，在一系列resolution参数中，当silhouette score取最大值时，对应的resolution为0.8（图2A）。基于此resolution值将细胞聚类为22个类。图2C、D分别展示了睾丸组织中体细胞和生殖系细胞的标志基因表达情况。基于此标注每个聚类类别的细胞类型，最终得到睾丸组织的单细胞转录组图谱（图2B）。共有13种细胞类型，每种细胞类型的全称、数量和标志基因详见表1。

CTG主要在生殖系细胞中表达，而在睾丸体细胞中不表达或低表达（图3A）。因此，CTG中包含的GEP主要是精子发生过程中的GEP。根据GEP鉴定的方法，在不同的resolution参数下，计算CTG聚类后的silhouette score。结果表明，在resolution = 0.5时，silhouette score达到最大值（图3B）。因此，以此参数对CTG进行聚类，将其分成7类，即7个GEP。每个GEP偏好表达于特定的细胞类型和精子发生时期（图3C）。

从染色体分布来看，GEP5包含的CTG绝大部分位于X染色体上（图3D）。这表明文献中提到的CT-X基因，在GEP的层面上主要集中在GEP5中。此外，图3C表明GEP5包含的基因在精子发生过程的早期阶段表达水平较高。

对7个GEP在每个细胞中的活跃程度进行定量计算，并用UMAP降维分析展示每个GEP在每个细胞中的活跃程度（图4A、B）。结果显示：GEP0主要活跃于Early primary SPC、Late primary SPC、Round ST中；GEP1几乎在每个细胞类型中都是活跃的，但在SSC、Differentiating SPG、Early primary SPC中的活跃程度较低；GEP2几乎在每个细胞类型中都是不活跃的；GEP3和GEP4主要活跃于Round ST、Elongating ST和Elongated ST中，GEP3在Round ST中最活跃，GEP4在Elongating ST中最活跃；GEP5主要活跃于SSC、Differentiating SPG、Early primary SPC；GEP6仅在Round ST中活跃。细胞类型缩写及含义见表1。

生存分析的结果汇总如图4C所示，整体来看，GEP5的活跃程度与患者的生存期有一定的关系。在大部分的肿瘤类型中，GEP5活跃程度高与患者的生存差有显著的关联，其中ACC（adrenocortical carcinoma）、BRCA（breast invasive carcinoma）、COAD（colon adenocarcinoma）、KIRC（kidney renal clear cell carcinoma）、KIRP（kidney renal papillary cell carcinoma）和UCEC（uterine corpus endometrial carcinoma）最显著（至少有4个风险比显著大于1）。BRCA、KIRC、KIRP、UCEC的疾病特异性生存（disease specific survival，DSS）曲线如图5所示，GEP5活跃程度高组的疾病特异性生存期与GEP5活跃程度低组比较，显著缩短。以上这些数据说明GEP5在多种肿瘤类型中的活跃程度与患者的生存期呈负相关。

CTG具有局限表达于睾丸且在肿瘤中高表达的特征，可作为肿瘤的治疗靶点或诊断标志物。目前，很多这方面的研究聚焦于单基因分析，缺乏对基因间协同作用的研究。本研究利用生物信息学计算分析技术，探索CTG的协同表达规律。

本研究通过基因表达特异性分析，筛选出1 271个TSG，结合TCGA数据库进一步确认917个CTG。在睾丸组织中，CTG的表达水平整体上高于其他基因的表达水平。收录在CTdatabase而没有鉴定到的CTG有110个，CTdatabase根据CTG的表达特征将其分为睾丸限定的（testis-restricted）、睾丸/大脑限定的（testis/brain-restricted）、睾丸选择性的（testis-selective）3类^［3］。我们推测，这110个CTG在本研究分析的睾丸组织表达特异性较差，故造成此差异。由于GTEx数据库和TCGA数据库更新，收录了更多的数据，本研究用这2个数据库的所有编码基因作为初始候选基因，而CTdatabase的初始候选基因是文献中记录的CTG^［3］。因此，本研究筛选到的CTG远多于CTdatabase中收录的CTG。

为了鉴定精子发生过程中CTG包含的GEP，本研究重新整理了睾丸单细胞转录组数据，并标注了细胞类型。通过睾丸单细胞转录组数据的分析，发现CTG主要表达于生殖系细胞中。睾丸的体细胞类型有SC、LC、PC、SMC、EC、TM。除了SC之外，其他的细胞类型或多或少都会在其他成体组织或器官中存在，故特异表达于睾丸的CTG不会表达于这些细胞中。

根据睾丸生殖系细胞的基因表达谱，对CTG进行聚类，鉴定出7个精子发生过程中的GEP。每个GEP有其偏好表达的细胞类型和精子发生时期。对GEP的活跃程度定量计算后，发现GEP5活跃于SSC、Differentiating SPG、Early primary SPC细胞。此外，研究还发现GEP5中包含了大量的CT-X基因^［5-6］。这说明在精子发生过程的前期，大部分的CT-X基因需要活跃的转录。X染色体上的基因受到减数分裂性染色体失活调控（meiotic sex chromosome inactivation，MSCI），其在减数分裂过程的前期活跃转录，在完成染色体联会后转录活性快速受到抑制，并形成XY体^［19］。这或许意味着在精子发生过程中，GEP5的活跃程度受到MSCI的调控。

GEP5活跃的细胞类型是SSC、Differentiating SPG、Early primary SPC，即精原细胞增殖期和减数分裂早期。之前的研究结果表明，减数分裂早期相关的CTG对于肺腺癌的发生和发展有重要作用^［4］。因此，本研究将GEP5在肿瘤中的活跃情况与患者的预后进行关联分析，以探究GEP5的活跃程度与患者预后之间的关系。结果显示，GEP5在ACC、BRCA、COAD、KIRC、KIRP和UCEC中的活跃程度高与患者的生存状况差有显著的关联。这说明GEP5在多种肿瘤类型中的活跃程度与患者的预后有关，具有成为多种肿瘤类型预后判断标志物的潜力。

综上所述，本研究通过多数据库联合分析，对CTG进行重新筛选和分类，并鉴定出7个GEP。其中，GEP5活跃于精子发生过程的前期，且富集CT-X基因。此外，在肿瘤中，GEP5的活跃程度越高，患者的生存状况越差，有望支持多种肿瘤类型的预后判断。