我们生信技能树收到一个粉丝的提问:为什么直接做 HNSC 的 TCGA 差异分析类似 GAPDH 这样的 housekeeper基因 差异会这么明显呢?我们先去看看 管家基因到底会不会差异表达吧?
人类 housekeeper基因有哪些
首先,我们来问一下人工智能大模型kimi:人类 housekeeper基因有哪些,介绍并生成一个R向量
"Housekeeper genes"(管家基因)是指在细胞中持续表达,对维持基本细胞功能至关重要的基因。这些基因的产物通常参与基本的细胞过程,如DNA复制、转录、翻译和细胞代谢等,它们的表达通常不受外部条件或环境因素的显著影响。
在人类基因组中,有许多基因可以被认为是管家基因。以下是一些常见的人类管家基因的例子:
ACTB:编码β-肌动蛋白(Beta-actin),是细胞骨架的重要组成部分。 GAPDH:编码甘油醛-3-磷酸脱氢酶(Glyceraldehyde-3-phosphate dehydrogenase),在糖酵解过程中起关键作用。 HSP90:编码热休克蛋白90(Heat shock protein 90),参与蛋白质折叠和稳定。 RPL13A:编码核糖体蛋白L13a,是蛋白质合成的必需成分。 TBP:编码TATA盒结合蛋白(TATA box-binding protein),在转录起始中起关键作用。 UBC:编码泛素C(Ubiquitin C),参与蛋白质降解和调控。 YWHAG:编码14-3-3蛋白gamma(14-3-3 protein gamma),参与信号转导和细胞周期调控。 以下是这些基因名称的R向量:
housekeeper_genes <- c("ACTB", "GAPDH", "HSP90", "RPL13A", "TBP", "UBC", "YWHAG")
请注意,管家基因的列表并不是固定的,不同的研究可能会根据实验条件和细胞类型识别出不同的管家基因。此外,管家基因的表达水平通常用于标准化基因表达数据,以消除实验中的技术变异。
其中 HSP90 是一个蛋白的名字,它对应的基因其实 有下面 三个:
这三个基因都是编码热休克蛋白90(Heat Shock Protein 90,简称HSP90)家族成员的基因。HSP90是一种高度保守的分子伴侣蛋白,对于维持蛋白质稳态、帮助蛋白质折叠、防止蛋白质聚集以及调控蛋白质降解等过程至关重要。下面是对这三个基因的具体解读:
HSP90AA1:
也被称为HSP90α(alpha),是HSP90家族中的一种主要亚型。 在细胞质中表达,参与广泛的细胞内蛋白质折叠和稳定。 在细胞应激反应中发挥重要作用,如热休克、氧化应激等。 与多种疾病相关,包括癌症、神经退行性疾病等。 HSP90AB1:
也被称为HSP90β(beta),与HSP90α共同存在于细胞质中。 在结构和功能上与HSP90α相似,但表达水平通常较低。 与HSP90α形成异源二聚体,共同参与蛋白质折叠和稳定。 在某些细胞类型和组织中可能具有特定的功能。 HSP90B1:
也被称为HSP90β1或grp94,主要存在于内质网中,而不是细胞质。 参与新合成蛋白质的折叠和质量控制,特别是在内质网中。 在抗原呈递、免疫反应以及某些癌症中发挥重要作用。 与HSP90AA1和HSP90AB1相比,HSP90B1在进化上更为古老,具有不同的调控机制。 总的来说,HSP90家族蛋白在细胞内发挥着关键的分子伴侣功能,对于维持蛋白质稳态和细胞功能至关重要。它们在多种生物学过程和疾病中发挥作用,是重要的研究和治疗靶点。
因此,得到 9个管家基因: c("ACTB", "GAPDH", "HSP90AA1", "HSP90AB1", "HSP90B1", "RPL13A", "TBP", "UBC", "YWHAG")
看 其 TCGA 数据库中HNSC 的差异情况
HNSC 是 "Head and Neck Squamous Cell Carcinoma" 的缩写,中文全称为头颈鳞状细胞癌。关于TCGA数据的下载以及整理,我们前面已经给大家介绍过两个帖子了:
1、HNSC 数据下载
根据上面的帖子,很容易就可以将 HNSC 这个癌症的基因表达count,fpkm等基因表达矩阵下载下来:
## download tcga data
## update: 2024-02-22
## Author: zhang juan
rm(list=ls())
## load packages
library(TCGAbiolinks)
library(SummarizedExperiment)
suppressPackageStartupMessages(library(tidyverse))
# 癌症类型,用 getGDCprojects()$project_id 查看所有
getGDCprojects()$project_id
pro <- "TCGA-HNSC"
# 设置query
query <- GDCquery(
project = "TCGA-HNSC", # 癌症类型,选择 TCGA-HNSC
data.category="Transcriptome Profiling", # 数据类别, 用getProjectSummary(project)查看所有类别
data.type ="Gene Expression Quantification", # 数据类型
workflow.type="STAR - Counts" # 工作流类型
)
## 下载数据
GDCdownload(query=query, files.per.chunk= 50, directory = "./")
## 整理数据并存储为 R对象
GDCprepare(query,save=T,save.filename=paste0(pro,".transcriptome.Rdata"), directory = "./")
## 再次加载
load("TCGA-HNSC/TCGA-HNSC.transcriptome.Rdata")
ls()
names(assays(data))
rowdata <- rowData(data)
table(rowdata$gene_type)
mrna <- data[rowdata$gene_type == "protein_coding",]
mrna_count <- assay(mrna,"unstranded") # mRNA的counts矩阵
mrna_fpkm <- assay(mrna,"fpkm_unstrand") # mRNA的fpkm矩阵
# 添加gene_symbol, 先提取gene_name
symbol_mrna <- rowData(mrna)$gene_name
head(symbol_mrna)
###### count值
mrna_count_symbol <- cbind(data.frame(symbol_mrna), as.data.frame(mrna_count))
# 去重复保留最大的那个
mrna_count_symbol1 <- mrna_count_symbol %>%
as_tibble() %>% # tibble不支持row name,我竟然才发现!
mutate(meanrow = rowMeans(.[,-1]), .before=2) %>%
arrange(desc(meanrow)) %>%
distinct(symbol_mrna,.keep_all=T) %>%
select(-meanrow) %>%
as.data.frame()
mrna_count_symbol1[1:4, 1:5]
rownames(mrna_count_symbol1) <- mrna_count_symbol1$symbol_mrna
mrna_count_symbol1 <- mrna_count_symbol1[,-1]
saveRDS(mrna_count_symbol1, file = "TCGA-HNSC/mrna_count_symbol.rds")
###### fpkm值
mrna_fpkm_symbol <- cbind(data.frame(symbol_mrna), as.data.frame(mrna_fpkm))
# 去重复保留最大的那个
mrna_fpkm_symbol1 <- mrna_fpkm_symbol %>%
as_tibble() %>% # tibble不支持row name,我竟然才发现!
mutate(meanrow = rowMeans(.[,-1]), .before=2) %>%
arrange(desc(meanrow)) %>%
distinct(symbol_mrna,.keep_all=T) %>%
select(-meanrow) %>%
as.data.frame()
mrna_fpkm_symbol1[1:4, 1:5]
rownames(mrna_fpkm_symbol1) <- mrna_fpkm_symbol1$symbol_mrna
mrna_fpkm_symbol1 <- mrna_fpkm_symbol1[,-1]
saveRDS(mrna_fpkm_symbol1, file = "TCGA-HNSC/mrna_fpkm_symbol.rds")
# "HSP90": "HSP90AA1" "HSP90AB1" "HSP90B1"
housekeeper_genes <- c("ACTB", "GAPDH", "HSP90AA1", "HSP90AB1", "HSP90B1", "RPL13A", "TBP", "UBC", "YWHAG")
mrna_count_symbol1[housekeeper_genes, 1:6]
mrna_fpkm_symbol1[housekeeper_genes, 1:6]
2、差异表达分析
随便选一个 差异算法 DESeq2吧:
rm(list = ls())
options(stringsAsFactors = F)
library(ggplot2)
library(ggstatsplot)
library(ggsci)
library(RColorBrewer)
library(patchwork)
library(ggplotify)
library(limma)
library(edgeR)
library(DESeq2)
options(scipen = 100)
symbol_matrix <- readRDS(file = 'TCGA-HNSC/mrna_count_symbol.rds')
symbol_matrix[1:4,1:4]
table(str_sub(colnames(symbol_matrix), 14,15))
table(str_sub(colnames(symbol_matrix), 14,16))
group_list <- str_sub(colnames(symbol_matrix), 14,15)
group_list <- ifelse(group_list=="11","normal","tumor")
group_list <- factor(group_list, levels = c("normal","tumor")) # 对照组在前
table(group_list)
levels(group_list)[2]
levels(group_list)[1]
exprSet <- symbol_matrix
(colData <- data.frame(row.names=colnames(exprSet), group_list=group_list) )
exprSet[, 1:ncol(exprSet)] <- lapply(exprSet[, 1:ncol(exprSet)], as.integer)
class(exprSet)
# 差异分析
dds <- DESeqDataSetFromMatrix(countData = exprSet, colData = colData, design = ~ group_list)
dds <- DESeq(dds)
res <- results(dds, contrast=c("group_list", levels(group_list)[2], levels(group_list)[1]))
res <- res[order(res$padj),]
res <- as.data.frame(res)
head(res)
DEG_deseq2 <- na.omit(res)
save(DEG_deseq2, file= 'TCGA-HNSC/DEG_deseq2.Rdata')
# "HSP90": "HSP90AA1" "HSP90AB1" "HSP90B1"
housekeeper_genes <- c("ACTB", "GAPDH", "HSP90AA1", "HSP90AB1", "HSP90B1", "RPL13A", "TBP", "UBC", "YWHAG")
temp <- DEG_deseq2[housekeeper_genes ,]
temp
火山图如下:
你可以看到,即使使用比较低的阈值进行筛选:FC>1.5, Pvalue<0.05,这些基因基本上都没有差异,除了有一个因为FC在临界点之外。
问题出现在哪里呢?你的分析中 存在 housekeeper基因 发生了显著差异表达吗?
如果你对TCGA等数据库挖掘感兴趣,但又不知道如何入门,也许你可以关注一下下面的课程: