方法文章

急性髓系白血病中化疗耐药特征的骨髓样本转录组学分析与生物信息学研究

DOI：

10.3791/70750

⸱

2026年8月4日

 ,  ,  ,  , 
1Department of Hematology, the Affiliated Ganzhou Hospital of Nanchang University, 2Department of Hematology, Affiliated Nanshan Hospital of Shenzhen University

* These authors contributed equally

本文内容

摘要

$$\rightleftharpoonup{xx}$$ $$\longleftharp{xx}$$, $$\longrightharp{xx}$$,

本方案介绍了一种标准化的生物信息学工作流程，用于分析急性髓系白血病（AML）中的转录组改变。目标是比较初诊和复发性骨髓样本，筛选出与化疗耐药及疾病进展相关的分子特征，以供后续研究优先开展。

摘要

$$\rightleftharpoonup{xx}$$ $$\longleftharp{xx}$$, $$\longrightharp{xx}$$,

急性髓系白血病（AML）是一种高度异质性的血液系统恶性肿瘤，其中复发和获得性化疗耐药仍是治疗失败的主要原因。本文介绍了一种针对骨髓穿刺液的转录组学生物信息学分析流程。该流程的主要目标是提供一个标准化的工作流程，用于识别与复发 AML 疾病进展和治疗耐药相关的分子特征。该分析流程详细描述了对非配对骨髓样本进行比较的计算步骤，示例数据来自 5 例初诊患者和 4 例复发患者的测序数据。该方法概述了处理 RNA 测序数据、进行差异基因表达分析以及开展下游功能评估的关键步骤。应用此工作流程在代表性数据集中鉴定出 2,025 个差异表达基因（DEGs），包括 FOXC1、HOXA11、HOXA11-AS 和 AXL，这些基因被列为与复发相关的候选转录本。功能分析和网络分析进一步筛选出与小 GTP 酶信号通路、炎症信号通路、细胞外基质相互作用以及 RNA 生物合成过程相关的关键基因集合和相互作用枢纽。总体而言，该方法提供了一个可重复的计算分析流程，可用于绘制与复发 AML 相关的转录组特征图谱，并生成需后续实验验证的科学假说。

引言

$$\rightleftharpoonup{xx}$$ $$\longleftharp{xx}$$, $$\longrightharp{xx}$$,

急性髓系白血病（AML）是一组起源于造血干细胞/祖细胞的克隆性恶性肿瘤，其特征是骨髓中未成熟髓系细胞异常增殖，并抑制造血分化1,2,3。尽管当前的标准诱导化疗方案（如阿糖胞苷联合蒽环类药物）可使大多数患者达到完全缓解，但复发率仍高达50%–70%，且复发患者的预后显著恶化4,5,6。获得性化疗耐药与AML治疗失败密切相关，因此有必要深入分析与此过程相关的分子特征，以改进治疗策略并提高患者生存率7,8。

在更广泛的研究文献中，已有研究表明，急性髓系白血病（AML）的化疗耐药性不仅限于药物外排泵的上调或药物代谢异常，还与白血病干细胞（LSCs）的存活维持、血液学微环境内类上皮-间质转化（EMT）表型的形成以及与骨髓微环境的相互作用相关9,10,11。例如，LSC群体表现出高度的自我更新能力和静息状态，这使其对细胞周期特异性化疗药物具有内在耐药性12。此外，受体酪氨酸激酶（如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）复发和化疗耐药相关的转录组变化仍未能被充分阐明，尤其是在临床样本中同时评估代谢、表观遗传以及骨髓微环境相关信号时。本工作流程通过优先筛选与从初诊到临床复发转变过程相关的差异表达基因（DEG）、通路及调控网络，弥补了这一不足。该方法整合了差异表达分析、基因集富集分析（GSEA）以及蛋白质-蛋白质相互作用（PPI）网络构建，以系统描绘全转录组重编程，并筛选出可用于后续机制验证的候选分子。

本方法的总体目标是提供一个标准化、可重复的生物信息学流程，用于比较新诊断与复发性急性髓系白血病（AML）骨髓样本的整体转录组。采用这种计算机模拟技术的依据在于其能够无偏倚地捕获全基因组范围的转录事件，突破单一通路分析的局限，从而优先识别复杂的多维调控网络。与微阵列或靶向多重qPCR面板等替代方法相比，该技术具有显著优势，包括更高的动态范围、检测新转录本的能力，以及在无需预设计探针限制的情况下对基因表达进行精确定量19,20。为判断该方法是否适用于其研究，读者应注意，此流程专为处理来自配对或非配对临床队列（如组织穿刺样本）的批量RNA测序数据的研究人员设计。该方法适用于识别广泛的耐药相关特征及候选调控网络，而若研究人员需要细胞类型特异性或空间分辨率，则需采用互补的单细胞或空间测序工作流程。最终，这一计算流程可帮助优先筛选出候选基因、通路和调控网络，供后续实验研究进一步验证。

方案

$$\rightleftharpoonup{xx}$$ $$\longleftharp{xx}$$, $$\longrightharp{xx}$$,

所有涉及人体组织取样的方法均遵循机构指南和《赫尔辛基宣言》（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. 采用标准的苯酚-胍裂解法立即处理采集的骨髓穿刺液，以保存RNA21。
    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 µg总RNA用于文库构建。利用连接有poly-T寡核苷酸的磁珠从总RNA中纯化mRNA，以富集带有polyA尾的转录本。
  3. 使用二价阳离子对富集的mRNA进行片段化。在5X第一链合成反应缓冲液中，于94 °C孵育15分钟。
  4. 使用随机六聚体引物和缺乏RNase H活性的逆转录酶合成第一链cDNA。
  5. 使用RNase H降解RNA模板链。在20 µL反应体系中，利用DNA聚合酶I和dNTPs合成第二链cDNA。
  6. 将第二链合成反应在16 °C孵育1小时。短暂离心（2,000 x g）反应混合物，使液体聚集于管底。
  7. 利用外切酶和聚合酶活性将剩余的突出末端转化为平末端。对DNA片段的3'端进行腺苷酸化，并连接带有发夹环结构的接头，以准备后续杂交。
  8. 使用磁性固相可逆固定化珠纯化文库片段，优先选择长度为370–420 bp的cDNA片段。
  9. 在珠纯化过程中进行乙醇洗涤。离心管子（2,000 x g，30秒），以收集并彻底去除最终洗脱前的任何残留乙醇。
  10. 使用高保真DNA聚合酶、通用PCR引物和样本特异性索引引物进行PCR扩增。
  11. 执行PCR热循环程序：98 °C初始变性30秒，随后进行12个循环（98 °C变性10秒，60 °C退火30秒，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 格式）进行原始读段质量控制与过滤，评估数据质量。将命令行参数记录在分析日志中。在本工作流程中，通过去除含有接头的读段、含有 poly-N 序列的读段以及低质量读段，采用统一的过滤设置在所有样本中生成干净读段（clean reads）。代表性双端测序命令见补充文件 1。
  2. 通过自动化预处理软件处理原始读段。通过去除含有接头的读段、含有 poly-N 序列的读段以及低质量读段，获得干净读段。对所有样本使用相同的过滤参数，并记录过滤后保留的读段数量、Q20、Q30 和 GC 含量。
  3. 计算干净数据的 Q20、Q30 和 GC 含量。在进行下游分析前，定义潜在的批次变量，包括样本采集日期、RNA 提取日期、文库制备批次、测序通道（lane）和测序运行批次（sequencing run）。
  4. 利用标准化表达值，通过主成分分析（PCA）和样本间相关性分析评估批次效应。若样本主要按技术性变量而非临床状态聚类，需记录受影响的变量，并将其作为协变量纳入差异表达分析的设计公式中，或在下游可视化前应用已确立的批次校正方法。
  5. 获取参考基因组（Homo sapiens, GRCh38）及对应的 Ensembl release 109 基因注释文件，用于读段比对。
  6. 使用 HISAT2 v2.0.5 构建参考基因组的索引。
  7. 使用 HISAT2 v2.0.5 将双端干净读段比对至参考基因组。采用此剪接位点感知（splice-aware）的比对方法，基于基因模型注释文件生成剪接位点数据库。

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

  1. 采用基于参考基因组的方法，使用 StringTie v1.3.3b 对每个样本的比对读段进行转录本组装。利用该工具组装并定量每个基因位点的全长转录本，以表征多种剪接变异体。
  2. 使用 featureCounts v1.5.0-p3 统计比对到每个基因的读段数量。将得到的原始整数读段计数矩阵作为下游差异表达分析的输入数据。
  3. 针对成对末端测序数据，使用成对末端选项（例如 -p）配置 featureCounts v1.5.0-p3。提供下载的 GRCh38 GTF 注释文件，以定义正确的基因组特征边界。
  4. 计算每个基因的每千碱基转录本每百万比对读段中的片段数（FPKM）。FPKM 值仅用于描述性可视化、主成分分析（PCA）、热图展示及探索性表达水平汇总；不得将 FPKM 值作为 DESeq2 差异表达检验的输入矩阵。

6. 差异基因表达分析

  1. 使用 R v3.5.0 和 DESeq2 R 软件包 v1.20.0 对新诊断组与复发组进行差异表达分析。将步骤 5.2 中生成的原始读段计数矩阵导入 R 环境，并仅保留 FPKM 值用于可视化和探索性分析。
  2. 构建分析软件包所需的专用数据集对象。执行特定命令（例如 DESeqDataSetFromMatrix()），将计数数据矩阵与相应的样本元数据表进行绑定。
  3. 在软件对象中定义实验设计公式。将临床状态（新诊断与复发）指定为主要比较变量（例如，design = ~ condition）。如果在步骤 4.3 中识别出技术批次变量，且该变量与临床状态未完全混杂，则将其纳入设计公式中（例如，design = ~ batch + condition）。
  4. 执行核心差异表达分析函数（例如 DESeq()）。允许软件自动进行规模因子估计、离散度估计以及负二项分布 Wald 检验拟合22。
  5. 利用结果提取函数（例如 results()）提取结果表格。通过指定 contrast 参数来定义确切的比较内容（复发组 versus 新诊断组）。
  6. 对得到的 P 值进行校正，以控制错误发现率。使用软件包自动应用的 Benjamini 和 Hochberg 综合校正程序23。
  7. 对提取的结果表格进行筛选，以鉴定出显著差异表达的基因（DEGs）。将任何校正后 P 值 < 0.05 且绝对 log2 倍数变化 > 1 的基因定义为显著差异表达基因。

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

  1. 使用 clusterProfiler v3.8.1 和 org.Hs.eg.db v3.6.0 对鉴定出的差异表达基因（DEGs）进行基因本体（GO）富集分析。输入对应于步骤 6.7 中鉴定出的显著差异表达基因的 Entrez 基因 ID 列表。
  2. 执行 GO 富集分析函数（例如 enrichGO()），指定必要的参数，包括合适的背景物种数据库（例如 OrgDb = org.Hs.eg.db）、特定的本体域（生物过程、细胞组分或分子功能）以及校正后的 P 值阈值（0.05）。
  3. 确保算法在富集计算过程中应用必要的校正。确认软件在内部对基因长度偏差进行校正，并采用 Benjamini 和 Hochberg 方法对 P 值进行校正24。
  4. 将校正后 P 值小于 0.05 的 GO 条目视为显著富集。利用该软件包内置的可视化功能生成点图或柱状图，以展示最显著富集的 GO 条目。

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

  1. 利用专注于解析高级生物系统功能的综合性数据库资源，识别失调的通路。准备与步骤 7.1 中所用相同的显著差异表达基因（DEG）Entrez ID 列表。
  2. 在功能注释 R 软件包中执行 KEGG 富集分析功能（例如 enrichKEGG()）。
  3. 在函数调用中定义关键参数。将物种代码严格设置为人类（例如，organism = 'hsa'），并指定 P 值校正方法（例如，pAdjustMethod = 'BH'）。
  4. 提取具有统计学显著性的 KEGG 通路。对结果进行筛选，仅保留校正后 P 值小于 0.05 的通路。
  5. 可视化富集程度最高的 KEGG 通路。利用内置绘图函数（例如 dotplot()）展示各通路相关的统计显著性及基因数量。

9. 基因集富集分析 (GSEA)

  1. 准备分析所需的预排序基因列表。使用差异表达分析得到的 log2 倍数变化的符号乘以带符号的 -log10(P 值) 来计算所有表达基因的排序指标。
  2. 启动本地安装的 Broad Institute GSEA 软件 v4.2.3。将新生成的预排序基因列表输入软件界面25。
  3. 下载所需的预定义基因集。从分子特征数据库（MSigDB，版本 7.5.1）获取基因本体（GO）和京都基因与基因组百科全书（KEGG）数据集26。
  4. 配置软件参数以执行统计富集检验。将置换次数设置为 1,000，并选择置换类型为“gene_set”。
  5. 运行分析算法，以确定预定义的基因集在新诊断和复发的生物学状态之间是否表现出具有统计学意义的一致性差异。
  6. 评估生成的富集谱的统计学显著性。使用严格的阈值定义显著基因集：标准化富集评分（NES）绝对值 > 1.0，名义 P 值 < 0.05，以及错误发现率（FDR）q 值 < 0.25。

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

  1. 访问 STRING 数据库以获取已知和预测的蛋白质-蛋白质相互作用。在此工作流程中，PPI 分析使用 STRING v11.527 进行。
  2. 将显著差异表达基因（在步骤 6.7 中鉴定出）的 Entrez 基因 ID 列表或官方基因符号输入数据库搜索界面。选择 Homo sapiens 作为目标生物体。
  3. 配置网络构建参数，以确保获取高质量的相互作用。将最低必需相互作用评分设置为高置信度阈值（评分 > 0.700）。
  4. 将生成的相互作用网络数据导出至本地目录。将相互作用图谱保存为标准表格文件（例如 TSV 格式）。
  5. 将导出的相互作用数据导入 Cytoscape v3.9.1，用于网络可视化与分析28。
  6. 对构建的网络进行过滤，以提高可视化清晰度并突出关键调控枢纽。移除任何未达到既定置信度阈值的连续相互作用的孤立节点或孤儿基因。

结果

$$\rightleftharpoonup{xx}$$ $$\longleftharp{xx}$$, $$\longrightharp{xx}$$,

临床队列与测序验证

上游RNA提取与文库构建方案的成功实施（图1通过测序产量和质量指标确认。在此代表性数据集中，来自五名新诊断急性髓系白血病（AML）患者和四名复发性AML患者的骨髓样本，每样本平均产生约6.0 GB的原始数据。质量控制评估（表 2) 确认碱基质量和测序深度达到下游生物信息学分析所需的阈值9。RNA完整性较低（例如，RIN < 6.0）、比对率低或转录本降解偏差高均表明输入材料质量不理想，会影响下游差异表达分析的可靠性。

全转录组变异性与主成分分析

为评估全局转录组变异并检验临床分组情况，对标准化表达数据进行了主成分分析（PCA）。在此代表性数据集中，新诊断组与复发组在二维空间中呈现出分离趋势（图2A)20，其中PC1和PC2分别占总方差的23.82%和18.75%。文氏图中的 图2B,C 为新诊断组和复发组中各样本检测到的基因提供额外的描述性总结，以支持在进行下游差异表达分析前的样本水平重复性检验。由于队列规模较小且样本未配对，主成分分析（PCA）的分离结果被解读为流程示例性输出，而非疾病状态特异性生物学的确定性证据。

差异表达基因（DEG）分析

将已建立的协议阈值（|log2FC| ≥ 1 且校正后 P 值 ≤ 0.05）应用于 DESeq2 的输出结果，共鉴定出 2,025 个差异表达基因（DEGs），其中包括在复发组中上调的 772 个基因和下调的 1,253 个基因（图 3A）。变异程度较高的候选转录本包括 FOXC1（log2FC = 7.55，P = 4.92 × 10-5）、HOXA11（log2FC = 7.76）、HOXA11-AS（log2FC = 7.23）和 AXL（log2FC = 3.50），以及下调的 RHOB、PTX3 和 CXCL8。已有文献表明其中多个基因与急性髓系白血病（AML）的干性、信号通路或治疗反应相关13,29；然而，本研究的工作流程仅将其识别为与复发相关的候选转录本。这些基因在临床耐药中是否具有明确的机制性作用，仍需后续独立的功能验证。

功能与通路富集分析（GO、KEGG 和 GSEA）

功能注释流程将差异表达基因（DEGs）映射到更广泛的生物系统中。基因本体（GO）分析发现，小GTP酶介导的信号转导、金属离子转运以及染色质组装相关的术语显著富集（图4A–C）。KEGG通路分析揭示了细胞外基质-受体相互作用和细胞因子-细胞因子受体相互作用的关联（图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_170男性FLT3-ITD (+)22已故
R_AML_229女性NPM1 (+)11存活
R_AML_340男性CEBPA (+)17存活
R_AML_455女性三重阴性*24已故

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

样本文库原始读段数原始碱基数干净读段数干净碱基数错误率Q20Q30GC含量百分比
AML_1FRAS25
0244891-1r		487050667.31G478075327.17G0.0199.3597.4847.48
AML_2FRAS25
0244896-1r		429699406.45G422379626.34G0.0199.3597.4446.74
AML_3FRAS2502
44906-1r		487383867.31G477444627.16G0.0199.3697.4847.28
AML_4FRAS250
244915-1r		487236507.31G476882407.15G0.0199.2997.2647.45
AML_5FRAS2502
44920-1r		495081987.43G477403087.16G0.0199.3797.5347.73
R_AML_1FRAS2502
44892-1r		478794087.18G466715847.0G0.0199.3997.4947.63
R_AML_2FRAS2502
70005-1r		476573787.15G469578827.04G0.0199.3997.4950.5
R_AML_3FRAS250
405722-1r		587547668.81G568671128.53G0.0199.3897.4246.52
R_AML_4FRAS2502
44902-1r		484911227.27G474693347.12G0.0199.2397.2146.43

表2：数据质量汇总。 表2报告了每个样本的测序质量指标，包括读段产量、碱基质量、GC含量以及用于判断样本是否适用于下游分析的比对相关质控信息。

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

figure-results-2
图 2：样本的定量分析。（A）进行主成分分析（PCA）以评估组间差异和组内样本的重复性。PCA基于所有样本的标准化基因表达值，采用线性代数方法进行。（B, C）维恩图显示在AML和R_AML组样本中检测到的基因。样本特有区域表示仅在单个样本中检测到的基因，而重叠区域表示在两个或多个样本中共同检测到的基因。请点击此处查看该图的放大版本。

figure-results-3
图3：差异基因表达分析。（A）柱状图显示通过DESeq2鉴定出的各比较组间差异表达基因（DEGs）的数量，筛选阈值为校正后P值 ≤ 0.05 且 |log2倍数变化| ≥ 1。（B）DEGs的火山图。x轴表示log2倍数变化值，y轴表示-log10(P值)。蓝色虚线表示用于筛选DEGs的阈值线。（C）DEGs的层次聚类热图。x轴表示样本名称，y轴表示DEGs的标准化表达值。请点击此处查看该图的放大版本。

figure-results-4
图4：差异表达基因的功能富集分析。（A）GO富集条形图。横轴表示GO术语，纵轴表示富集显著性，以-log10(padj)表示。颜色代表BP（生物过程）、CC（细胞组分）和MF（分子功能）。（B）GO富集气泡图。横轴表示注释到每个GO术语的差异表达基因（DEGs）占DEGs总数的比例，纵轴表示GO术语。气泡大小对应于注释基因的数量，颜色梯度表示富集显著性。（C）KEGG富集条形图。横轴表示KEGG通路，纵轴表示富集显著性。（D）KEGG富集气泡图。气泡大小表示注释基因的数量，颜色梯度反映富集显著性。 请点击此处查看该图的放大版本。

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

讨论

$$\rightleftharpoonup{xx}$$ $$\longleftharp{xx}$$, $$\longrightharp{xx}$$,

实验方案中的关键步骤

成功执行该生物信息学工作流程依赖于多个关键步骤。首先，骨髓穿刺液的立即速冻和适当裂解（步骤 1.6）至关重要，因为骨髓微环境富含核糖核酸酶，可迅速降解转录组完整性30。在计算分析阶段，正确设置 DESeq2 软件包中的实验设计公式（步骤 6.3）对于实现准确的差异表达分析尤为关键，尤其是在比较临床状态（初诊 versus 复发）的同时考虑潜在混杂变量的情况下。最后，在基因集富集分析（GSEA）过程中应用严格的错误发现率（FDR）阈值（步骤 9.6）是一个关键的统计学质控点，用以避免对假阳性功能网络的过度解读。

修改与故障排除

该方法中一个常见的挑战是批次效应的存在，当临床样本在较长时间内分批收集和测序时，常会出现此类问题。应在分析前定义批次变量，包括样本采集日期、RNA提取日期、文库制备批次、测序通道和测序运行批次。如果主成分分析（PCA）或样本相关性分析显示样本根据测序日期或其他技术变量聚类，而非基于临床表型，则用户应在统计可行的情况下，在差异表达分析的设计公式中纳入批次变量，或在可视化前应用批次校正算法（如 ComBat 或 SVA）进行调整31。若将本方案应用于全血样本而非骨髓穿刺样本，一个关键的修改是在文库制备过程中增加珠蛋白 mRNA 耗竭步骤，以防止高丰度的珠蛋白转录本占据过多的测序读长深度。代表性工作流程所使用的软件版本及主要参数补充如下：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 bp 测序、GSEA v4.2.3（采用 1,000 次基因集置换检验）、MSigDB v7.5.1、STRING v11.5（使用高置信度相互作用）、Cytoscape v3.9.1。代表性命令行参数和分析设置详见补充文件1。

方法的局限性

尽管该方案具有全面性，但仍存在固有的方法学局限性。首先，该方案采用批量RNA测序技术，仅能获取整个骨髓穿刺样本的平均转录组谱，缺乏单细胞水平的空间分辨率。因此，该工作流程无法确定上调的复发相关特征是来源于白血病干细胞、基质细胞、免疫细胞，还是细胞类型组成的变化32。其次，代表性数据集样本量较小（n = 9）且为非配对设计，限制了统计学的稳健性，并妨碍得出明确的因果推论。第三，该工作流程完全基于生物信息学分析，虽可识别出候选的调控枢纽和信号通路，但在缺乏正交的体外或体内实验验证的情况下，无法独立证实这些通路在化疗耐药中的功能必要性。

近期的单细胞及单细胞基因组学研究通过更高分辨率解析了急性髓系白血病（AML）的细胞状态异质性、克隆结构以及与治疗相关的演化过程，从而拓展了AML的参考研究框架33,34,35,36。这些方法与本文所述的批量RNA测序（bulk RNA-seq）工作流程相辅相成：批量测序为队列水平的转录组特征提供了实用且具有成本效益的筛查策略，而单细胞及多组学方法则可用于后续研究，将候选信号归因于特定的恶性细胞或微环境细胞群体。

相对于现有方法的意义

尽管存在这些局限性，该转录组学分析流程相较于其他诊断和分析技术仍具有优势。急性髓系白血病（AML）复发的传统临床评估通常依赖于靶向多重qPCR检测面板或常规流式细胞术。虽然这些靶向方法在快速诊断方面具有实用价值，但受限于预定义的探针，只能评估已知的耐药标志物19。通过采用无偏倚的全基因组转录组测序并结合网络分析，本方案能够识别出新颖的转录本以及系统范围内的关联，而这些可能是现有靶向方法所忽略的。

重要性与潜在应用

本方案中所述的方法学与转化血液学及个体化医学相关，因其可优先筛选出与复发相关的转录组特征，供进一步研究。一个潜在的后续应用是鉴定在复发过程中出现的表面抗原或免疫逃逸通路。此类候选分子可为未来验证性研究的设计提供依据，若经实验确认，可能有助于开发新一代免疫疗法，包括CAR-T或CAR-NK细胞治疗策略37。这些转化性应用应被视为基于当前数据集生成的假设，而非已确立的结论。

披露

$$\rightleftharpoonup{xx}$$ $$\longleftharp{xx}$$, $$\longrightharp{xx}$$,

作者声明无利益冲突。

致谢

$$\rightleftharpoonup{xx}$$ $$\longleftharp{xx}$$, $$\longrightharp{xx}$$,

本研究由赣州市科学技术局（2022—ZD1368）资助。

材料

本文使用的材料清单
姓名公司目录编号评论
Agilent 2100 生物分析仪Agilent Technologies, Santa Clara, CA, USARRID:SCR_019389G2939BA
AMPure XP 系统Beckman Coulter, Brea, CA, USARRID:SCR_008452A63881
cBot 簇生成系统Illumina, San Diego, CA, USA-SY-301-2002 或机构特定系统编号
clusterProfiler（软件）BioconductorRRID:SCR_016884v3.8.1
Cytoscape（软件）Cytoscape 联盟RRID:SCR_003032v3.9.1
DESeq2（软件）BioconductorRRID:SCR_015687v1.20.0
DNA 聚合酶 INew England Biolabs (NEB, Ipswich, MA, USA-M0209L
dNTP 溶液混合物New England Biolabs (NEB, Ipswich, MA, USA-N0447L
edgeR（软件）BioconductorRRID:SCR_012802v3.22.5
fastp（软件）OpenGeneRRID:SCR_016962v0.23.2
featureCounts / Subread（软件）沃尔特和伊丽莎·霍尔研究所RRID:SCR_012919featureCounts v1.5.0-p3
GRCh38 参考基因组基因组参考联盟 / Ensembl-GRCh38; Ensembl release 109
GSEA 软件Broad 研究所RRID:SCR_003199v4.2.3
HISAT2（软件）约翰斯·霍普金斯大学RRID:SCR_015530v2.0.5
M-MuLV 逆转录酶（RNase H-）New England Biolabs (NEB, Ipswich, MA, USA-M0253L
MSigDB 基因集Broad 研究所RRID:SCR_016863v7.5.1
NEBNext Ultra II 定向 RNA 文库构建试剂盒（适用于 Illumina）New England Biolabs (NEB, Ipswich, MA, USA-E7760L/E7765L 或实验室特定试剂盒
NEBNext Ultra II RNA 文库构建试剂盒（适用于 Illumina）New England Biolabs (NEB, Ipswich, MA, USA-E7770L
NovaSeq 测序平台Illumina, San Diego, CA, USARRID:SCR_016387NovaSeq 系统；服务提供商仪器编号
org.Hs.eg.db（注释包）Bioconductor-v3.6.0
Phusion 高保真 DNA 聚合酶Thermo Fisher Scientific, Waltham, MA, USARRID:AB_2756816F530L
Qubit 2.0 荧光仪Thermo Fisher Scientific, Waltham, MA, USARRID:SCR_018095Q32866
Qubit dsDNA HS 检测试剂盒Thermo Fisher Scientific, Waltham, MA, USA-Q32851
R 软件统计计算基金会 RRRID:SCR_001905v3.5.0
随机六聚体引物Thermo Fisher Scientific, Waltham, MA, USA-SO142
RNA 6000 Nano 试剂盒Agilent Technologies, Santa Clara, CA, USA-5067-1511
RNase HNew England Biolabs (NEB, Ipswich, MA, USA-M0297L
RNA-seq 文库构建试剂盒 / 测序服务诺禾致源，中国北京-项目编号 X101SC25054246-Z01-J003
STRING 数据库STRING 联盟RRID:SCR_005223v11.5
StringTie（软件）约翰斯·霍普金斯大学 / 计算生物学中心RRID:SCR_016323v1.3.3b
TruSeq PE 簇生成试剂盒 v3-cBot-HSIllumina, San Diego, CA, USA-PE-401-3001

参考文献

$$\rightleftharpoonup{xx}$$ $$\longleftharp{xx}$$, $$\longrightharp{xx}$$,
  1. Döhner H, et al. Diagnosis and management of AML in adults: 2017 ELN recommendations from an international expert panel. Blood. 2017;129(4):424-47.
  2. Papaemmanuil E, et al. Genomic classification and prognosis in acute myeloid leukemia. N Engl J Med. 2016;374(23):2209-21.
  3. Arber DA, et al. The 2016 revision to the World Health Organization classification of myeloid neoplasms and acute leukemia. Blood. 2016;127(20):2391-405.
  4. Burnett AK, et al. Optimization of chemotherapy for younger patients with acute myeloid leukemia: results of the Medical Research Council AML15 trial. J Clin Oncol. 2013;31(27):3360-8.
  5. Ravandi F, et al. Phase I/II study of combination therapy with sorafenib, idarubicin, and cytarabine in younger patients with acute myeloid leukemia. J Clin Oncol. 2010;28(11):1856-62.
  6. Thol F, Schlenk RF, Heuser M, Ganser A. How I treat refractory and early relapsed acute myeloid leukemia. Blood. 2015;126(3):319-27.
  7. Zahreddine H, Borden KL. Mechanisms and insights into drug resistance in cancer. Nat Rev Drug Discov. 2013;12(4):314-28.
  8. Meads MB, Gatenby RA, Dalton WS. Environment-mediated drug resistance: a major contributor to minimal residual disease. Nat Rev Cancer. 2009;9(9):665-74.
  9. Ishikawa F, et al. Chemotherapy-resistant human AML stem cells home to and engraft within the bone-marrow endosteal region. Nat Biotechnol. 2007;25(11):1315-21.
  10. Pollyea DA, et al. Targeting acute myeloid leukemia stem cells: a review and principles for the development of clinical trials. Haematologica. 2014;99(8):1277-84.
  11. Shafat MS, Gnaneswaran B, Bowles KM, Rushworth SA. The bone marrow microenvironment-home of the leukemic blasts. Blood Rev. 2017;31(5):277-86.
  12. Lagadinou C, et al. BCL-2 inhibition targets oxidative phosphorylation and selectively eradicates quiescent human leukemia stem cells. Cell Stem Cell. 2013;12(3):329-41.
  13. Park IK, et al. Receptor tyrosine kinase Axl is required for resistance of leukemic cells to FLT3-targeted therapy in acute myeloid leukemia. Leukemia. 2015;29(12):2382-9.
  14. Ben-Batalla I, et al. Axl, a prognostic and therapeutic target in acute myeloid leukemia mediates paracrine crosstalk of leukemia cells with bone marrow stroma. Blood. 2013;122(14):2443-52.
  15. Farge T, et al. Chemotherapy-resistant human acute myeloid leukemia cells are not enriched for leukemic stem cells but require oxidative metabolism. Cancer Discov. 2017;7(7):716-35.
  16. Bosc C, et al. Mitochondrial inhibitors circumvent survival of leukemic stem cells. Cancer Cell. 2017;32(5):724-5.
  17. Li S, et al. Distinct evolution and dynamics of epigenetic and genetic heterogeneity in acute myeloid leukemia. Nat Med. 2016;22(7):792-9.
  18. Stevens AM, et al. Increased STAT3 activity in acute myeloid leukemia is driven by the bone marrow microenvironment. Cancer Res. 2017;77(18):5053-62.
  19. Wang Z, Gerstein M, Snyder M. RNA-Seq: a revolutionary tool for transcriptomics. Nat Rev Genet. 2009;10(1):57-63.
  20. Conesa A, et al. A survey of best practices for RNA-seq data analysis. Genome Biol. 2016;17(1):13.
  21. Chomczynski P, Sacchi N. Single-step method of RNA isolation by acid guanidinium thiocyanate-phenol-chloroform extraction. Anal Biochem. 1987;162(1):156-9.
  22. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.
  23. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B Stat Methodol. 1995;57(1):289-300.
  24. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284-7.
  25. Subramanian A, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545-50.
  26. Liberzon A, et al. Molecular signatures database (MSigDB) 3.0. Bioinformatics. 2011;27(12):1739-40.
  27. Szklarczyk D, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019;47(D1):D607-13.
  28. Shannon P, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498-504.
  29. Somerville TD, et al. Frequent derepression of the mesenchymal transcription factor gene FOXC1 in acute myeloid leukemia. Cancer Cell. 2015;28(3):329-42.
  30. Stark R, Grzelak M, Hadfield J. RNA sequencing: the teenage years. Nat Rev Genet. 2019;20(11):631-56.
  31. Leek JT, et al. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. 2012;28(6):882-3.
  32. Hwang B, Lee JH, Bang D. Single-cell RNA sequencing technologies and bioinformatics pipelines. Exp Mol Med. 2018;50(8):1-14.
  33. Döhner H, et al. Diagnosis and management of AML in adults: 2022 recommendations from an international expert panel on behalf of the ELN. Blood. 2022;140(12):1345-77.
  34. Khoury JD, et al. The 5th edition of the World Health Organization classification of haematolymphoid tumours: myeloid and histiocytic/dendritic neoplasms. Leukemia. 2022;36(7):1703-19.
  35. Morita K, et al. Clonal evolution of acute myeloid leukemia revealed by high-throughput single-cell genomics. Nat Commun. 2020;11:5327.
  36. Miles LA, et al. Single-cell mutation analysis of clonal evolution in myeloid malignancies. Nature. 2020;587(7834):477-82.
  37. Kenderian SS, et al. CD33-specific chimeric antigen receptor T cells exhibit potent preclinical activity against human acute myeloid leukemia. Leukemia. 2015;29(8):1637-47.

重印与许可

申请许可以重复使用本 JoVE 文章的文本或图表

申请许可

标签

234 234

相关文章