差异分析（二）

发表于 更新于 分类于 ， 阅读次数： Waline：

原始数据处理

这篇文章旨在提供一个更加方便准确的差异分析流程，同样以GSE115269为例。

首先是原始数据的下载，GEO 项目页面中的 Series Matrix File(s) 是文章作者所提供处理好的芯片表达数据，往往经过了一些标准化，或者其他的数据格式化的方式，如果你想要比较不同项目之间的样本，除了要进行去批次外还要考虑到每个矩阵数据处理方式不同，故，如果是比较不同项目之间的样本，那么应当从 [GSE115269_RAW.tar]( (http)) 开始处理。当然从RAW开始的处理流程目前也有非常成熟的包。但这不是本教程的核心，这里同上篇文章一样，用的也是文章作者处理好的芯片表达数据。

你可以像上篇文章一样手动在网页下载 表达数据 和 平台信息 (这个是soft格式，只是格式不一样，下其他格式的是一样的)

同样的你也可以直接在R里面用代码直接从网页导入数据，在开始代码前我们先做一些简单的准备工作

# ====准备工作====

# 初始化环境
rm(list = ls())		
pacman::p_unload(pacman::p_loaded(), character.only = TRUE)

# 设置工作目录
setwd(dirname(rstudioapi::getSourceEditorContext()$path))
getwd()

将R中所有环境变量去除，卸除所有已载入的包（因为有许多包含有同名函数），然后设置工作目录

然后直接用 GEOquery 包中的函数获取GEO数据

如果没有安装本文章中的包都可以通过下面的命令安装

# 方式1，安装的是R通过官方审查的包
install.packages('[包名]')

# 方式2，安装的是BiocManage审核过的包，同理也有许多其他的非官方的包管理器，但生信领域用的比较多的是BiocManage
BiocManager::install('[包名]')

# 方式3，未通过任何人审核，发在github上的包
devtools::install_github('[仓库链接]')

现在我默认本文中所有包你已经安装好

1.数据下载

现在我们开始数据下载，先介绍如何用R直接下载数据

library(GEOquery)

# 获取数据
origin_data <- getGEO(
    GEO = 'GSE115269',
    getGPL = TRUE
)

<- 的默认快捷键是 Alt + -

当然因为网络的原因并不是总能成功，如果非要批量下数据，推荐还是用python写爬虫

然后是通过网页下载的方法，之前的文章已经提到过如何下载。这里我直接给出下载链接

GSE115269_series_matrix.txt.gz

GPL18233_family.soft.gz

将这两个文件下载到R代码的同目录下

library(GEOquery)
library(dplyr)		# 这是一个用于数据处理和管处理(逐步处理原始传参)的包 %>% 需要载入dplyr

# 读取原始数据
origin_data <- getGEO(
    filename = './GSE115269_series_matrix.txt.gz',
    getGPL = F
)

# 读取平台信息
GPL18233 <- getGEO(filename = './') %>% Table()

请注意 *.family.soft.gz 文件比较大，事实上，*.annot格式的文件也可以按照上面的代码进行处理，只不过本例平台中没有 *.annont 格式的文件，在 GPL10558 这个平台中

image-20230331195525388

Annotation SOFT table 中可以下载到*.annot格式文件，且比较小

同样的下载txt格式的数据也完全可以，下载按钮是这个image-20230403083323240

# 读取平台信息
GPL18233 <- read.table('./GPL18233-21178.txt', header = TRUE, comment.char = '#', sep = '\t')

2.数据处理

第二步是将数据处理成我们想要的样子，具体的来讲

表达矩阵原来的数据是 行为探针ID，列为样本的矩阵。我们需要把探针ID转化为基因名或者是一些基因的编号

我们先来看我们的平台注释信息，这里里面有 ID（探针ID）信息和基因名 gene_assignment 信息，这里我们只要基因名的信息，其他暂时不用。对这个平台进行一个处理提取出ID和基因名的对应信息

image-20230401195334996

这里我们可以用正则表达式将中间这个基因提取出来表达式为 \b[A-Z][a-z0-9]+ 

\b[A-Z]: 单词开头为A-Z之间任意一个值
[a-z0-9]+: 后续有若干个a-z以及0-9的字符

在R中实现

GPL18233$gene_name <- str_match(GPL18233$gene_assignment, "\\b[A-Z][a-z0-9]+")[, 1]

R中用正则需要将 \ 进行转义，所以用了两个\符号，每种语言的正则表达式有些细微的不同，需要微调一下。最后的 [, 1] 是因为 str_match 函数匹配的结果是一个矩阵，会把不同的组放到不同列中，我们的正则匹配中没有组的信息所以只有一组，之接用[, 1] 就行

现在我们得到了探针ID其基因的信息，但是还有个问题，就是这些基因不一定都是表达基因，有些可能不会表达成蛋白，这些非编码的基因是我们所不需要的，需要去除。

从 genecode 或者 NCBI 上下载小鼠的全基因注释信息，读取到R

# 读取参考基因组
annotation <- read.table('gencode.vM32.annotation.gtf.gz',  header = F,  sep = '\t')

数据格式如下图，其中V9列有用的部分是 gene_type 和 gene_name，我们需要的是所有gene_tpye 为 protein_coding 的 基因名

image-20230403084509750

# 参考基因组进行处理，提取有效信息
annotation$gene_type <- annotation$V9 %>% 
  stringr::str_extract('gene_type.*?;') %>% 
  stringr::str_replace('gene_type ', '') %>% 
  stringr::str_replace(';', '')

annotation$gene_name <- annotation$V9 %>% 
  stringr::str_extract('gene_name .+?;') %>% 
  stringr::str_replace('gene_name ', '') %>% 
  stringr::str_replace(';', '')

annotation <- annotation %>% 
  dplyr::filter(
    gene_type == 'protein_coding',
    V3 == 'transcript'
  )

# 去重
annotation <- annotation[!duplicated(annotation$gene_name), ]

最后的数据格式如下图：

image-20230403085006747

然后对平台数据和注释信息取一个交集，提取出平台中所有protenin_coding的基因和探针ID的对应关系

# 生成新的平台信息(仅含有编码基因)
mRNA_annot <- merge(GPL18233, annotation)

因为此时 GPL18233 和 annotation 只有 gene_name 是相同的列，故会自动按 gene_name 取交集，若有多个相同的列，可以用by参数自行指定

得到了注释信息后，我们可以开始注释芯片表达矩阵了，在注释过程中需要知道一个问题：探针和基因并不是一一对应的，会出现一个探针对应多个基因或者一个基因对应多个探针的情况。这个问题出现的本质是基因是条长的ATGC序列，而探针是一条短的子序列，不同基因有可能子片段有部分ATGC重复，同样的同基因的不同位置也可能有重复所以导致了一探针对多基因，和一基因对多探针。

这两种情况也需要两种处理方式

  1. 一探针对多基因：这种情况为了尽可能多的保留基因的信息我们需要保留每一个探针
  2. 一基因对多探针：这种情况我们只需要一个表达值，可以采用取均值或中值的方式处理。

最后目的是下图数据

image-20230403091720955

处理成

image-20230403091801752

# 注释表达矩阵
expr_origin_data <- exprs(origin_data) %>% 
  as.data.frame() %>% 
  tibble::rownames_to_column('ID') %>% 
  merge(mRNA_annot[, c(1, 2)], .) %>% 
  dplyr::select(-'ID') %>%    # 删除第一列 probe_id
  aggregate(.~ gene_name, data = ., median) %>%  # 以gene_name作为因子水平，相同的取中位数
  tibble::column_to_rownames(var = 'gene_name')

这里可以注意到，表达矩阵中的值分布上下差太大，可以选择进一步处理一下数据，比如取一个log2值，做个标准化等，但是如果取log2值会出现负值，导致差异分析最后的fold change 的正负号表达含义不准确，因此需要整个数据右移到0值前，这里就不处理了，知道就行。不进行标准化出来的fold change会比较大，难看一点

最后一步将样本的注释信息从origin_data中提取出来数据处理阶段就算是完成了

# 样本信息
pd_origin_data <- pData(origin_data) %>% 
  dplyr::select(
    sample = geo_accession,
    group = 'genotype:ch1'
  )

pd_origin_data$group <- pd_origin_data$group %>% 
  stringr::str_replace('DJ-1\\+/\\+', 'case') %>%    # 这里默认正则匹配“+”号有其含义所有需要转义
  stringr::str_replace('DJ-1-/-', 'control') %>% 
  factor(levels = c('case', 'control'), ordered = T) # 设置因子顺序，这步影响不大，加不加都行

这组数据研究的是 DJ-1 这个基因敲除后对小鼠有什么影响。我们如果要做差异表达那么其含义就是找出 DJ-1 这个基因敲除前后小鼠编码基因的变化情况。

3.差异表达

差异表达用的是limma包

library(limma)

limma差异分析的代码基本都是一样的，显示构造比较的矩阵，然后进行limma分析最后 expDiff_result 中的数据就是差异分析的结果

# ====差异表达====
# 构造设计矩阵
design <- model.matrix(~0 + factor(pd_origin_data$group))
colnames(design) <- levels(factor(pd_origin_data$group))
rownames(design) <- rownames(pd_origin_data)

# 模型拟合及分析
contrast_matrix <- makeContrasts("control-case", levels = design)    # 这步影响fold change的正负，+：后面那组高，-: 前面那组高
expDiff_result <- limma::lmFit(expr_origin_data, design) %>%
  contrasts.fit(contrast_matrix) %>%
  eBayes() %>%
  topTable(coef = 1, n = Inf) %>%
  na.omit()

image-20230403102902244

  • logFC： 指的是control 和 case 的倍数变化的log值，makeContrasts 函数影响他的正负，一般是log2
  • AveExpr：所有样本的均值，一般不看，没什么用
  • t：t 值表示基因的显著性，是基于t分布计算的统计量，其公式为 (logFC - 0) / SE，其中 SE 是标准误差。 t 值的绝对值越大，则表示基因在两组样本之间的差异越显著。 t 值还可以用来计算基因的调整 p 值。
  • P.Value：指的是单个结果的显著性
  • adj.P.Value：指定是矫正后的p值，因为每个基因的p值基于一个检验假设，我们如果有n个基因，那么就有n个假设，即使每个假设发生的概率都很小，但是基因数量太多会放大这个错误发生的概率。所以需要矫正