现在的单细胞数据确实很全面，但是几乎所有的单细胞数据都是不带随访的信息的。这样的话，要在单细胞层面实现基因的预后分析不太现实。因此，我们还是需要“温故”，TCGA是不二选择啦~
假设现在我已经根据单细胞分析获得了一批基因，这些基因可能是与巨噬细胞或者中性粒细胞有关的特异性基因。接下来我想利用这些基因看看它们在TNBC中的预后意义。
于是：
获取TNBC的临床数据
library(cBioPortalData)
library(survival)
library(survminer)
library(TCGAbiolinks)
library(dplyr)
library(ggplot2)
cbio
hostname = "www.cbioportal.org",
protocol = "https",
api. = "/api" # 更新 API 路径
)
clinical "brca_tcga")
# 筛选 TNBC 样本
tnbc_samples
clinical$ER_STATUS_BY_IHC == "Negative" &
clinical$PR_STATUS_BY_IHC == "Negative" &
clinical$IHC_HER2 == "Negative",
]
# 查看 TNBC 样本数量
# 去除 patientId 为 NA 的行
clin_data
save(clin_data,file = "TNBC_clin.Rdata")
获取TNBC的表达矩阵
# 创建查询对象
query
project = "TCGA-BRCA",
data.category = "Transcriptome Profiling",
data.type = "Gene Expression Quantification",
workflow.type = "STAR - Counts")
# 启用多线程下载
GDCdownload(query,
files.per.chunk = 10,
directory = "GDCdata")
data
expre_matrix
save(expre_matrix,file = "BRCA_expre_matrix.Rdata")
筛选基因表达矩阵
# 过滤出 TNBC 样本的表达数据
# 提取简化的 patientId
colnames(expre_matrix) 1, 12)
# 查看转换后的列名
colnames(expre_matrix)[1:4]
# 筛选出在 clin_data 中存在的样本
matched_samples in% clin_data$patientId
exp_matrix
# 查看筛选后的样本数量
print(paste("匹配的样本数量:", ncol(exp_matrix)))
exp_mat ##去除重复样本
这里有一步额外的基因转换：
# 加载 biomaRt 包
library(biomaRt)
# 连接到 Ensembl 数据库
ensembl "ensembl", dataset = "hsapiens_gene_ensembl")
# 获取基因 ID 列表
rownames(exp_mat) "\\.[0-9]+$", "", rownames(exp_mat))
gene_ids
# 定义需要获取的属性
attributes "ensembl_gene_id", "hgnc_symbol")
# 使用 getBM() 进行转换
gene_info
attributes = attributes,
filters = "ensembl_gene_id",
values = gene_ids,
mart = ensembl
)
# 将 gene_info 与 exp_mat 合并
exp_mat
exp_mat$ensembl_gene_id
# 使用 left_join() 合并数据
exp_mat "ensembl_gene_id")
exp_mat=exp_mat[!duplicated(exp_mat$hgnc_symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果s
# 去除 hgnc_symbol 列中的 NA 值
exp_mat
rownames(exp_mat)
# 去除多余的列
exp_mat in% c("ensembl_gene_id", "hgnc_symbol")]
# 查看转换后的表达矩阵
exp_mat[1:4, 1:4]
# 对表达量矩阵取 log2
log_exp_mat 1) # 加 1 是为了避免对 0 取对数
# 查看转换后的数据
range(log_exp_mat)
save(log_exp_mat,file = "gene_exp_mat.Rdata")
# 假设你的基因列表是一个字符向量
