housekeeper基因在肿瘤与正常样本中会发生显著差异表达吗？

housekeeper基因在肿瘤与正常样本中会发生显著差异表达吗？

我们生信技能树收到一个粉丝的提问：为什么直接做 HSC 的 TCGA 差异分析类似 GAPDH 这样的 housekeeper基因 差异会这么明显呢？我们先去看看 管家基因到底会不会差异表达吧？

首先，我们来问一下人工智能大模型kimi：人类 housekeeper基因有哪些，介绍并生成一个R向量

"Housekeeper genes"（管家基因）是指在细胞中持续表达，对维持基本细胞功能至关重要的基因。这些基因的产物通常参与基本的细胞过程，如DA复制、转录、翻译和细胞代谢等，它们的表达通常不受外部条件或环境因素的显著影响。

在人类基因组中，有许多基因可以被认为是管家基因。以下是一些常见的人类管家基因的例子：

1. ACTB：编码β-肌动蛋白（Beta-actin），是细胞骨架的重要组成部分。
2. GAPDH：编码甘油醛--磷酸脱氢酶（Glyceraldehyde--phosphate dehydrogenase），在糖酵解过程中起关键作用。
3. HSP90：编码热休克蛋白90（Heat shock protein 90），参与蛋白质折叠和稳定。
4. RPL1A：编码核糖体蛋白L1a，是蛋白质合成的必需成分。
5. TBP：编码TATA盒结合蛋白（TATA box-binding protein），在转录起始中起关键作用。
6. UBC：编码泛素C（Ubiquitin C），参与蛋白质降解和调控。
7. YWHAG：编码14--蛋白gamma（14-- protein gamma），参与信号转导和细胞周期调控。

以下是这些基因名称的R向量：

housekeeper_genes <- c("ACTB", "GAPDH", "HSP90", "RPL1A", "TBP", "UBC", "YWHAG")

请注意，管家基因的列表并不是固定的，不同的研究可能会根据实验条件和细胞类型识别出不同的管家基因。此外，管家基因的表达水平通常用于标准化基因表达数据，以消除实验中的技术变异。

其中 HSP90 是一个蛋白的名字，它对应的基因其实 有下面 三个：

这三个基因都是编码热休克蛋白90（Heat Shock Protein 90，简称HSP90）家族成员的基因。HSP90是一种高度保守的分子伴侣蛋白，对于维持蛋白质稳态、帮助蛋白质折叠、防止蛋白质聚集以及调控蛋白质降解等过程至关重要。下面是对这三个基因的具体解读：

1. HSP90AA1：
   • 也被称为HSP90α（alpha），是HSP90家族中的一种主要亚型。
   • 在细胞质中表达，参与广泛的细胞内蛋白质折叠和稳定。
   • 在细胞应激反应中发挥重要作用，如热休克、氧化应激等。
   • 与多种疾病相关，包括癌症、神经退行性疾病等。
2. HSP90AB1：
   • 也被称为HSP90β（beta），与HSP90α共同存在于细胞质中。
   • 在结构和功能上与HSP90α相似，但表达水平通常较低。
   • 与HSP90α形成异源二聚体，共同参与蛋白质折叠和稳定。
   • 在某些细胞类型和组织中可能具有特定的功能。
3. HSP90B1：
   • 也被称为HSP90β1或grp94，主要存在于内质网中，而不是细胞质。
   • 参与新合成蛋白质的折叠和质量控制，特别是在内质网中。
   • 在抗原呈递、免疫反应以及某些癌症中发挥重要作用。
   • 与HSP90AA1和HSP90AB1相比，HSP90B1在进化上更为古老，具有不同的调控机制。

   总的来说，HSP90家族蛋白在细胞内发挥着关键的分子伴侣功能，对于维持蛋白质稳态和细胞功能至关重要。它们在多种生物学过程和疾病中发挥作用，是重要的研究和靶点。

因此，得到 9个管家基因： c("ACTB", "GAPDH", "HSP90AA1", "HSP90AB1", "HSP90B1", "RPL1A", "TBP", "UBC", "YWHAG")

HSC 是 "Head and eck Squamous Cell Carcinoma" 的缩写，中文全称为头颈鳞状细胞癌。关于TCGA数据的下载以及整理，我们前面已经给大家介绍过两个帖子了：

• TCGA样本barcode详细介绍
• TCGA数据库| 如何将表达矩阵与样本临床数据进行合并？

1、HSC 数据下载

根据上面的帖子，很容易就可以将 HSC 这个癌症的基因表达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-HSC"

# 设置query 
query <- GDCquery(
  project = "TCGA-HSC",  # 癌症类型，选择 TCGA-HSC
  ="Transcriptome Profiling",     # 数据类别, 用getProjectSummary(project)查看所有类别
   ="Gene Expression Quantification", # 数据类型
  ="STAR - Counts"                # 工作流类型
)

## 下载数据
GDCdownload(query=query, files.= 50, directory = "./")


## 整理数据并存储为 R对象
GDCprepare(query,save=T,save.filename=paste0(pro,".transcriptome.Rdata"), directory = "./")

## 再次加载
load("TCGA-HSC/Rdata")
ls()
names(assays(data))
rowdata <- rowData(data)
table(rowdata$gene_type)

mrna <- data[rowdata$gene_type == "protein_coding",]
mrna_count <- assay(mrna,"unstranded")   # mRA的counts矩阵
mrna_fpkm <- assay(mrna,"fpkm_unstrand") # mRA的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-HSC/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-HSC/mrna_fpkm_symbol.rds")

# "HSP90": "HSP90AA1" "HSP90AB1" "HSP90B1"
housekeeper_genes <- c("ACTB", "GAPDH", "HSP90AA1", "HSP90AB1", "HSP90B1", "RPL1A", "TBP", "UBC", "YWHAG")

mrna_count_symbol1[housekeeper_genes, 1:6]
mrna_fpkm_symbol1[housekeeper_genes, 1:6]

2、差异表达分析

随便选一个 差异算法 DESeq2吧：

rm(list = ls()) 
opti(stringsAsFactors = F)
library(ggplot2)
library(ggstatsplot)
library(ggsci)
library(RColorBrewer)
library(patchwork)
library(ggplotify)
library(limma)
library(edgeR)
library(DESeq2)
opti(scipen = 100)

symbol_matrix <- readRDS(file = 'TCGA-HSC/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(=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 <- (res)
save(DEG_deseq2, file= 'TCGA-HSC/DEG_deseq2.Rdata') 

# "HSP90": "HSP90AA1" "HSP90AB1" "HSP90B1"
housekeeper_genes <- c("ACTB", "GAPDH", "HSP90AA1", "HSP90AB1", "HSP90B1", "RPL1A", "TBP", "UBC", "YWHAG")

temp <- DEG_deseq2[housekeeper_genes ,]
temp

火山图如下：

你可以看到，即使使用比较低的阈值进行筛选：FC>1.5, Pvalue<0.05，这些基因基本上都没有差异，除了有一个因为FC在临界点之外。

问题出现在哪里呢？你的分析中 存在 housekeeper基因 发生了显著差异表达吗？