从TCGA数据下载到基因表达相关性分析的完整指南 To "square tiger"

在本文中，我将详细向方虎介绍如何从 TCGA（The Cancer Genome Atlas） 中下载数据，整理和处理这些数据，然后使用 R 语言进行基因表达相关性分析，并绘制散点图以展示结果。并希望他能成为如何在R中进行从数据下载到可视化的完整流程的人。

1. 下载TCGA转录组数据

我们首先需要使用 TCGAbiolinks 包从 TCGA 数据库下载所需的转录组数据。这个包可以直接从 Bioconductor 中获取。以下是下载数据的具体步骤。虽然你喜欢在网页下载，但我还是告诉你在命令行操作显得比较专业，想想别人都在gui页面点点点，你打开r就是一段噼里啪啦，然后进行参数调整，这多装逼

1.1 安装和加载必要的包

在开始之前，请确保安装并加载了以下 R 包：

if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}
BiocManager::install("TCGAbiolinks")
BiocManager::install("SummarizedExperiment")

library(TCGAbiolinks)
library(SummarizedExperiment)

1.2 下载TCGA-BRCA转录组数据

下面的代码示例展示了如何下载乳腺癌（BRCA）数据：如果是lich那你应该也知道替换哪些变量，从project,category,data.type,workflow.type进行修改

查询并下载数据

query <- GDCquery(project = "TCGA-BRCA",
                  data.category = "Transcriptome Profiling",
                  data.type = "Gene Expression Quantification",
                  workflow.type = "STAR - Counts")
GDCdownload(query)

数据准备

data <- GDCprepare(query)

在这一步，我们获取了 TCGA-BRCA 项目的 RNA-seq 数据，并使用 STAR - Counts 作为工作流类型,要注意支持哪些工作流进行选择，一开始不知道选错了看了半天无语😓。

2. 整理和处理数据

下载后的数据需要进一步整理，以便进行相关性分析。在 TCGA 数据中，基因名常用 Ensembl ID，而且通常带有版本号。我们需要去掉这些版本号，并将 Ensembl ID 转换为 基因符号。非常之重要，我一开始直接拿你给我的基因名称去找，发现怎么也找不到，我想这方虎也太不靠谱了，这病里也没这些基因怎么分析啊，真水

2.1 处理 Ensembl ID 和 基因名 转换

去掉 Ensembl ID 的版本号

rownames(data) <- stringr::str_split(rownames(data), "\\.", simplify = TRUE)[, 1]

确保基因名唯一

data <- data[!duplicated(rownames(data)), ]

2.2 提取感兴趣的基因数据

假设我们感兴趣的是 NR1D1 和自噬通路相关的基因，我们需要将这些基因的数据提取出来。将其存储到你定义的字符串中，便于后续处理，怕你不懂告诉你一下

使用 biomaRt 包将基因符号转换为 Ensembl ID

if (!requireNamespace("biomaRt", quietly = TRUE)) {
    BiocManager::install("biomaRt")
}
library(biomaRt)

mart <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")

genes_of_interest <- c("NR1D1", "ATG5", "ATG7", "ATG12", "ATG16", "ULK1", "ULK2", "ULK3",
                       "SQSTM1", "MAP1LC3A", "MAP1LC3B", "BECN1", "XBP1")

gene_mapping <- getBM(attributes = c("hgnc_symbol", "ensembl_gene_id"),
                      filters = "hgnc_symbol",
                      values = genes_of_interest,
                      mart = mart)

提取感兴趣基因的表达数据

genes_filtered <- gene_mapping$ensembl_gene_id[gene_mapping$ensembl_gene_id %in% rownames(data)]
expression_subset <- data[genes_filtered, ]

2.3 数据的 log2 转换

为了处理 RNA-seq 数据，我们通常对表达量进行 log2(x + 1) 转换，这样可以减小极端值的影响，并使数据更符合正态分布。这个应该知道哎

对表达数据进行 log2 转换

expression_subset_log2 <- log2(expression_subset + 1)

转置数据框，使基因为列，样本为行 你用wps做这一步了吗，前几行全是sample，后面又全是基因，人家是矩阵分析，你是人肉文档搜索啊

expression_df <- as.data.frame(t(expression_subset_log2)) `非常重要之一步，到这里前戏就做好了，后面就是用此矩阵开始分析啦`

3. 相关性分析

现在我们有了整理后的数据，我们可以开始进行相关性分析。我们将 NR1D1 的表达与所有自噬通路相关基因的平均表达量进行比较。算通路就算平均量吧，感觉够用了

3.1 计算自噬通路基因的平均表达 三步走，提取，过滤，计算

提取自噬通路基因的表达数据

autophagy_genes <- c("ATG5", "ATG7", "ATG12", "ATG16", "ULK1", "ULK2", "ULK3",
                     "SQSTM1", "MAP1LC3A", "MAP1LC3B", "BECN1", "XBP1")

过滤存在于数据中的基因

autophagy_genes_filtered <- autophagy_genes[autophagy_genes %in% colnames(expression_df)]
autophagy_expression <- expression_df[, autophagy_genes_filtered]

计算自噬通路基因的平均表达值

autophagy_avg_expression <- rowMeans(autophagy_expression, na.rm = TRUE)

4. 绘制相关性散点图

我们将 NR1D1 的表达值与自噬通路的平均表达值进行比较，绘制散点图并计算相关系数。

4.1 绘制散点图

我们使用 ggplot2 和 ggpubr 包来绘制散点图，并添加线性拟合线和相关性统计信息。

安装并加载 ggplot2 和 ggpubr 包

if (!requireNamespace("ggplot2", quietly = TRUE)) {
    install.packages("ggplot2")
}
if (!requireNamespace("ggpubr", quietly = TRUE)) {
    install.packages("ggpubr")
}

library(ggplot2)
library(ggpubr)

创建一个新的数据框

data_for_plot <- data.frame(
  NR1D1 = expression_df$NR1D1,
  Autophagy_Avg = autophagy_avg_expression
)

绘制散点图 这里你可以自我调整绘图，找一个你们常用的格式，我这绘图不一定符合你们文章绘图格式要求

ggplot(data_for_plot, aes(x = Autophagy_Avg, y = NR1D1)) +
  geom_point(color = "blue", alpha = 0.4, size = 2) +
  geom_smooth(method = "lm", se = FALSE, color = "red", size = 1) +
  stat_cor(method = "pearson", label.x.npc = "left", label.y.npc = "top", size = 5) +
  theme_bw() +
  labs(title = "Correlation between Autophagy Pathway Average Expression and NR1D1 in BRCA Samples",
       x = "Autophagy Pathway Average Expression (log2 transformed)",
       y = "NR1D1 Expression (log2 transformed)") +
  theme(
    plot.title = element_text(hjust = 0.5, size = 14, face = "bold"),
    axis.title = element_text(size = 12),
    axis.text = element_text(size = 10)
)

5. 解读结果

在绘制的散点图中，红色的回归线表示了 NR1D1 与自噬通路基因平均表达之间的线性关系。图中还显示了相关系数 r 和 p 值，用于评估两者之间的相关性强度及其统计显著性。

• r 值 = 0.4：表示中等强度的正相关。
• p 值 < 2.2e-16：表示在统计学意义上具有高度显著的相关性，但在解读时要谨慎，可能受样本量影响较大。

总结

在这篇指南中，我们从 TCGA 下载了 BRCA 数据，进行了数据的预处理，包括 基因名转换 和 log2 转换，然后计算了 NR1D1 与自噬通路相关基因的平均表达的相关性，最后通过散点图展示了结果。通过这种分析流程，你可以了解特定基因与生物学通路之间的潜在关系。学会这五步，走上无敌路