上下调基因数量不平衡？再探！

生信菜鸟团

发布于 2023-09-09 17:02:39

6750

发布于 2023-09-09 17:02:39

文章被收录于专栏：生信菜鸟团生信菜鸟团

引言

这周曾老师给我分享了一篇文章，TCGA-STAD队列肿瘤样本EBV分型后的差异表达基因出现了上下调数量不平衡，想让我看看是不是样本数量的问题

文章背景：

EBVaGCs vs EBVnGCs

We explored the differences in cell populations between EBVaGC and EBVnGC using a combination of bulk tumor and single cell RNA sequencing data.

We identifified three unique immune cell subpopulations in EBVaGC actively dividing T and B cells and B cells that also expressed T cell features.

这些与EBV相关的GCs（EBVaGCs）在分子、组织病理学和临床上都与EBV阴性的GCs（EBVnGCs）不同。虽然EBVaGCs中的病毒基因参与了癌变过程，但与EBVnGCs相比，病毒蛋白也代表了外源抗原，可以触发增强的免疫应答。

eb病毒（EBV）是一种伽马疱疹病毒，已知可导致黏膜上皮细胞和B淋巴细胞的终身感染，有超过90%的世界人口感染了这种病原体

总的来说，EBV感染占全球人类癌症病例的1.5%，占全球癌症死亡的1.8%

本研究的目的是探索EBVaGCs的TME中独特的细胞亚群，目的是识别在EBVaGCs中出现过多的亚群，并在基因表达水平上描述这些细胞亚群及其对TME的贡献。从Zhang等人获得的单细胞RNA测序（scRNA-seq）GC数据集确定了细胞亚型。对以EBVaGC样本中存在的细胞为主导的所有细胞亚群进行了深入的表征。来自癌症基因组图谱（TCGA）胃癌队列的大规模RNA测序肿瘤数据也被使用，以帮助突出EBV阳性和EBV阴性GCs之间一致差异表达的感兴趣的基因。

我这里转录组专辑还是主要看bulk

关于转录组差异分析中上下调基因数量不平衡现象，我们之前也有篇推文讨论了一些可能的解释

若出现差异表达基因中上下调极度不平衡现象可以从以下几个方面进行检查、进一步分析和解释：

实验操作上对RNA样本进行质控查看RIN值
确定自己的数据来源是微阵列还是测序，是否存在许多新基因没有被微阵列检测到或者测序结果中新基因剔除了
检查自己的上游分析和数据处理，如是否没有去除批次效应或者去除了PCR重复
在差异表达分析中检查过滤策略（根据数据分布），对结果调整fold-change和p值的阈值，查看差异表达基因数量变化
根据自己的科学问题通过不同的方法（如对差异基因进行进一步分析、借助外部数据集）联系生物学背景进行解释

获取数据

TCGA STAD Stomach adenocarcinoma

我之前都是在GDC Data portal下

作者下载：

Level 3 mRNA expression data for the TCGA STAD datasets were obtained from the Broad Genome Data Analysis Center’s Firehose server (https://gdac.broadinstitute.org/, accessed on 2 March 2017), with the data normalized using the RSEM algorithm.

应该是这个：

给了MD5值，看一下下载数据的完整性：

获取TCGA数据比较：

GDC Data portal

Xena

这篇文献使用的第三方处理数据库

https://gdac.broadinstitute.org/

https://broadinstitute.atlassian.net/wiki/spaces/GDAC/pages/844334036/FAQ

利用firehose下载TCGA数据 https://zhuanlan.zhihu.com/p/125401980 “不同于TCGA官网下载的数据，这种方法往往数据更新很慢，已经2020年了，数据居然还停留在2016_01_28。但目测样本数和数据量与官网版本区别不太大” 确实是这样的

关于TCGA（癌症基因组数据集）使用指南 https://www.jianshu.com/p/829c3e311e54

确实，整合起来了已经，比较方便

数据挖掘专题 | GDAC Firehose下载 TCGA 数据(https://mp.weixin.qq.com/s?src=11&timestamp=1689388208&ver=4651&signature=vOeFqv5AQ1Pw72IA4PYb8kI5RkMx0eiNjLN9J2Uzw7VvltfLDyfhagi77CDjZmjctppmxLm76cS5HoGTFMHyrjyDmSEYB5tjiqa4d2kaqgdk1aEvuSEqD0s6eZmn&new=1) Broad GDAC Firehose—TCGA数据分析中心(https://www.omicsclass.com/article/1565) TCGA数据分析在线工具总(https://mp.weixin.qq.com/s?src=11&timestamp=1689386976&ver=4651&signature=jrhxwbSF4VnM7gaTrckmhFfOhit4c6H853AWhHGl7ZQbxPWVirGhLRFN0eTEceJhAQmJGLWshbA53PNvRJUJylfspTWerW56KHGQ9hJkShl-07CU2phQoK3m8x96ok&new=1) TCGA数据库的挖掘工具(https://www.cnblogs.com/modaidai/p/16640759.html)

读取数据

见鬼，为什么明明有文件R却检测不到？--①

修改默认下载文件的名称：

我去，该都不能改 --②

将文件复制到上一级目录下：

发现①②操作都可以正常进行

师姐告诉我会不会是路径太长（深度+文件名长度）

联系她的那篇“走近科学”

这bug有《走近科学》那味儿了

发现无论是将这个文件放到上一级目录还是放到同级文件名较短目录下①②操作都可以进行

所以以后在下载使用TCGA的数据时，这一点也需要注意，TCGA默认下载文件名都比较长

TCGA样本id

https://docs.gdc.cancer.gov/Encyclopedia/pages/TCGA_Barcode/

dat <- read.table("STAD.rnaseqv2__illuminahiseq_rnaseqv2__unc_edu__Level_3__RSEM_genes_normalized__data.data.txt",
                  sep = "\t",
                  header = TRUE)
TCGA_id<- colnames(dat)[-1]
library(tidyverse)
library(stringr)
split_list <- str_split_fixed(TCGA_id,"\\.",7)
apply(split_list,2,table)
# [[1]] project
# TCGA 
# 450 
# [[4]] sample type
# 01A 01B 11A 
# 413   2  35
# 说明有35个normal对照样本 415个tumor样本

# participants
table(split_list[,3]) %>% table()
# 1   2 
# 386  32 
# 说明有32个参与者被收集了两种组织

接下来如何筛选所需EBV分型样本？

寻找路径：

GDC今天2023-7-15在维护
看看https://gdac.broadinstitute.org/ 提供的样本信息，找不到

https://gdac.broadinstitute.org/runs/stddata__latest/samples_report/STAD.html

Google检索“How to classify TCGA-STAD samples as infected with EBV”

Exploration of the Tumor Immune Landscape and Identification of Two Novel Immunotherapy-Related Genes for Epstein-Barr virus-associated Gastric Carcinomavia Integrated Bioinformatics Analysis

这篇文章也没说怎么获取的EBVaGC/EBVnGC TCGA-STAD

但他的实验分组和本文文献数量上差不多，可以看到DEGs也是不平衡的

Classification of gastric cancer by EBV status combined with molecular profiling predicts patient prognosis

A formal molecular classification based on TCGA data in 2014 demarcated EBV-associated tumors from EBV-negative tumors.^4^ In the present study, we integrated EBV infection with TMB and LGI status, yielding a novel four-subtype molecular classification system.

提到2014年对TCGA-STAD tumors进一步分型的研究

In the present study, we integrated EBV infection with TMB and LGI status, yielding a novel four-subtype molecular classification system. This approach indicated a significantly different OS for each subtype of gastric cancer.

这篇文章主要是在此基础上联合其他表型进一步分类鉴定标记基因和相关变异

Comprehensive molecular characterization of gastric adenocarcinoma 2014年对TCGA-STAD tumors进一步分型

https://tcga-data.nci.nih.gov/docs/publications/stad_2014/

链接失效

手动检索 https://gdc.cancer.gov/about-data/publications#/

https://gdc.cancer.gov/about-data/publications/stad_2014

还是跳转到GDC

只能等恢复了再看看

（2023-7-17 gdc 恢复）

https://api.gdc.cancer.gov/data/017f7658-4b39-493e-ad16-e00739a56118

和原文描述对的上

但本文收集的样本 30 EBVa vs 353 EBVn

其它相关文章&资料： https://xenabrowser.net/datapages/?dataset=TCGA-STAD.GDC_phenotype.tsv&host=https%3A%2F%2Fgdc.xenahubs.net&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443 https://gdc-hub.s3.us-east-1.amazonaws.com/download/TCGA-STAD.GDC_phenotype.tsv.gz tcga队列的胃癌里面，病人，走转录组上游，分子分型 https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7050242/ （2022年，我上面这篇分子分型是2014年，本文文献获取的数据是2016年版TCGA）从上游看转录组测序fq文件里面的ebv情况 https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7050242/ （跟前面那篇差不多，都自己走上游然后分析）

没有办法了，太菜了，找不到

看看这篇2014年的TCGA EBV分子分型，28EBVa vs 267EBVn

如果DEGs仍表现出下调基因数量远多于上调那么就可以拿来继续研究

使用本文的差异表达基因筛选标准

Starting with 20,531 genes featured in the TCGA dataset, the log 2-fold change (log2FC) in gene expression between EBVaGCs and EBVnGCs was calculated for every gene. p-Values were calculated with R’s built-in wilcox.test function. **q-Values were then calculated for every p-value using the qvalue package **(version 3.16). Genes with the magnitude of log2FC less than 0.4 or q-values greater than 0.05 for differential gene expression were filtered out.

先看看这个2014年的数据集样本和本数据集样本的关系：

EBV_groups_dat <- readxl::read_xlsx("STAD_Master_Patient_Table_20140207.xlsx")
EBV_groups_dat$group <- ifelse(EBV_groups_dat$`CIMP Category`=="GASTRIC-EBV-CIMP","EBVa","EBVn")
table(EBV_groups_dat$group)
# EBVa EBVn 
# 28  267
head(EBV_groups_dat)
# # A tibble: 6 × 47
# `TCGA barcode` `Lauren Class` `Intestinal Type Subclass` `Signet Ring` `WHO Class`     `Pathologic T`
# <chr>          <chr>          <chr>                      <chr>         <chr>           <chr>   
#   1 TCGA-B7-5816   Diffuse        NA                         0             Poorly_Cohesive T4a   
# 2 TCGA-B7-5818   Intestinal     NA                         NA            NA              T2    
# 3 TCGA-BR-4183   Intestinal     Tubular                    NA            Tubular         T4a   
# 4 TCGA-BR-4184   Intestinal     Tubular                    NA            Tubular         T4a   
# 5 TCGA-BR-4187   Diffuse        NA                         1             Poorly_Cohesive TX    
# 6 TCGA-BR-4188   Diffuse        NA                         1             Poorly_Cohesive TX   

new_TCGA_id <- paste(split_list[,1],split_list[,2],split_list[,3],sep = "-")
head(new_TCGA_id)
# "TCGA-3M-AB46" "TCGA-3M-AB47" "TCGA-B7-5816" "TCGA-B7-5818" "TCGA-B7-A5TI" "TCGA-B7-A5TJ"
colnames(dat) <- c("gene_id",new_TCGA_id)
dat <- dat[-1,]
rownames(dat) <- dat$gene_id
dat$gene_id <- NULL

# 看看2014的样本信息在不在这里
table(EBV_groups_dat$`TCGA barcode` %in% new_TCGA_id)
# FALSE  TRUE 
# 15   280 
# 看起来大部分都在
# 重点看EBVa在不在
table(EBV_groups_dat$`TCGA barcode`[EBV_groups_dat$group=="EBVa"] %in% new_TCGA_id)
# FALSE  TRUE 
# 1    27 
# 太好了大部分都在

可以发现2014年这篇EBV分型的样本，特别是EBVa，在本文文献使用的数据集中

那就可以拿来差异表达分析看看

# 提取对应数据
# 先确保是肿瘤样本的
sub_dat <- dat[,ifelse(split_list[,4]=="11A",FALSE,TRUE)]
dim(sub_dat)  # 20531   415
table(duplicated(colnames(sub_dat)))
# FALSE 
# 415
# 获取2014EBV分型样本
EBVa <- EBV_groups_dat$`TCGA barcode`[EBV_groups_dat$group=="EBVa"]
EBVn <- EBV_groups_dat$`TCGA barcode`[EBV_groups_dat$group=="EBVn"]
EBVa_dat <- sub_dat[,which(colnames(sub_dat) %in% EBVa)]  # 27
EBVn_dat <- sub_dat[,which(colnames(sub_dat) %in% EBVn)]  # 250
EBV_dat <- cbind(EBVa_dat,EBVn_dat)
# 发现读进来的表达量值为character
# numeric化
EBV_dat_num <- apply(EBV_dat, 2, as.numeric)
rownames(EBV_dat_num) <- rownames(EBV_dat)

######差异表达分析######

# 作者直接使用的R’s built-in wilcox.test function
# Wilcoxon秩和检验进行非参数检验

# 输入EBV_dat_num 和d
d <- 27
n <- length(colnames(EBV_dat_num))
deg_pvalue <- apply(EBV_dat_num, 1, function(x){
  # wilcoxon
  # print(x)
  res <- wilcox.test(x[1:d],x[(d+1):n],
                     exact = FALSE)
  return(res$p.value)
})
deg_log2fc <- apply(EBV_dat_num, 1, function(x){
  fc <- mean(x[1:d]) / mean(x[(d+1):n])
  return(log2(fc))
})
# > table( qvalue(deg_pvalue)$qvalues < 0.05)
# FALSE  TRUE 
# 12235  8039 
deg_qvalue <- qvalue(deg_pvalue)
# 放一起
dat_to_draw <- data.frame(deg_log2fc,deg_qvalue$qvalues)
# 去掉全部是0的基因
dat_to_draw <-  na.omit(dat_to_draw)

dat_to_draw$regulation <- ifelse(
  dat_to_draw$deg_log2fc<(-0.4) & dat_to_draw$deg_qvalue.qvalues<0.05, "down",
  ifelse(dat_to_draw$deg_log2fc>0.4 & dat_to_draw$deg_qvalue.qvalues<0.05,"up","normal")
)
table(dat_to_draw$regulation)
# down normal     up 
# 4416  14518   1340 

#####绘图######

# 火山图
library(EnhancedVolcano)
EnhancedVolcano(dat_to_draw,x='deg_log2fc',y='deg_qvalue.qvalues',
                lab=rownames(dat_to_draw),
                pCutoff = 0.05,
                FCcutoff = 0.4)

volcano_p <- ggplot(data=dat_to_draw, aes(x=deg_log2fc, y=-log10(deg_qvalue.qvalues),color=regulation)) +  # 按regulated分组颜色
  geom_point(alpha=0.5, size=1.8) +  # alpha控制点的透明度 size大小
  theme_set(theme_set(theme_bw(base_size=20))) + 
  xlab("log2FC") + ylab("-log10(Pvalue)") +
  scale_colour_manual(values = c('blue','black','red'))+theme_bw()
volcano_p

DEGs仍表现出下调基因数量远多于上调

所以就拿这个数据集来做好了

在取子集的时候如果通过索引来取需要注意：

> length(EBV_dat_num[1,][d+1:n])
[1] 277
> length(EBV_dat_num[1,][(d+1):n])
[1] 250

# 整理到函数里
degs_analysis_draw <- function(mat,d){
  print(dim(mat))
  # 输入EBV_dat_num 和d
  n <- length(colnames(mat))
  deg_pvalue <- apply(mat, 1, function(x){
    # wilcoxon
    # print(x)
    res <- wilcox.test(x[1:d],x[(d+1):n],
                       exact = FALSE)
    return(res$p.value)
  })
  deg_log2fc <- apply(mat, 1, function(x){
    fc <- mean(x[1:d]) / mean(x[(d+1):n])
    return(log2(fc))
  })
  # > table( qvalue(deg_pvalue)$qvalues < 0.05)
  # FALSE  TRUE 
  # 12235  8039 
  deg_qvalue <- qvalue(deg_pvalue)
  # 放一起
  dat_to_draw <- data.frame(deg_log2fc,deg_qvalue$qvalues)
  # 去掉全部是0的基因
  dat_to_draw <-  na.omit(dat_to_draw)
  
  dat_to_draw$regulation <- ifelse(
    dat_to_draw$deg_log2fc<(-0.4) & dat_to_draw$deg_qvalue.qvalues<0.05, "down",
    ifelse(dat_to_draw$deg_log2fc>0.4 & dat_to_draw$deg_qvalue.qvalues<0.05,"up","normal")
  )
  print(table(dat_to_draw$regulation))
  
  #####绘图######
  
  # 火山图
  library(EnhancedVolcano)
  EnhancedVolcano(dat_to_draw,x='deg_log2fc',y='deg_qvalue.qvalues',
                  lab=rownames(dat_to_draw),
                  pCutoff = 0.05,
                  FCcutoff = 0.4)
  
  volcano_p <- ggplot(data=dat_to_draw, aes(x=deg_log2fc, y=-log10(deg_qvalue.qvalues),color=regulation)) +  # 按regulated分组颜色
    geom_point(alpha=0.5, size=1.8) +  # alpha控制点的透明度 size大小
    theme_set(theme_set(theme_bw(base_size=20))) + 
    xlab("log2FC") + ylab("-log10(Pvalue)") +
    scale_colour_manual(values = c('blue','black','red'))+theme_bw()
  volcano_p
  
  print("----------end-------------")
  
}

抽样测试

sample size

#########抽样测试########
# 先看看函数能不能用
degs_analysis_draw(mat = EBV_dat_num,d = 27)
# 27 vs 250
# sample size
index <- seq(28,277)
for (sample_size in seq(50,250,by=50)){
  # EBVn 50 100 .. 200 250
  print(sample_size)
  cur_mat <- EBV_dat_num[,c(seq(1,27),sample(index, sample_size))]
  degs_analysis_draw(cur_mat, d)
  # break
}

[1] 50
[1] 20531    77

  down normal     up 
  3823  15037   1357 
[1] "----------end-------------"
[1] 100
[1] 20531   127

  down normal     up 
  4005  14826   1411 
[1] "----------end-------------"
[1] 150
[1] 20531   177

  down normal     up 
  4400  14501   1355 
[1] "----------end-------------"
[1] 200
[1] 20531   227

  down normal     up 
  4320  14553   1387 
[1] "----------end-------------"
[1] 250
[1] 20531   277

  down normal     up 
  4416  14518   1340 
[1] "----------end-------------"

# 散点图
num_degs <- data.frame(
  down = c(3823,4005,4400,4320,4416),
  up = c(1357,1411,1355,1387,1340)
)
# basic scatterplot
scatter_plot <- ggplot(num_degs, aes(x=down, y=up)) + 
  geom_point()
scatter_plot

除了第一个EBVn sample size为50，其他更大的sample size的结果中，DEGs总数保持一个相对稳定的状态

sample distribution

# sample distribution
index <- seq(28,277)
for (test in seq(1,5)){
  # EBVn 50 50 .. 50 50
  # 选小一点，这样可以把样本异质性凸显出来
  # print(sample_size)
  sample_size <- 50
  cur_mat <- EBV_dat_num[,c(seq(1,27),sample(index, sample_size))]
  degs_analysis_draw(cur_mat, d)
  # break
}
# 散点图
num_degs <- data.frame(
  down = c(2939,2932,2996,3609,3396),
  up = c(1199,1197,1217,1329,1234)
)
# basic scatterplot
scatter_plot <- ggplot(num_degs, aes(x=down, y=up)) + 
  geom_point()
scatter_plot