转录组学综合性使用范例

转录组学分析在实际研究中有着广泛的应用场景，从病原感染研究到癌症生物学，从环境适应到发育调控，不同的研究问题需要采用不同的分析策略和工具组合。本节通过五个具有代表性的实际案例，展示如何将前面介绍的转录组学分析方法整合应用于具体的科学问题。这些案例涵盖了差异表达分析、功能富集、可变剪接检测、单细胞转录组、从头组装和非编码RNA研究等多个方面，帮助你理解如何根据研究目标设计合适的分析流程。

范例一：病原感染后宿主转录组响应（差异表达 + 功能富集 + 可变剪接）

在病毒感染研究中，理解宿主细胞如何响应病原体入侵是开发抗病毒策略的基础。一个典型的应用场景是将人类细胞系暴露于流感病毒，在感染后6小时、12小时和24小时收集RNA进行双端链特异性RNA-seq测序。这种时间序列实验设计可以捕捉宿主抗病毒反应的动态变化过程，从早期信号传导到晚期效应分子表达的完整图景。

整个分析流程从原始数据的质量控制开始。使用fastp对双端测序数据进行接头修剪和质量过滤，去除低质量碱基和接头序列以获得干净的读段。随后使用STAR比对器将这些读段比对到人类参考基因组hg38上，STAR的优势在于它可以准确处理跨越剪接位点的读段，这对于后续的可变剪接分析尤为重要。比对完成后，使用featureCounts对每个基因进行定量统计，生成基因水平的表达计数矩阵用于下游的差异表达分析。

差异表达分析采用DESeq2方法，设计公式中包含时间点作为主要因素。通过对比24小时感染组与0小时未感染对照组，可以鉴定出约2000个显著差异表达的基因。这些基因代表了宿主细胞对流感病毒感染的全面转录响应，包括上调的抗病毒效应分子和下调的与正常细胞功能相关的基因。

为了理解这些差异基因的功能含义，使用clusterProfiler进行GO功能富集分析和KEGG通路富集分析。GO分析结果显示干扰素信号通路和炎症反应等免疫相关过程被显著激活，这符合宿主对抗病毒感染的预期反应。KEGG富集分析进一步揭示RIG-I-like受体通路是关键的抗病毒信号传导途径，该通路在识别病毒RNA并启动先天免疫反应中发挥核心作用。

除了基因表达水平的变化外，可变剪接也是病毒感染过程中重要的调控层面。使用rMATs工具可以系统检测感染过程中的差异剪接事件。在这个案例中检测到200个显著的差异剪接事件，其中特别值得关注的是ISG15基因的外显子跳跃事件。ISG15是一个重要的干扰素刺激基因，其蛋白产物具有广泛的抗病毒活性。外显子跳跃可能导致ISG15产生截短的蛋白异构体，从而削弱其正常的抗病毒功能，这可能是病毒用来逃避宿主防御的一种机制。

# 病原感染宿主转录组响应分析流程
# 1. 差异表达分析 (DESeq2)
library(DESeq2)

# 读取计数数据
count_data <- read.csv("gene_counts_matrix.csv", row.names = 1)
col_data <- data.frame(
    sample = colnames(count_data),
    condition = factor(c(rep("0h", 3), rep("6h", 3), rep("12h", 3), rep("24h", 3)))
)
rownames(col_data) <- col_data$sample

# 创建DESeqDataSet对象
dds <- DESeqDataSetFromMatrix(countData = count_data,
                              colData = col_data,
                              design = ~ condition)
# 预过滤低表达基因
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep, ]
# 运行DESeq2标准化和差异分析
dds <- DESeq(dds)
# 对比24h vs 0h
res <- results(dds, contrast = c("condition", "24h", "0h"))
# 筛选显著差异基因 (FDR < 0.05, |log2FC| > 1)
sig_genes <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
cat(sprintf("鉴定到 %d 个显著差异基因\n", nrow(sig_genes)))

# 2. 功能富集分析 (clusterProfiler)
library(clusterProfiler)
library(org.Hs.eg.db)

# 转换为Entrez ID
gene_list <- rownames(sig_genes)
entrez_ids <- bitr(gene_list, fromType = "SYMBOL",
                   toType = "ENTREZID", OrgDb = org.Hs.eg.db)

# GO富集分析
go_result <- enrichGO(gene = entrez_ids$ENTREZID,
                      OrgDb = org.Hs.eg.db,
                      ont = "BP",
                      pAdjustMethod = "BH",
                      pvalueCutoff = 0.05,
                      qvalueCutoff = 0.10)
head(go_result)

# KEGG通路富集
kegg_result <- enrichKEGG(gene = entrez_ids$ENTREZID,
                          organism = 'hsa',
                          pAdjustMethod = "BH",
                          pvalueCutoff = 0.05)
head(kegg_result)

# 3. 可变剪接分析 (rMATs)
# 构建rMATs输入文件列表
python rmats.py --b1 control.txt --b2 infected_24h.txt \
    --gtf /path/to/gencode.v44.annotation.gtf \
    -t paired --readLength 150 --nthread 12 \
    --od rmats_output --tmp tmp_dir

# 4. 可视化: 火山图、热图、富集气泡图
# 使用R的EnhancedVolcano, pheatmap, ggplot2包进行可视化

可视化部分对于结果的展示和解释至关重要。火山图可以直观地展示所有基因的表达变化幅度和统计显著性，帮助快速识别最显著的上下调基因。热图则能够显示差异基因在不同样本间的表达模式，揭示可能的共表达聚类。富集气泡图将GO或KEGG富集结果以图形化方式呈现，圆圈大小代表基因数量，颜色代表富集显著性。对于可变剪接事件，Sashimi plot是一种专门的可视化方法，它可以在基因组浏览器中展示特定基因在不同条件下的剪接模式差异，直观地显示外显子包含或跳跃的情况。

这个综合分析案例的意义在于它揭示了病毒与宿主之间复杂的互作机制。通过整合基因表达和可变剪接两个层面的信息，我们不仅可以发现哪些抗病毒基因被诱导表达，还可以了解病毒是否通过操纵宿主的剪接机制来逃避免疫攻击。ISG15的外显子跳跃事件就是一个很好的例子，它提示了潜在的新治疗靶点——如果能够阻止这种有害的剪接事件，可能增强宿主的抗病毒能力。这类多组学整合分析策略正在成为感染性疾病研究中的标准方法。

范例二：癌症亚型分类与生物标志物发现（单细胞转录组 + 差异表达）

肿瘤组织内部的细胞异质性是癌症研究和治疗面临的核心挑战之一。三阴性乳腺癌作为一种缺乏明确治疗靶点的侵袭性乳腺癌亚型，其肿瘤微环境中存在多种细胞类型的复杂相互作用。单细胞RNA-seq技术为解析这种异质性提供了前所未有的分辨率，使我们能够在单个细胞的水平上识别不同的细胞亚群、推断细胞状态转换轨迹并预测患者对治疗的反应。

单细胞RNA-seq数据的处理从Cell Ranger流程开始，这是10X Genomics平台配套的分析软件套件。Cell Ranger负责将原始的fastq文件比对到参考基因组、生成特征-条码矩阵并进行初步的细胞 calling。输出结果包括每个细胞的基因表达矩阵，其中行是基因，列是细胞条码，值是UMI计数。这个原始表达矩阵需要经过严格的质量控制才能用于下游分析。

质量控制是单细胞数据分析中最关键的一步之一，因为低质量细胞或空液滴会引入噪声并影响后续的所有分析结果。常用的过滤标准包括：每个细胞检测到的基因数应在合理范围内（通常200-8000个），过低的基因数可能是空液滴或死细胞，过高则可能是双细胞；线粒体基因 reads 的比例应低于10-20%，高比例表明细胞可能处于应激状态或已经死亡；每个细胞的总UMI数也应在合理范围内。这些阈值需要根据具体的数据集进行调整，可以通过小提琴图可视化各指标的分布来确定合适的 cutoff 值。

归一化和批次校正是处理多个样本时必不可少的步骤。sctransform是一种基于负二项广义线性模型的归一化方法，它相比传统的log normalization更能有效去除技术噪声并保留生物学变异。当实验涉及多个患者样本或不同批次时，还需要使用Harmony、Seurat的Integration功能或Liger等方法进行批次校正，以确保不同来源的细胞可以在同一个空间中进行比较和分析。

降维和聚类是发现细胞亚群的核心步骤。主成分分析（PCA）首先将高维的基因表达数据降维到主要的几个主成分上，这些主成分捕捉了数据中最主要的变异来源。然后使用UMAP或t-SNE等非线性降维方法进一步将数据投影到二维或三维空间中进行可视化。Louvain或Leiden算法则在构建的K近邻图上进行社区检测，将相似的细胞聚集成簇。在三阴性乳腺癌的这个案例中，聚类分析识别出了15个不同的细胞簇，每个簇代表一种特定的细胞类型或状态。

标记基因注释是将计算得到的细胞簇赋予生物学意义的步骤。通过对每个簇进行差异表达分析，找到在该簇中特异性高表达的标记基因，然后利用已知的细胞类型标记数据库来确定每个簇的身份。在这个案例中，恶性细胞簇高表达上皮细胞标志物EPCAM和细胞角蛋白基因如KRT8、KRT18；巨噬细胞簇表达CD68和CD163；T细胞簇表达CD3D、CD4或CD8A。这种注释结果可以帮助我们理解肿瘤微环境的组成和各细胞成分的比例。

拷贝数变异推断是区分恶性细胞和正常细胞的重要手段。InferCNV算法通过比较肿瘤细胞与正常参照细胞之间的基因表达模式来推断大规模的染色体拷贝数变化。恶性肿瘤细胞通常表现出明显的染色体扩增或缺失模式，而浸润的正常细胞则保持正常的二倍体状态。这种方法特别适用于没有已知肿瘤特异性标记基因的情况，可以从计算角度识别恶性细胞群体。

拟时序分析旨在重建细胞状态转变的动态过程。Monocle3、Slingshot或RNA velocity等方法可以根据细胞的基因表达相似性构建轨迹图，推断细胞从一个状态向另一个状态转变的路径和方向。在这个三阴性乳腺癌案例中，Monocle3分析预测部分恶性细胞正在经历上皮-间质转化过程，向间充质状态转变，这些过渡态细胞高表达VIM、ZEB1等间质标志物。上皮-间质转化与肿瘤转移和耐药性密切相关，因此识别处于这一转变过程中的细胞具有重要的临床意义。

细胞通讯分析揭示了肿瘤微环境中不同细胞类型之间的相互作用网络。CellPhoneDB是一种常用的细胞通讯推断工具，它基于已知配体-受体相互作用数据库，评估每种细胞类型组合之间是否存在显著的配体-受体信号交流。分析结果显示，肿瘤相关巨噬细胞通过分泌TGF-β等因子作用于恶性细胞，诱导其发生上皮-间质转化。这一发现揭示了肿瘤微环境中免疫抑制机制的一个具体环节，提示靶向巨噬细胞-TGF-β轴可能有助于阻断肿瘤的恶性进展。

# 单细胞转录组分析流程 (Seurat v5)
library(Seurat)
library(harmony)
library(monocle3)
library(InferCNV)

# 1. 创建Seurat对象并QC过滤
sc_data <- Read10X(data.dir = "filtered_feature_bc_matrix/")
seu <- CreateSeuratObject(counts = sc_data, project = "TNBC", min.cells = 3, min.features = 200)

# 计算线粒体基因比例
seu[["percent.mt"]] <- PercentageFeatureSet(seu, pattern = "^MT-")
# 过滤低质量细胞
seu <- subset(seu, subset = nFeature_RNA > 200 & nFeature_RNA < 7500 &
              percent.mt < 10 & nCount_Rna > 500)

# 2. sctransform归一化
seu <- SCTransform(seu, verbose = FALSE)

# 3. 批次校正 (Harmony)
seu <- RunPCA(seu, verbose = FALSE)
seu <- RunHarmony(seu, group.by.vars = "patient_id")

# 4. 降维与聚类
seu <- RunUMAP(seu, reduction = "harmony", dims = 1:30)
seu <- FindNeighbors(seu, reduction = "harmony", dims = 1:30)
seu <- FindClusters(seu, resolution = 0.8)  # 获得15个簇

# 5. 标记基因鉴定
markers <- FindAllMarkers(seu, only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0.25)
head(markers)

# 6. 细胞类型注释
cell_type_markers <- list(
    Epithelial = c("EPCAM", "KRT8", "KRT18"),
    Macrophage = c("CD68", "CD163"),
    T_cell = c("CD3D", "CD4", "CD8A")
)

这个综合性案例的意义在于它展示了如何利用单细胞转录组学技术深入解析肿瘤的复杂性。通过识别耐药或转移相关的细胞亚群（如正在经历上皮-间质转化的恶性细胞），我们可以更好地理解肿瘤进展的分子机制。同时，构建的肿瘤微环境模型揭示了免疫抑制性细胞通讯网络，这为开发新的联合治疗策略提供了理论基础。例如，如果巨噬细胞通过TGF-β促进肿瘤恶性进展，那么将免疫检查点抑制剂与TGF-β通路抑制剂联合使用可能会产生协同效果。这类精细的单细胞水平分析正在改变我们对癌症的认识和治疗方式。

范例三：非模式昆虫的转录组从头组装与环境适应（de novo 组装 + 功能注释）

许多生态学和进化生物学研究的对象是非模式生物，它们没有高质量的参考基因组可用。在这种情况下，转录组从头组装成为了研究这些物种基因表达的主要手段。一个典型的应用场景是研究淡水湖中某种蜻蜓幼虫在重金属污染水域与清洁水域之间的转录组差异。由于这种蜻蜓没有现成的参考基因组，所有分析都必须基于从头组装的转录本序列进行。

实验设计采用标准的RNA-seq方案，每个处理条件设置三个生物学重复以保证统计功效。测序使用Illumina平台的双端125bp模式，这种读长对于从头组装来说是一个较好的平衡点——足够长以提供足够的重叠信息用于组装，又不会因为过长而增加错误率。每个样本的测序量应足以覆盖大多数中等至高表达的转录本，通常建议至少20-30 million read pairs。

转录组从头组装的第一步是使用Trinity程序对各样本分别进行组装。Trinity是目前最广泛使用的转录组组装工具之一，它采用三步策略：首先将reads组装成contig（Inchworm步骤），然后将相关的contig连接成组件（Chrysalis步骤），最后构建转录本并量化其丰度（Butterfly步骤）。Trinity的优点在于它不需要参考基因组，可以直接从短读段数据重建全长转录本序列。然而，当有多个样本时，分别组装后再合并的策略往往能获得更完整的结果，因为它可以捕获样本特异的转录本。

由于分别组装会产生冗余，下一步是使用EvidentialGene工具将各样本的组装结果整合起来。EvidentialGene采用了一套复杂的算法来评估和选择最佳的代表转录本，它会比较来自不同组装的转录本，去除冗余版本并保留最长最完整的序列。这个过程类似于对组装结果进行去重和精炼，最终产生一个非冗余的转录本集合用于后续分析。

进一步的冗余去除可以使用CD-HIT-EST工具完成。CD-HIT-EST基于序列相似度对转录本进行聚类，默认设置95%的相似度阈值意味着序列一致性超过95%的转录本会被聚为一类，只保留最长的一条作为代表。这一步可以有效减少后续功能注释和定量分析的计算量，同时保持转录本集合的信息完整性。需要注意的是，过于激进的聚类（如使用更高的相似度阈值）可能会将真正的旁系同源转录本错误地合并在一起。

功能注释是理解从头组装转录本生物学意义的关键步骤。Blastx搜索将转录本的蛋白质序列与Swiss-Prot数据库进行比对，基于序列相似性推断转录本的可能功能。KEGG注释可以将转录本映射到代谢通路中，帮助理解参与特定生物学过程的基因集合。InterProScan则通过蛋白质结构域和家族的预测来提供互补的功能信息。这三种方法的结合可以获得较为全面的功能注释覆盖率，但即使如此，非模式生物的转录本通常仍有相当一部分无法获得任何功能注释，这是该领域的一个普遍挑战。

在没有参考基因组的情况下，表达定量必须基于组装的转录本本身进行。Salmon是一种快速准确的准映射定量工具，它不需要传统的比对步骤，而是直接将reads映射到转录本序列上并估计其丰度。首先需要用Salmon index命令为组装的转录本建立索引，然后对每个样本运行Salmon quant来获得TPM（Transcripts Per Million）和估计计数值。Salmon的优势在于速度快且内存效率高，适合处理大量样本的定量任务。

差异表达分析同样可以使用DESeq2来完成，但需要先通过tximport包将Salmon在转录本水平上的定量结果汇总到基因水平。这是因为DESeq2的设计是基于基因水平的计数数据，而Salmon默认输出的是转录本水平的定量。tximport可以根据转录本与基因的对应关系将转录本水平的估计值求和，同时正确传播不确定性。之后的标准DESeq2流程与有参考基因组的情况完全相同，包括归一化、离散度估计、假设检验和多重检验校正等功能富集分析的结果直接回答了这个环境适应性研究的核心问题。在重金属污染条件下，蜻蜓幼虫的金属硫蛋白（Metallothionein, MT）基因显著上调，金属硫蛋白是一类富含半胱氨酸的小分子蛋白质，能够结合重金属离子从而降低其毒性。谷胱甘肽-S-转移酶（Glutathione S-transferase, GST）的上调则反映了解毒系统的激活，GST参与谷胱甘肽结合反应，是II相解毒反应的关键酶类。这些分子标志物的发现不仅从机理上解释了蜻蜓对重金属污染的耐受能力，也为环境监测提供了潜在的生物标志物——通过检测野生种群中这些基因的表达水平，可以评估水体重金属污染的程度及其生态影响。

# 非模式生物转录组从头组装流程

# 1. Trinity组装 (对每个样本分别执行)
Trinity --seqType fq --max_memory 50G \
    --left sample_1_R1.fq.gz --right sample_1_R2.fq.gz \
    --CPU 12 --sample_d sample_1 \
    --output trinity_sample_1

# 2. EvidentialGene整合多个组装结果
perl tr2aacds.pl -mrnast sample1.fasta sample2.fasta sample3.fasta -cdhit 0.95

# 3. CD-HIT-EST去冗余 (95%相似度)
cd-hit-est -i combined_transcripts.fasta -o nr_transcripts.fasta \
    -c 0.95 -n 10 -M 0 -T 12

# 4. Salmon索引和定量
salmon index -t nr_transcripts.fasta -i salmon_index --type quasi -k 31
salmon quant -i salmon_index -l A \
    -1 clean_1_R1.fq.gz -2 clean_1_R2.fq.gz \
    -p 8 --validateMappings -o salmon_quant/sample_1

# R语言分析流程

# 5. tximport汇总到基因水平 + DESeq2差异表达
library(tximport)
library(readr)
library(DESeq2)

# 读取Salmon定量结果
files <- list.files("salmon_quant", pattern = "quant.sf$", full.names = TRUE, recursive = TRUE)
names(files) <- sub("/quant.sf", "", basename(dirname(files)))

tx2gene <- read.csv("transcript_to_gene_mapping.csv")  # 转录本-基因对应关系
txi <- tximport(files, type = "salmon", tx2gene = tx2gene)

# DESeq2分析
sample_info <- data.frame(
    condition = factor(c(rep("clean", 3), rep("polluted", 3)))
)
dds <- DESeqDataSetFromTximport(txi, sample_info, design = ~ condition)
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "polluted", "clean"))

# 6. 功能富集分析
sig_up <- subset(res, padj < 0.05 & log2FoldChange > 1)
# 对上调基因进行GO/KEGG富集，重点关注MT、GST等解毒相关基因

这个案例的意义在于它展示了如何在没有参考基因组的情况下开展完整的转录组学研究。从环境科学的角度看，这项研究从分子水平解析了水生生物对重金属污染的耐受机制，发现的MT和GST等分子标志物可以作为环境监测的生物指标。从方法论的角度看，这个工作流程（Trinity组装-EvidentialGene整合-CD-HIT去冗余-Salmon定量-DESeq2分析）已经成为非模式生物转录组学的标准范式，可以被推广应用到其他无参考基因组的物种研究中去。

范例四：lncRNA 在心肌分化中的功能研究（lncRNA 分析 + 共表达网络）

长链非编码RNA（lncRNA）是一类长度超过200个核苷酸但不编码蛋白质的RNA分子，近年来被发现参与了多种重要的生物学过程，包括细胞分化、发育和疾病发生。在人诱导多能干细胞向心肌细胞分化的过程中，动态变化的lncRNA可能扮演着关键调控角色。这个案例的研究目标是系统鉴定分化过程中差异表达的lncRNA，并通过共表达网络分析预测其功能和作用机制。

实验设计采用时间序列采样策略，在iPSC心肌分化过程的四个时间点收集样本进行RNA-seq。这种设计可以捕捉lncRNA表达的时间动态模式，区分早期应答型和晚期效应型的lncRNA。为了确保能够检测到lncRNA信号，测序文库制备采用了总RNA提取加核糖体RNA去除的方案，而不是polyA富集方案，因为许多lncRNA缺少polyA尾，会被polyA富集步骤丢失。

比对和定量的第一步是使用STAR将clean reads比对到人类参考基因组hg38上。与mRNA分析不同的是，lncRNA分析需要使用包含完整lncRNA注释的GTF文件。Ensembl数据库提供的GTF注释包含了biotype字段，可以用来筛选lncRNA类型的转录本，包括lincRNA（基因间区lncRNA）、antisense（反义lncRNA）、processed_transcript等多种亚型。featureCounts或HTSeq-count可以基于这个扩展的GTF文件同时对mRNA和lncRNA进行定量。

lncRNA差异表达分析使用DESeq2方法，与mRNA分析的流程基本相同。在这个心肌分化案例中，鉴定到了300个显著差异表达的lncRNA。对这些lncRNA进行聚类分析可以发现不同的表达模式：有些lncRNA在分化早期达到峰值然后下降，有些随着分化进程持续上升，有些呈现中间态高峰的模式。总共可以归纳为6种主要的表达模式类别，每一类可能对应着分化过程中不同的调控阶段或功能角色。

加权基因共表达网络分析（WGCNA）是挖掘lncRNA功能的强大工具。WGCNA的基本思想是将表达模式相似的基因聚集到同一个模块中，假设同一模块内的基因可能参与相同的生物学过程或受到共同的调控。通过构建lncRNA和mRNA的混合共表达网络，可以识别出与lncRNA共表达的mRNA模块，进而通过模块内mRNA的功能富集来推断关联lncRNA的潜在功能。此外，还可以计算模块特征向量与表型性状（如分化时间点）之间的相关性，找出与特定生物学过程最相关的模块。

模块功能富集分析是解读共表达网络的关键步骤。在这个案例中，某些lncRNA所在的共表达模块显著富集了心脏发育相关的通路，如NOTCH信号通路和WNT信号通路。NOTCH和WNT都是胚胎发育和细胞命运决定中高度保守的信号通路，在心脏发育过程中发挥着重要作用。这些富集结果强烈暗示对应的lncRNA可能通过调控这些通路参与心肌分化的过程。

基于上述分析，研究人员锁定了一个名为LINC00882的心脏特异性lncRNA作为候选功能分子。这个lncRNA的特征包括：在心肌分化过程中显著上调，在成熟心肌细胞中高表达，位于与心脏发育相关的WGCNA模块中，且其表达模式与关键心肌转录因子（如NKX2-5、GATA4）呈正相关。这些证据使其成为一个值得深入研究的功能候选者。

体外功能验证实验是确认lncRNA功能的最终步骤。通过siRNA或shRNA介导的敲降实验降低LINC00882在分化中的iPSC中的表达水平，然后检测心肌标志物（如TNNT2、MYH6、MYH7）的变化。如果LINC00882确实在心肌分化中起正向调控作用，敲降后应该观察到心肌标志物的下调以及分化效率的降低。这种实验验证将计算预测转化为实验事实，确立lncRNA的真实功能角色。

# lncRNA分析与共表达网络分析流程

# 1. lncRNA注释提取 (从Ensembl GTF)
library(rtracklayer)
gtf <- import.gtf("gencode.v44.annotation.gtf")
lncrna_gtf <- gtf[gtf$gene_type %in% c("lncRNA", "3prime_overlapping_ncRNA",
                                        "antisense", "lincRNA",
                                        "processed_transcript")]
write.table(as.data.frame(lncrna_gtf), "lncrna_annotation.gtf",
            sep="\t", quote=FALSE, row.names=FALSE)

# 2. DESeq2 lncRNA差异表达分析
library(DESeq2)
count_mat <- read.csv("lncrna_counts_matrix.csv", row.names = 1)
col_info <- data.frame(
    timepoint = factor(c(rep("D0", 3), rep("D3", 3), rep("D7", 3), rep("D14", 3)))
)
dds <- DESeqDataSetFromMatrix(count_mat, col_info, design = ~ timepoint)
dds <- DESeq(dds)
res_lnc <- results(dds, name = "timepoint_D14_vs_D0")
sig_lncrna <- subset(res_lnc, padj < 0.05)
cat(sprintf("鉴定到 %d 个差异表达lncRNA\n", nrow(sig_lncrna)))

# 3. WGCNA共表达网络分析
library(WGCNA)
expr_data <- vst(assay(dds))  # 变换后的表达矩阵
datExpr <- t(expr_data)

# 选择软阈值
sft <- pickSoftThreshold(datExpr, powerVector = c(1:20))
soft_power <- sft$powerEstimate

# 构建共表达网络
net <- blockwiseModules(datExpr, power = soft_power,
                        TOMType = "unsigned", minModuleSize = 30,
                        reassignThreshold = 0, mergeCutHeight = 0.25,
                        numericLabels = TRUE, pamRespectsDendro = FALSE,
                        verbose = 3)

# 模块与分化时间的关联
module_colors <- labels2colors(net$colors)
MEs <- net$MEs
time_trait <- as.numeric(col_info$timepoint)
module_trait_cor <- cor(MEs, time_trait, use = "p")

# 4. 找到LINC00882所在模块并做功能富集
target_module <- which(net$colors == net$colors["LINC00882"])
module_genes <- names(which(net$colors == target_module))

这个案例的意义在于它提供了一个从计算预测到实验验证的完整lncRNA功能研究范式。通过整合差异表达分析、共表达网络和功能富集等多层次信息，研究者可以从数千个lncRNA中高效筛选出最有可能具有重要功能的候选分子。LINC00882的例子说明，即使在非编码区域也存在对细胞分化至关重要的调控元件。这类研究不仅增进了我们对基因表达调控网络复杂性的认识，还为心血管疾病的机制研究和潜在治疗靶点发现提供了新思路。更重要的是，这套分析框架可以推广应用到其他发育过程或疾病状态的lncRNA功能筛选中去。

范例五：胶质瘤中异常可变剪接与药物敏感性（可变剪接 + 差异剪接 + 功能关联）

可变剪接是真核基因表达调控的重要层面，越来越多的证据表明异常剪接与癌症的发生发展及治疗反应密切相关。胶质母细胞瘤是最具侵袭性的原发性脑肿瘤，替莫唑胺是其标准化疗药物，但耐药性的产生严重限制了治疗效果。这个案例利用TCGA公共数据库中的胶质母细胞瘤RNA-seq数据和临床随访信息，寻找与替莫唑胺耐药相关的可变剪接事件，旨在发现新的预后生物标志物和潜在的治疗靶点。

数据获取是这个分析的第一步。TCGA（The Cancer Genome Atlas）项目提供了大规模的癌症组学数据，包括基因表达、突变、拷贝数变异和临床信息等。对于胶质母细胞瘤，可以从Genomic Data Commons (GDC) 数据门户下载RNA-seq的BAM格式比对文件和对应的临床 metadata。临床数据中包含了患者的治疗信息和生存随访数据，特别是替莫唑胺治疗反应和总生存期等信息，这些是后续关联分析的关键变量。下载时需要注意选择匹配的样本，确保RNA-seq数据和临床数据来自同一批患者。

剪接事件定量是从RNA-seq数据中提取剪接信息的核心步骤。MAJIQ（Modeling Alternative Junction Inclusion Quantification）是一种专门用于检测和定量可变剪接事件的工具。MAJIQ的独特之处在于它不是简单地计算已知的剪接事件，而是通过构建剪接图谱（splicing graphs）来建模所有观测到的剪接连接，然后估算每个样本中各种剪接形式的PSI（Percent Spliced In）值。PSI表示某个剪接变异体在某个基因的所有转录本中所占的比例，范围从0到100%。MAJIQ可以直接从BAM文件中读取比对信息，无需重新比对原始reads，这在处理大量公共数据时非常方便。

差异剪接事件的筛选需要将患者按照治疗反应分组。根据临床数据，将患者分为响应组（对替莫唑胺敏感，生存期较长）和耐药组（对替莫唑胺不敏感，生存期较短）。对每个剪接事件，在两组之间进行统计学检验，常用的是t检验或Wilcoxon秩和检验，然后进行多重检验校正（如Benjamini-Hochberg FDR校正）。在这个案例中，经过筛选获得了150个与替莫唑胺耐药性显著相关的差异剪接事件，这些事件在耐药组和响应组之间存在一致的PSI差异。

事件功能注释是将剪接事件与其宿主基因联系起来的步骤。SUPPA2（Speedy Processing of Annotations for Alternative Splicing Analysis）是一个快速的事件定量和注释工具，它可以将检测到的剪接事件映射到对应的基因上，并提供事件的类型分类（如外显子跳跃、内含子保留、可变供体位点、可变受体位点等）。通过SUPPA2的注释，研究者发现MGMT基因的内含子保留事件在替莫唑胺耐药组中显著富集。MGMT（O6-甲基鸟嘌呤-DNA甲基转移酶）是一个已知的DNA修复酶，它的表达水平与替莫唑胺的疗效密切相关，因为MGMT可以修复替莫唑胺造成的DNA损伤从而使癌细胞存活。

蛋白截断预测是评估剪接事件功能后果的重要环节。内含子保留事件通常会引入过早终止密码子（premature termination codon, PTC），导致产生的mRNA被无义介导的降解（Nonsense-Mediated Decay, NMD）途径识别和降解。这意味着虽然MGMT基因可能在转录水平上看起来正常，但由于异常剪接导致的内含子保留，实际能够翻译成功能性蛋白的mRNA大大减少了。有趣的是，这与MGMT启动子甲基化的效果类似——两者都导致功能性MGMT蛋白的减少，但通过完全不同的分子机制。然而在这个案例中，内含子保留似乎与耐药性相关而非敏感性，这可能提示存在更为复杂的调控关系或其他补偿机制。

生存分析是评估剪接事件临床价值的关键步骤。将患者按照某个剪接事件的PSI值分为高风险组和低风险组（如以中位数作为分界点），然后绘制Kaplan-Meier生存曲线并计算log-rank检验p值。高风险PSI组（即具有更多异常剪接的患者）显示出显著缩短的总生存期，这证明该剪接事件不仅是与耐药性相关的分子标志物，还具有重要的预后价值。Cox比例风险回归模型可以进一步量化风险比（hazard ratio），并在调整其他临床混杂因素后评估剪接事件的独立预后价值。

剪接因子关联分析旨在寻找可能调控这些耐药相关剪接事件的RNA结合蛋白。剪接因子是一类调节可变剪接过程的蛋白质，它们的异常表达或突变可能导致全局性的剪接模式改变。通过相关性分析或机器学习方法（如随机森林），可以鉴定出与耐药相关剪接事件PSI值变化最相关的剪接因子表达量。在这个案例中，RBFOX2（RNA Binding Fox-1 Homolog 2）被鉴定为一个候选的调控因子。RBFOX2是一个组织特异性的剪接因子，在神经组织和肌肉组织中高表达，它已被报道参与多种癌症的剪接重组。如果RBFOX2确实调控了这些耐药相关的剪接事件，那么靶向RBFOX2或其下游效应分子可能成为克服替莫唑胺耐药的新策略。

# TCGA胶质母细胞瘤可变剪接与耐药性分析

# 1. 读取MAJIQ输出的PSI矩阵
psi_data <- read.csv("majiq_psi_matrix.csv", row.names = 1)
clinical <- read.csv("tcga_gbmlgg_clinical.csv")

# 匹配样本并分组
common_samples <- intersect(colnames(psi_data), clinical$sample_id)
psi_subset <- psi_data[, common_samples]
clinical_subset <- clinical[clinical$sample_id %in% common_samples, ]

# 根据治疗反应分组 (响应 vs 耐药)
clinical_subset$response_group <- ifelse(
    clinical_subset$survival_months > median(clinical_subset$survival_months),
    "responder", "non_responder"
)

# 2. 差异剪接事件筛选 (t检验 + FDR校正)
diff_splicing_results <- apply(psi_subset, 1, function(event_psi) {
    responder_psi <- event_psi[clinical_subset$response_group == "responder"]
    non_responder_psi <- event_psi[clinical_subset$response_group == "non_responder"]
    t_test <- t.test(responder_psi, non_responder_psi)
    return(c(mean_diff = mean(responder_psi) - mean(non_responder_psi),
             p_value = t_test$p.value))
})
diff_df <- as.data.frame(t(diff_splicing_results))
diff_df$padj <- p.adjust(diff_df$p_value, method = "BH")
sig_events <- subset(diff_df, padj < 0.05 & abs(mean_diff) > 0.1)
cat(sprintf("鉴定到 %d 个显著差异剪接事件\n", nrow(sig_events)))

# 3. 生存分析 (Kaplan-Meier + Cox回归)
library(survival)
library(survminer)

event_name <- "MGMT_intron_retention"
event_psi <- psi_subset[event_name, ]
risk_group <- ifelse(event_psi > median(event_psi), "high_PSI", "low_PSI")

surv_obj <- Surv(clinical_subset$survival_months, clinical_subset$status) ~ risk_group
fit <- survfit(surv_obj, data = clinical_subset)
ggsurvplot(fit, pval = TRUE, risk.table = TRUE,
           title = paste(event_name, "- 生存分析"))

# Cox回归
cox_model <- coxph(surv_obj, data = clinical_subset)
summary(cox_model)

# MAJIQ可变剪接定量 (从BAM文件)
# 1. 构建剪接图谱
majiq build -gff /path/to/gencode.v44.annotation.gtf \
    -j bam_list.txt -c config.ini -o majiq_output -junc 50 -read 20 -l 200

# 2. PSI定量
majiq quantitate -grp grp1:resp_bam_list.txt -grp grp2:nonresp_bam_list.txt \
    -o majiq_output -config config.ini -n 10000

# 3. 差异剪接检测
majiq deltapsi -grp grp1 -grp grp2 -junc 50 -read 20 -prob 0.2 \
    -o majiq_output/deltapsi -config config.ini

这个案例的意义体现在多个层面。从临床角度看，发现的MGMT内含子保留事件和其他耐药相关剪接事件可以作为预测替莫唑胺疗效的生物标志物，帮助医生在治疗前就判断患者可能的反应，从而制定个体化的治疗方案。从机制角度看，剪接因子RBFOX2的鉴定为理解耐药性的分子基础提供了新的切入点，靶向RBFOX2或其调控的剪接事件可能成为克服耐药的联合治疗策略。从方法论角度看，这个工作流程展示了如何充分利用公共组学数据资源进行转化医学研究，将计算生物学发现与临床问题紧密联系起来。可变剪接作为一个新兴的研究维度，正在为癌症生物学和精准医学带来全新的视角和机遇。