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

学术   2025-01-17 18:32   广东  

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

人类 housekeeper基因有哪些

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

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

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

  1. ACTB:编码β-肌动蛋白(Beta-actin),是细胞骨架的重要组成部分。
  2. GAPDH:编码甘油醛-3-磷酸脱氢酶(Glyceraldehyde-3-phosphate dehydrogenase),在糖酵解过程中起关键作用。
  3. HSP90:编码热休克蛋白90(Heat shock protein 90),参与蛋白质折叠和稳定。
  4. RPL13A:编码核糖体蛋白L13a,是蛋白质合成的必需成分。
  5. TBP:编码TATA盒结合蛋白(TATA box-binding protein),在转录起始中起关键作用。
  6. UBC:编码泛素C(Ubiquitin C),参与蛋白质降解和调控。
  7. 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是一种高度保守的分子伴侣蛋白,对于维持蛋白质稳态、帮助蛋白质折叠、防止蛋白质聚集以及调控蛋白质降解等过程至关重要。下面是对这三个基因的具体解读:

  1. 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等数据库挖掘感兴趣,但又不知道如何入门,也许你可以关注一下下面的课程:


    生信技能树
    生物信息学学习资料分析,常见数据格式及公共数据库资料分享。常见分析软件及流程,基因检测及癌症相关动态。
     最新文章