转录组测序来源的表达量矩阵数据分析大家都耳熟能详了，从基因表达数据的差异分析到功能注释和通路分析的整个分析流程是理解基因表达变化及其生物学意义的重要步骤。以下是详细步骤的介绍：
基因表达量矩阵的准备：
首先，需要获得全部基因的mRNA表达量矩阵，这通常来自于RNA测序（RNA-Seq）或蛋白质组学数据。
对原始数据进行质量控制、标准化和过滤，以确保数据的可靠性和可比性。
使用统计方法（如t检验、ANOVA、Limma、DESeq2等）比较两组样本（例如，治疗组与对照组）之间的基因表达差异，以识别差异表达基因（DEGs）。
根据差异表达分析的结果，使用合适的统计学阈值去确定上调（表达量增加）和下调（表达量减少）的基因。
使用富集分析工具（如DAVID、clusterProfiler等）来确定差异表达基因是否在特定的GO术语或KEGG通路中富集，这有助于揭示基因表达变化背后的生物学过程。
基因集富集分析（GSEA）是一种方法，用于评估预定义的基因集（如KEGG通路或GO术语）在排序列表中的富集情况，这有助于识别与特定生物学状态或条件相关的基因集。在GSEA中，基因的重要性排序通常基于它们的差异表达程度和统计显著性。
使用实验方法（如qPCR、Western blot、细胞实验等）来验证差异表达分析和GSEA的结果。
利用图表（如火山图、热图、富集图等）来直观展示差异表达基因和富集分析的结果。
将基因表达数据与其他类型的数据（如表观遗传学、蛋白质互作网络等）整合，以获得更全面的生物学见解。
整个流程涉及从原始数据到生物学解释的多个步骤，每一步都对最终的生物学结论至关重要。通过这些分析，研究人员可以更好地理解基因表达变化背后的分子机制，并为进一步的实验研究提供方向。
一个简单的转录组测序来源的表达量矩阵大概如下所示：
Sample1 Sample2 Sample3
Gene1 10 15 20
Gene2 5 8 12
Gene3 18 22 25
Gene4 7 9 11
Gene5 3 6 9
这个表达量矩阵可以用作演示如何进行差异表达分析、数据可视化或其他生物信息学分析的起点。在实际应用中，表达量矩阵可能会包含成千上万的基因和样本，并且数据通常存储在文件中，如CSV、TXT或专用的生物信息学文件格式。
可以使用一个简单的r代码，作为非常简单的表达量矩阵示例，它包含了5个基因在3个样本中的表达量数据。这个矩阵可以用作差异表达分析或其他类型的基因表达分析的基础。
# 创建一个简单的表达量矩阵示例
expression_matrix <- matrix(c(
10, 15, 20, # 基因1在样本1、2、3中的表达量
5, 8, 12, # 基因2在样本1、2、3中的表达量
18, 22, 25, # 基因3在样本1、2、3中的表达量
7, 9, 11, # 基因4在样本1、2、3中的表达量
3, 6, 9 # 基因5在样本1、2、3中的表达量
), nrow = 5, byrow = TRUE, dimnames = list(c("Gene1", "Gene2", "Gene3", "Gene4", "Gene5"), c("Sample1", "Sample2", "Sample3")))
# 查看表达量矩阵
print(expression_matrix)
这段代码创建了一个5行3列的矩阵，其中行代表基因，列代表样本。nrow = 5 指定了矩阵的行数，byrow = TRUE 表示矩阵是按行填充的，dimnames 为矩阵的行和列提供了名称。不过，绝大部分情况下，我们是直接下载网络的表达量矩阵，或者读取自己的txt或者csv文件啦，不需要上面的每个基因在每个样品的手动创建！比如这个学徒作业：不要简单的相信作者提供的表达量矩阵，针对这个GSE13904数据集里面的找到脓毒症和正常对照，这两个分组的样品，就是简单的从网络下载。
肿瘤外显子测序后的突变与否矩阵
上面的转录组测序表达量矩阵，如果是最原始的count值，每个基因在每个样品的值动态范围可以很大很大。但是在肿瘤外显子测序（Exome Sequencing）中，突变与否矩阵通常记录了特定基因位点在样本中的突变状态。这种矩阵可以用二进制值（0和1）来表示，其中0可能表示野生型（未突变），1表示突变型。
以下是一个简单的R代码示例，用于创建一个肿瘤外显子测序后的突变矩阵：
# 创建一个简单的突变矩阵示例
# 假设我们有5个样本和3个基因位点
mutation_matrix <- matrix(c(
0, 1, 0, # 样本1的基因位点1、2、3的突变状态
1, 0, 1, # 样本2的基因位点1、2、3的突变状态
0, 1, 0, # 样本3的基因位点1、2、3的突变状态
1, 1, 0, # 样本4的基因位点1、2、3的突变状态
0, 0, 1 # 样本5的基因位点1、2、3的突变状态
), nrow = 3, byrow = TRUE, dimnames = list(c("Gene1", "Gene2", "Gene3"), c("Sample1", "Sample2", "Sample3", "Sample4", "Sample5")))
# 查看突变矩阵
print(mutation_matrix)
这段代码创建了一个3行5列的矩阵，其中行代表基因位点，列代表样本。nrow = 3 指定了矩阵的行数，byrow = TRUE 表示矩阵是按行填充的，dimnames 为矩阵的行和列提供了名称。
输出结果将如下所示：
Sample1 Sample2 Sample3 Sample4 Sample5
Gene1 0 1 0 1 0
Gene2 1 0 1 1 0
Gene3 0 1 0 0 1
在这个矩阵中，0 表示该基因位点在相应样本中未检测到突变，1 表示检测到突变。这种类型的矩阵可以用于后续的突变分析，如突变负荷分析、突变谱分析或与临床数据的相关性分析。在实际应用中，突变矩阵可能包含更多的样本和基因位点，并且数据通常来源于生物信息学分析流程，如变异检测工具（如GATK、MuTect等）的输出。
而且，因为绝大部分基因在绝大部分样品其实都是野生型，根据TMB概念，每个肿瘤病人的2万多个编码基因发生somatic突变的基因大多数在100到1000之间，也就是说其它的19000个基因都是野生型的，这样的话上面的肿瘤外显子测序后的突变与否矩阵就太稀疏了。所以肿瘤外显子测序（Exome Sequencing）后得到的突变数据通常保存在特定的文件格式中，以便于存储和分析。MAF（Mutation Annotation Format）和VCF（Variant Call Format）是两种常用的突变数据文件格式：
MAF（Mutation Annotation Format）：
MAF是一种用于存储突变信息的文本文件格式，它由TCGA（The Cancer Genome Atlas）项目推广使用。 MAF文件包含了突变的详细注释信息，包括染色体位置、突变类型、氨基酸变化、基因名称、样本信息等。 这种格式便于在不同研究之间共享和比较突变数据，因为它提供了一种标准化的方式来记录突变信息。
VCF是一种用于存储基因组变异信息的文本文件格式，被广泛用于高通量测序数据的变异检测结果。 VCF文件包含了关于变异（如SNPs、插入、缺失等）的详细信息，包括变异的位置、参考和替代等位基因、质量分数、过滤状态等。 VCF格式支持多种类型的变异注释，并且可以包含多个样本的变异数据。
MAF和VCF通常是记录了哪些基因在哪个样品有突变即可，没有突变的基因就压根不会记录。
肿瘤外显子测序后的突变与否矩阵也可以差异分析吗
大部分情况下， 大家都会从公司那边拿到了，肿瘤外显子测序后的突变与否信息，可能是MAF和VCF格式。如果万一大家拿到了的是fastq格式的测序数据，就需要服务器走下面的流程了：
肿瘤外显子数据处理系列教程（九）找 Somatic Mutations的 Signature 图谱 肿瘤外显子数据处理系列教程（九）拷贝数变异分析（不同软件的比较） 肿瘤外显子数据处理系列教程（九）拷贝数变异分析（主要是GATK） 肿瘤外显子数据处理系列教程（八）不同注释软件的比较(下)：可视化比较maf文件 肿瘤外显子数据处理系列教程（八）不同注释软件的比较（中）：注释后转成maf文件 肿瘤外显子数据处理系列教程（八）不同注释软件的比较（上）：安装及使用 肿瘤外显子数据处理系列教程（七）maftools可视化 肿瘤外显子数据处理系列教程（六）vcf文件的注释及ANNOVAR的使用 肿瘤外显子数据处理系列教程（五）GATK的最佳实践 肿瘤外显子数据处理系列教程（番外篇）bam文件载入igv可视化 肿瘤外显子数据处理系列教程（四）比对结果的质控 肿瘤外显子数据处理系列教程（三）比对 肿瘤外显子数据处理系列教程（二）质控与去接头 肿瘤外显子数据处理系列教程（一）读文献并且下载测序数据
如果是，比如这个2024的文章：《PICK1 inhibits the malignancy of nasopharyngeal carcinoma and serves as a novel prognostic marker》，就 Landscape of somatic mutation profiles in 13 NPC samples. 分成了 NDM (no distant metastasis) and DM (distant metastasis) groups. 如下所示：
就是一个很简单的差异分析，然后略微有点莫名其妙的就定位到了PICK1基因。而如果是要实现上面的图超级简单，把上面的两分组各自独立做成maf格式文件然后分开读取，使用这个函数即可：https://rdrr.io/bioc/maftools/man/mafCompare.html
primary.apl <- system.file("extdata", "APL_primary.maf.gz", package = "maftools")
relapse.apl <- system.file("extdata", "APL_relapse.maf.gz", package = "maftools")
primary.apl <- read.maf(maf = primary.apl)
relapse.apl <- read.maf(maf = relapse.apl)
pt.vs.rt <- mafCompare(m1 = primary.apl, m2 = relapse.apl, m1Name = 'Primary',
m2Name = 'Relapse', minMut = 5)
然后就可以看出来了：
## $results
## Hugo_Symbol Primary Relapse pval or ci.up ci.low
## <char> <num> <int> <num> <num> <num> <num>
## 1: PML 1 11 1.529935e-05 0.03537381 0.2552937 0.000806034
## 2: RARA 0 7 2.574810e-04 0.00000000 0.3006159 0.000000000
## 3: RUNX1 1 5 1.310500e-02 0.08740567 0.8076265 0.001813280
## 4: FLT3 26 4 1.812779e-02 3.56086275 14.7701728 1.149009169
## 5: ARID1B 5 8 2.758396e-02 0.26480490 0.9698686 0.064804160
## 6: WT1 20 14 2.229087e-01 0.60619329 1.4223101 0.263440988
## 7: KRAS 6 1 4.334067e-01 2.88486293 135.5393108 0.337679367
## 8: NRAS 15 4 4.353567e-01 1.85209500 8.0373994 0.553883512
## 9: ARID1A 7 4 7.457274e-01 0.80869223 3.9297309 0.195710173
## adjPval
## <num>
## 1: 0.0001376942
## 2: 0.0011586643
## 3: 0.0393149868
## 4: 0.0407875250
## 5: 0.0496511201
## 6: 0.3343630535
## 7: 0.4897762916
## 8: 0.4897762916
## 9: 0.7457273717
##
## $SampleSummary
## Cohort SampleSize
## <char> <num>
## 1: Primary 124
## 2: Relapse 58
但是，这一切的结果的合理性的前提都是样品数量。关于上面的对比结果，其实文档里面 有很多 对比的可视化函数：
9.5 Comparing two cohorts (MAFs) 9.5.1 Forest plots 9.5.2 Co-onco plots 9.5.3 Co-bar plots
