骨髓样本的转录组分析和生物信息学分析，以识别急性骨髓性白血病中的化疗耐药性特征

摘要

该方案提供了一套标准化的生物信息学工作流程，用于分析急性髓性白血病（AML）中的转录组变化。目标是比较新诊断和复发的骨髓样本，并优先考虑与化疗耐药性和疾病进展相关的分子特征，以便后续调查。

摘要

急性骨髓性白血病（AML）是一种高度异质性的血液恶性肿瘤，复发和获得性化疗耐药仍是治疗失败的主要原因。本文提出了一种用于骨髓抽吸转录组分析的生物信息学方案。该方案的主要目标是提供标准化工作流程，以识别与复发急性髓性白血病（AML）疾病进展和治疗耐药性的分子特征。该流程详细介绍了比较未配对骨髓样本的计算程序，通过五例新诊断病例和四例复发病例的测序数据进行了演示。该方法概述了处理RNA测序数据、进行差异基因表达分析以及进行下游功能评估的基本步骤。应用该工作流程识别出2025个差异表达基因（DEG），包括FOXC1、HOXA11、HOXA11-AS和AXL，作为与复发相关的候选转录本。功能和网络分析进一步优先考虑与小GTP酶信号、炎症信号、细胞外基质相互作用及RNA生物合成过程相关的基因集和交互中心。总体而言，该方法为复发AML相关的转录组特征绘制图谱及产生需要后续实验验证的假设提供了可重复的计算流程。

引言

急性髓系白血病（AML）是一类起源于造血干细胞/祖细胞的克隆恶性肿瘤，其特征是骨髓中未成熟骨髓系细胞异常增殖，造血分化受抑制，^1,2,3。尽管目前标准的诱导化疗（如西他菌碱与菽环素联合使用）可诱导大多数患者完全缓解，但复发率仍高达50%–70%，复发患者的预后显著受损，为^4,5,6。获得性化疗耐药性与AML治疗失败相关，因此需要深入分析与该过程相关的分子特征，以提升治疗策略和患者^生存率7,8。

在更广泛的文献中，现有研究表明，AML中的化疗耐药性不仅限于药物排流泵的上调或异常代谢，还与白血病干细胞（LSCs）的存活维持、血液生态位内上皮-间充质转变（EMT）样表型的形成以及与骨髓^微环境的相互作用有关^9，10,11。例如，LSC群体表现出较高的自我更新能力和静止性，这与对细胞周期特异化疗剂的固有耐药性有关¹²。此外，受体酪氨酸激酶（如AXL）的上调与FLT3-ITD+ AML的耐药性相关，同时激活PI3K/AKT和MAPK通路，增强抗凋亡能力^13,14。代谢重编程和表观遗传重塑也被认识到是抗药形成中重要的调控轴。证据表明，AML细胞在复发期间可能通过增强氧化磷酸化（OXPHOS）活性、调节的NAD/NADH比率以及改变组蛋白修饰来适应化疗诱导的氧化应激和DNA损伤，15,16,17。骨髓微环境中的炎症因子，如IL-6和CXCL8，也与LSC存活和化疗耐药性相关，通常与STAT3/NF-κB信号通路的激活同步发生（^11,18）。

尽管这些机制已被认可，但与AML复发和化疗耐药相关的转录组变化仍未完全被表征，尤其是在临床样本中同时评估代谢、表观遗传和骨髓微环境相关信号时。该工作流程通过优先考虑与初诊到临床复发过渡相关的候选DEGs、通路和调控网络来满足这一需求。该方法整合了差异表达分析、基因集富集分析（GSEA）和蛋白质-蛋白质相互作用（PPI）网络构建，以绘制系统范围的转录重编程图谱，并推荐候选方案进行后续机制验证。

该方法的总体目标是提出一个标准化、可重复的生物信息学流程，用于比较新诊断与复发AML骨髓样本的整体转录组。采用这种计算机模拟技术的理由在于其能够捕捉无偏的全基因组转录事件，超越单一途径分析的局限，优先关注复杂多维调控网络。该技术相比微阵列或靶向多路复用qPCR面板等其他方法具有显著优势，提供更高的动态范围、检测新转录本的能力以及基因表达的精确定量，且不受预设探针限制^19,20。为判断该方法是否适合其应用，读者应注意该流程专为研究者处理配对或未配对临床队列（如组织抽吸）的大规模RNA测序数据设计。它适合识别广泛的抗药性相关特征和候选调控网络，而需要细胞类型特异性或空间分辨率的研究人员则需要使用互补的单细胞或空间测序工作流程。最终，该计算协议使候选基因、通路和调控网络能够优先排序，以便后续的实验研究。

方案

所有涉及人体组织采样的方法均遵循机构指南和赫尔辛基宣言（2013年修订）进行。临床骨髓样本经机构伦理委员会批准（批准号TY-ZKY2024-116-01和TY-ZKY2024-116-02）。

1. 临床样本采集与患者分类

1. 根据世界卫生组织（WHO）分类标准，从正式诊断为急性骨髓性白血病（AML）患者的骨髓抽吸样本中选取。
2. 在患者选择时应用具体的纳入和排除标准，以确保队列的同质性和可重复性。纳入原发性AML的成人患者，排除继发性AML、急性骨髓细胞白血病或既往有其他恶性肿瘤病史的患者（见表1）。
3. 如果患者在初次临床诊断时出现此前未治疗的AML，则应将采集的样本分配给新诊断的组别。
4. 如果患者在完全缓解后外周血中白血细胞复发，或骨髓细胞发生超过5%的白血病性细胞，则将采集的样本归入复发组。
5. 在常规临床骨髓抽吸手术后，立即采集去识别的残留骨髓样本。
   注：在该特定方案的演示中，从2024年9月到2025年9月采集了九个连续样本（五个新诊断，四个复发）。由于该项目仅使用去标识的残留临床样本，伦理委员会免除了书面知情同意的要求。
6. 立即处理采集的骨髓样本，使用标准的酚-鸟嘌呤裂^解法保存RNA。
   1. 将新鲜采集的骨髓抽吸转移到含有抗凝剂的收集管中。摇晃管子，彻底混合吸吸剂和抗凝剂。
   2. 提取一定体积的抗凝样品，并将其加入商业化的酚-镍裂解试剂中。保持3份裂解试剂对1份样品的体积比。
      注意：酚-鸟苷溶解试剂含有有毒和腐蚀性化学物质，可能导致严重烧伤和组织损伤。在化学熏罩内进行所有试剂操作，同时佩戴适当的个人防护装备。
   3. 用力摇晃管子，使样品和裂解试剂完全均质化。确保混合物完全混合，并确认溶液中没有可见的血块残留。
   4. 通过将管子浸泡在液氮中，立即对均质样品进行速冻。
      注意：液氮极冷，接触时可能引发严重冻伤。处理液氮时请戴低温手套和全罩式防护罩。
7. 将快冻样本转移到-80°C的冷冻库进行长期保存，然后进行下游RNA分离和转录组测序流程。这代表了一个安全点，可以暂停实验，之后再重新开始。
   注意：本协议所呈现的工作流程完全专注于生成计算阻力签名。未对通过该特定管道识别的关键差异表达基因进行独立实验验证，如实时定量PCR（RT-qPCR）。

2. RNA质量控制与文库制备

1. 使用微流控毛细管电泳系统评估RNA完整性。在这一代表性工作流程中，包含RNA完整性数（RIN）≥6.0，A260/280比值介于1.8至2.1之间，且无可见降解峰的RNA样本。在准备文库前，记录每个样品的测得的RIN和纯度比。
2. 每个样品输入1微克总RNA用于文库制备。利用多寡糖附着磁珠从总RNA中纯化mRNA，以富集多A尾转录本。
3. 利用二价阳离子对富集mRNA进行片段处理。将混合物在94°C下培养15分钟，置于5X第一链合成反应缓冲液中。
4. 使用随机六聚体引物和一种缺乏RNase H活性的逆转录酶合成第一股cDNA。
5. 利用RNase H降解RNA模板链。利用DNA聚合酶I和dNTP在20微升反应系统中合成第二股cDNA。
6. 在16°C下孵育第二链合成反应1小时。在2000 x g 下短暂离心反应混合物，收集管底液体。
7. 利用外切酶和聚合酶活性将剩余的悬垂转化为钝端。对DNA片段的3'端进行腺苷酸化处理，并用发夹环结构连接接合器，为杂交做准备。
8. 使用磁性固相可逆固定珠纯化文库片段，优先选择长度为370–420碱基对的cDNA片段。
9. 在珠子净化过程中进行乙醇洗涤。在2000 x g 下离心管30秒，收集并彻底去除残留乙醇，然后再进行最终洗脱。
10. 使用高精度DNA聚合酶、通用PCR引物和样本特异性指数引物进行PCR扩增。
11. 在98°C下进行初始变性，执行PCR热剖面，持续30秒。接着进行12个循环，分别是98°C持续10秒、60°C30秒和72°C持续30秒，最后延长为72°C持续5分钟。
12. 再次用磁珠纯化PCR产物。应用第2.9步相同的离心参数，获得最终的文库。
13. 用荧光仪定量初始文库浓度。将最终文库稀释至1.5 ng/μL浓度。
14. 把稀释后的库彻底搅拌均匀。在4°C下，将混合物以10,000 x g 离心1分钟，以去除残留杂质，然后进行最终分析。
15. 利用微流控毛细电泳系统评估插入件库的尺寸。
16. 确认插入件尺寸符合预期后，通过实时定量PCR（qRT-PCR）准确定量文库的有效浓度。确保浓度高于1.5 nM，以确保文库稳定性和测序质量。
    注意：这是一个可以暂停实验的安全点。制备的文库可在-20°C下保存，直到聚类和测序。

3. 聚类与转录组测序

1. 在自动聚类生成系统上对索引编码样本进行聚类。按照厂家说明使用商用的成对端集群套件。
2. 在成功生成簇后，在高通量测序平台上对文库制备进行测序。生成150个碱基对（bp）对端读段。

4. 数据质量控制与读取映射

1. 使用fastp v0.23.2评估原始数据（FASTQ格式）的质量，进行原始读取质量控制和过滤。将命令行参数记录在分析日志中。在该工作流程中，通过使用相同的样本过滤设置，移除含适配器的读段、包含多重N序列的读段和低质量读段，生成了干净的读段。补充 文件1中提供了代表性的对对端命令。
2. 原始读数通过自动预处理软件处理。通过去除含有适配器的读段、包含多元素N序列的读段和低质量读段，获得干净的读段。对所有样品使用相同的滤波参数，并在滤波后记录保留的读数、Q20、Q30和气谱含量。
3. 计算干净数据中的Q20、Q30和GC含量。在下游分析前定义潜在的批次变量，包括样本采集日期、RNA提取日期、文库准备批次、测序通道和测序运行。
4. 通过PCA和样本间相关分析，使用归一化表达值评估批量效应。如果样本主要按技术变量而非临床状态聚集，记录受影响变量并将其作为协变量纳入差分表达设计公式，或在下游可视化前应用已建立的批量调整方法。
5. 获取参考基因组（智人GRCh38）及相应的Ensembl 109版本基因注释文件以进行读序列比对。
6. 使用HISAT2 v2.0.5构建参考基因组索引。
7. 使用HISAT2 v2.0.5将配对端的干净读段与参考基因组比对。使用这种基于基因模型注释文件的剪接比对方法，生成一个连接连接数据库。

5. 新颖转录本预测与基因表达定量

1. 使用StringTie v1.3.3b以基于参考的方法组装每个样本的映射读段。利用该工具组装和定量代表每个基因位点多个剪接变异的完整转录本。
2. 使用featureCounts v1.5.0-p3统计每个基因被映射的读取次数。使用所得的原始整数读计数矩阵作为下游微分表达式分析的输入。
3. 使用配对端选项（例如 -p）配置 featureCounts v1.5.0-p3 以实现配对端测序数据。提供下载的GRCh38 GTF注释文件以定义正确的基因组特征边界。
4. 计算每个基因每百万个映射读段的转录本片段数（FPKM）。仅在描述性可视化、主成分分析（PCA）、热图显示和探索性表达式摘要中使用FPKM值;不要将FPKM值作为DESeq2差分表达式测试的输入矩阵。

6. 差异基因表达分析

1. 使用R v3.5.0和DESeq2 R 套件v1.20.0，对新诊断组和复发组进行差异表达分析。将步骤5.2生成的原始读数矩阵导入R环境，并仅保留FPKM值用于可视化和探索性分析。
2. 构建分析软件包所需的专用数据集对象。执行特定命令（例如DESeqDataSetFromMatrix（））将计数数据矩阵与对应的样本元数据表绑定。
3. 在软件对象中定义实验设计公式。指定临床状态（新诊断与复发）作为比较的主要变量（例如，设计 = ~ 病况）。如果在第4.3步中识别出技术性批次变量且与临床状态不完全混淆，应将其纳入设计公式（例如，设计 = ~ 批次 + 条件）。
4. 执行核心差分表达式分析函数（例如，DESeq（））。允许软件自动执行尺寸因子估计、色散估计和负二项Wald检验拟合²²。
5. 利用结果提取函数（例如，results（））提取结果表。具体说明对比论证以定义具体对比（复发与新诊断）。
6. 调整所得的P值以控制假发现率。利用软件包²³自动应用的Benjamini和Hochberg集成程序。
7. 过滤提取的结果表，以分离显著的差异表达基因（DEGs）。将任何调整后的P值<0.05且绝对log₂ 倍变化>1的基因，都指定为显著表达差异。

7. 基因本体（GO）富集分析

1. 使用clusterProfiler v3.8.1和org.Hs.eg.db v3.6.0。输入对应第6.7步中识别出的重要DEG的Entrez基因ID列表。
2. 执行GO富集函数（例如enrichGO（））。指定所需参数，包括相应的背景生物数据库（例如，OrgDb = org。Hs.eg.db）、特定本体域（生物过程、细胞成分或分子功能）以及调整后的P值截止值（0.05）。
3. 确保算法在浓缩计算过程中应用必要的修正。确认软件内部是否能校正基因长度偏差，并使用Benjamini和Hochberg方法调整P值²⁴。
4. 将校正P值小于0.05的Go项视为显著富集。利用软件包集成的可视化功能生成点阵图或条形图，显示顶部丰富的 GO 术语。

8. 京都基因与基因组百科全书（KEGG）途径富集分析

1. 利用一个专门理解高层次生物系统功能的综合数据库资源，以识别失调的通路。准备与第7.1步中使用的相同重要DEG入口ID清单。
2. 在函数注释R包中执行KEGG富集函数（例如enrichKEGG（））。
3. 在函数调用中定义关键参数。将生物体代码严格设置为人类（例如，生物体 = 'has'），并定义P值调整方法（例如，pAdjustMethod = 'BH'）。
4. 提取具有统计意义的KEGG通路。过滤输出，只保留修正P值小于0.05的通路。
5. 想象富集最丰富的KEGG通路。利用集成绘图函数（例如点图）来绘制每条通路相关的统计显著性和基因计数。

9. 基因集富集分析（GSEA）

1. 准备分析所需的预先排序基因列表。计算所有表达基因的排名指标，使用带符号的-log10（P值）乘以差异表达分析得出的log₂ 倍变化符号。
2. 启动Broad Institute GSEA软件的本地安装，v4.2.3。将新生成的预排序基因列表输入软件界面²⁵。
3. 下载所需的预定义基因集。获取分子签名数据库（MSigDB，版本7.5.1）²⁶中的基因本体（GO）和京都基因与基因组百科全书（KEGG）数据集。
4. 配置软件参数以执行统计富集测试。将排列数设置为1000，并选择排列类型为“gene_set”。
5. 执行分析算法，确定预定义的基因组是否在新诊断与复发的生物状态之间表现出统计显著且一致的差异。
6. 评估生成的富集曲线的统计显著性。用严格阈值定义显著基因集：归一化富集（NES）绝对值>1.0，名义P值<0.05，假发现率（FDR）q值<0.25。

10. 蛋白质间相互作用（PPI）网络分析

1. 访问STRING数据库，了解已知和预测的蛋白质-蛋白质相互作用。在此工作流程中，PPI分析使用STRING v11.5²⁷进行。
2. 将Entrez基因ID或显著差异表达基因（在第6.7步识别）的官方基因符号列表输入数据库搜索界面。选择智人作为目标生物。
3. 配置网络构建参数，确保获取高质量交互。将最低要求交互得分设为高置信阈值（得分>0.700）。
4. 将生成的交互网络数据导出到本地目录。将交互映射保存为标准表格文件（例如TSV格式）。
5. 将导出的交互数据导入 Cytoscape v3.9.1 进行网络可视化和分析²⁸.
6. 过滤构建网络以提升可视化清晰度，突出关键监管枢纽。去除任何未连续相互作用达到既定置信阈的断开节点或孤儿基因。

结果

临床队列与测序验证

通过测序产量和质量指标，确认了上游RNA提取和文库制备协议（图1）的成功执行。在该具代表性数据集中，来自五名新诊断AML患者和四名复发AML患者的骨髓样本，平均每个样本获得了约6.0GB的原始数据。质量控制评估（表2）确认碱基质量和读段深度符合下游生物信息分析所需的阈值9。低RNA完整性（例如RIN < 6.0）、低定位率或高转录本降解偏倚将代表输入质量不佳，并降低下游差异表达分析的可靠性。

全局转录组方差与PCA

为评估全球转录组方差并检查临床分组情况，对归一化表达数据进行了PCA。在该代表性数据集中，新诊断组和复发组在二维空间中表现出分离（见图2A）^20，PC1和PC2分别占总方差的23.82%和18.75%。图2B，C中的维恩图提供了新诊断和复发组样本中检测到的基因的额外描述性总结，支持在下游差异表达分析前进行样本层级的重复性检查。由于队列规模小且未配对，PCA分离被视为一个说明性工作流输出，而非疾病状态特异生物学的确凿证据。

差异表达基因（DEG）分析

将既定的协议阈值（|log₂FC| ≥ 1，调整后的P值≤0.05）应用于DESeq2输出，发现复发组中有2,025个DEG基因，包括772个上调基因和1,253个下调基因（见图3A）。变异较大的候选转录包括FOXC1（log₂FC = 7.55，P = 4.92 x ^10-5）、HOXA11（log₂FC = 7.76）、HOXA11-AS（log₂FC = 7.23）和AXL（log₂FC = 3.50），以及下调的RHOB、PTX3和CXCL8。现有文献将这些基因中的若干与AML干性、信号传导或治疗反应联系起来^13,29;然而，目前的工作流程仅将其识别为与复发相关的候选转录本。任何明确的临床耐药机制作用都需要后续独立的功能验证。

功能与途径丰富（通路、KEGG和GSEA）

功能注释协议将DEGs映射到更广泛的生物系统。抗性酶分析发现了与小GTP酶介导信号转导、金属离子转运和染色质组装相关术语的富集现象（见图4A–C）。KEGG通路定位发现了与ECM受体相互作用及细胞因子-细胞因子受体相互作用的关联（见图4D）。GSEA显示复发组RNA生物合成过程富集，新诊断组能量代谢通路富集（见图5A）。这些富集结果提供了改变基因集的描述性路线图，应被解读为假设生成的关联，而非被证实的复发驱动因素。

蛋白质间相互作用（PPI）网络构建

最初的STRING网络包含56个节点和193个交互。在移除断开或孤立节点后，显示的Cytoscape子网络包含42个节点和136个交互（见图5B）。网络模块分析优先将TP53、CCL2、CXCL8和IL6作为交互次数最多的中心数学枢纽。由于PPI网络依赖数据库预测的相互作用评分（例如ATF3评分：0.982），枢纽识别应被视为未来实证研究的目标优先级，而非p53介导的凋亡规避或其他耐药机制的直接证据。

该协议生成的原始RNA测序数据已存入Figshare存储库，并通过以下DOI公开访问：https://doi.org/10.6084/m9.figshare.30655814。处理过的数据及相关分析文件包含在文章及其补充材料中。用于重现计算工作流程的代表性命令行参数和分析设置作为 补充文件1提供了。支持本研究发现的所有数据均可自由获取。

患者识别	年龄（年龄）	性	分子突变	存活/随访（数月）	临床状况
R_AML_1	70	男性	FLT3-ITD （+）	22	已故
R_AML_2	29	女性	NPM1 （+）	11	还活着
R_AML_3	40	男性	CEBPA （+）	17	还活着
R_AML_4	55	女性	三重否定*	24	已故

表1：复发急性髓性白血病（R_AML）组患者的人口统计和临床特征。 表1总结了代表性分析中使用的复发AML队列的人口统计和临床特征，包括与转录组工作流程解释相关的患者层面临床特征。

样本	图书馆	Raw_reads	Raw_bases	Clean_reads	Clean_bases	Error_rate	Q20	Q30	GC_pct
AML_1	FRAS25 0244891-1r	48705066	7.31G	47807532	7.17G	0.01	99.35	97.48	47.48
AML_2	FRAS25 0244896-1r	42969940	6.45G	42237962	6.34G	0.01	99.35	97.44	46.74
AML_3	FRAS2502 44906-1r	48738386	7.31G	47744462	7.16G	0.01	99.36	97.48	47.28
AML_4	FRAS250 244915-1r	48723650	7.31G	47688240	7.15G	0.01	99.29	97.26	47.45
AML_5	FRAS2502 44920-1r	49508198	7.43G	47740308	7.16G	0.01	99.37	97.53	47.73
R_AML_1	FRAS2502 44892-1r	47879408	7.18G	46671584	7.0G	0.01	99.39	97.49	47.63
R_AML_2	FRAS2502 70005-1R	47657378	7.15G	46957882	7.04G	0.01	99.39	97.49	50.5
R_AML_3	FRAS250 405722-1R	58754766	8.81G	56867112	8.53G	0.01	99.38	97.42	46.52
R_AML_4	FRAS2502 44902-1R	48491122	7.27G	47469334	7.12G	0.01	99.23	97.21	46.43

表2：数据质量总结。 表2报告了每个样品的测序质量指标，包括读取率、碱性质量、GC含量以及用于判断样品是否适合下游分析的定位相关质量控制信息。

图1：协议工作流程。 该工作流程总结了主要的实验和计算阶段，包括临床样本采集、RNA质量控制、文库制备与测序、读取处理与比对、转录本定量、差分表达分析、 GO/KEGG 富集、GSEA以及PPI网络构建。 请点击此处查看该图的放大版本。

图2：样本定量分析。 （A） 进行了主成分分析（PCA）以评估组间差异及组内样本的可重复性。主体分析采用基于所有样本归一化基因表达值的线性代数方法进行。（B，C）维恩图分别显示了AML和R_AML组样本中检测到的基因。样本限制区域表示在单个样本中检测到的基因，而重叠区域则表示在两个或多个样本中常见的基因。 请点击此处查看该图的放大版本。

图3：差异基因表达分析。 （A） 条形图显示了对照组间差异表达基因（DEGs）数量，由DESeq2识别，调整后P值阈值≤0.05，|log₂折叠变化|≥ 1.（B）DEGs火山图。横轴表示 log₂FoldChange 值，纵轴表示 -log10（P 值）。蓝色虚线表示用于DEG选择的阈值线。（C） DEGs的分层聚类热图。横轴表示样本名称，纵轴表示DEG的归一化表达式值。 请点击此处查看该图的放大版本。

图4：差异表达基因的功能富集分析。 （A） GO 丰富酒吧地块。横轴表示 GO 项，纵轴表示富集显著性，表示为 -log10（padj）。颜色代表BP（生物过程）、CC（细胞成分）和MF（分子功能）。（B）GO浓缩气泡图。横轴表示每个GO项注释DEG与总DEG数量的比例，y轴表示GO项。气泡大小对应注释基因的数量，颜色渐变代表富集的重要性。（C）KEGG丰富酒吧地块。横轴表示KEGG通路，纵轴表示富集重要性。（D）KEGG富化气泡图。气泡大小表示注释基因数量，颜色渐变反映富集显著性。 请点击此处查看该图的放大版本。

图5：GSEA富集与蛋白质-蛋白质相互作用（PPI）网络分析。 （A） 显示选定显著基因集的规范化富集（NES）条形图。NES值为正表示R_AML组富集，NES值为负表示新诊断AML组富集。（B） 蛋白质-蛋白质相互作用（PPI）网络。每个节点代表一个蛋白质，每条边表示连接蛋白质之间的相互作用。 请点击此处查看该图的放大版本。

讨论

协议的关键步骤

成功执行这一生物信息学工作流程依赖于若干关键步骤。首先，骨髓抽吸物的即时快速冷冻和正确溶解（步骤1.6）至关重要，因为骨髓微环境中富含核糖核酸酶，能迅速降解转录组完整性³⁰。在计算阶段，DESeq2包（步骤6.3）中实验设计公式的正确配置对于准确的差异表达至关重要，尤其是在临床状态（新诊断与复发）进行对比并考虑潜在混杂变量时。最后，在基因集富集分析（GSEA）（步骤9.6）中应用严格的伪发现率（FDR）阈值，是防止对假阳性功能网络过度解读的关键统计检查点。

改装与故障排除

该方法的一个常见挑战是批次效应的存在，这在临床样本收集和测序时，在较长时间内频繁出现。分析前应定义批次变量，包括样本采集日期、RNA提取日期、文库制备批次、测序通道和测序运行。如果PCA或样本相关分析发现基于测序日期或其他技术变量而非临床表型的聚类，用户应在统计可行的情况下将批量变量纳入差分表达设计公式，或在可视化前应用批量修正算法，如ComBat或SVA，修改方案³¹.如果将该方案应用于全血而非骨髓抽取，一个关键改进是在文库制备过程中加入一个globin mRNA耗尽步骤，以防止大量globin转录本垄断测序读段深度。软件版本和代表性工作流程的主要参数补充如下：fastp v0.23.2、HISAT2 v2.0.5、StringTie v1.3.3b、featureCounts v1.5.0-p3、R v3.5.0、DESeq2 v1.20.0、clusterProfiler v3.8.1、org。Hs.eg.db v3.6.0，配对端150碱基对测序，GSEA v4.2.3，包含1000个基因集排列，MSigDB v7.5.1，STRING v11.5，高置信度交互，以及Cytoscape v3.9.1。代表性的命令行参数和分析设置见 补充文件1。

方法的局限性

虽然该协议内容全面，但存在固有的方法论局限性。首先，它采用了体型RNA测序，能够捕捉整个骨髓抽吸物的平均转录组结构，但缺乏单细胞空间分辨率。因此，工作流程无法确定上调的复发相关特征是否源自白血病干细胞、基质细胞、免疫细胞，或细胞类型成分的变化³²。其次，代表性数据集体积小（n = 9）且未配对，限制了统计鲁棒性，阻碍了确定的因果推断。第三，工作流程纯粹是计算机模拟。它生成候选调控枢纽和信号通路，但若无正交体外或体内实验验证，无法独立验证其在化疗抵抗中的功能必要性。

近期的单细胞和单细胞基因组学研究通过解析细胞状态异质性、克隆结构及治疗相关进化，扩展了AML参考框架，分辨率^{更高，分辨率为33,34,35,36}。这些方法与本文所述的体RNA测序工作流程相辅相成：体型测序为队列级转录组签名提供了实用且经济的筛选策略，而单细胞和多组学方法则可用于后续研究，将候选信号分配给特定的恶性或微环境细胞群体。

相对于现有方法的重要性

尽管存在这些局限性，该转录组管线相较于其他诊断和分析技术具有优势。传统的AML复发临床评估通常依赖靶向多重qPCR面板或标准流式细胞术。虽然这些针对性方法适合快速诊断，但受预定义探针限制，只能评估已知电阻标记¹⁹。通过结合无偏的全基因组转录组测序与网络分析，该方案能够推荐现有靶向方法可能忽视的新转录本和系统范围关联。

重要性与潜在应用

本方案中提出的方法论对转化血液学和个体化医学具有重要意义，因为它可以优先考虑与复发相关的转录组特征，供进一步研究。一个潜在的后续应用是选择复发期间出现的表面抗原或免疫逃避通路。这些候选药物可能为未来验证研究设计提供参考，如果实验确认，可能有助于下一代免疫疗法的发展，包括CAR-T或CAR-NK细胞策略³⁷。这些转化应用应被视为假设生成，而非从现有数据集中确立的结论。

材料

本文使用的材料清单

姓名	公司	目录编号	评论
Agilent 2100 Bioanalyzer	Agilent Technologies, Santa Clara, CA, USA	RRID:SCR_019389	G2939BA
AMPure XP system	Beckman Coulter, Brea, CA, USA	RRID:SCR_008452	A63881
cBot Cluster Generation System	Illumina, San Diego, CA, USA	-	SY-301-2002 or institution-specific system ID
clusterProfiler (Software)	Bioconductor	RRID:SCR_016884	v3.8.1
Cytoscape (Software)	Cytoscape Consortium	RRID:SCR_003032	v3.9.1
DESeq2 (Software)	Bioconductor	RRID:SCR_015687	v1.20.0
DNA Polymerase I	New England Biolabs (NEB, Ipswich, MA, USA	-	M0209L
dNTP Solution Mix	New England Biolabs (NEB, Ipswich, MA, USA	-	N0447L
edgeR (Software)	Bioconductor	RRID:SCR_012802	v3.22.5
fastp (Software)	OpenGene	RRID:SCR_016962	v0.23.2
featureCounts / Subread (Software)	The Walter and Eliza Hall Institute	RRID:SCR_012919	featureCounts v1.5.0-p3
GRCh38 reference genome	Genome Reference Consortium / Ensembl	-	GRCh38; Ensembl release 109
GSEA software	Broad Institute	RRID:SCR_003199	v4.2.3
HISAT2 (Software)	Johns Hopkins University	RRID:SCR_015530	v2.0.5
M-MuLV Reverse Transcriptase (RNase H-)	New England Biolabs (NEB, Ipswich, MA, USA	-	M0253L
MSigDB gene sets	Broad Institute	RRID:SCR_016863	v7.5.1
NEBNext Ultra II Directional RNA Library Prep Kit for Illumina	New England Biolabs (NEB, Ipswich, MA, USA	-	E7760L/E7765L or laboratory-specific kit
NEBNext Ultra II RNA Library Prep Kit for Illumina	New England Biolabs (NEB, Ipswich, MA, USA	-	E7770L
NovaSeq sequencing platform	Illumina, San Diego, CA, USA	RRID:SCR_016387	NovaSeq system; service-provider instrument ID
org.Hs.eg.db (Annotation package)	Bioconductor	-	v3.6.0
Phusion High-Fidelity DNA Polymerase	Thermo Fisher Scientific, Waltham, MA, USA	RRID:AB_2756816	F530L
Qubit 2.0 Fluorometer	Thermo Fisher Scientific, Waltham, MA, USA	RRID:SCR_018095	Q32866
Qubit dsDNA HS Assay Kit	Thermo Fisher Scientific, Waltham, MA, USA	-	Q32851
R software	R Foundation for Statistical Computing	RRID:SCR_001905	v3.5.0
Random Hexamer Primer	Thermo Fisher Scientific, Waltham, MA, USA	-	SO142
RNA 6000 Nano Kit	Agilent Technologies, Santa Clara, CA, USA	-	5067-1511
RNase H	New England Biolabs (NEB, Ipswich, MA, USA	-	M0297L
RNA-seq Library Prep Kit / Sequencing Service	Novogene, Beijing, China	-	Project No. X101SC25054246-Z01-J003
STRING database	STRING Consortium	RRID:SCR_005223	v11.5
StringTie (Software)	Johns Hopkins University / Center for Computational Biology	RRID:SCR_016323	v1.3.3b
TruSeq PE Cluster Kit v3-cBot-HS	Illumina, San Diego, CA, USA	-	PE-401-3001