IOBR2：转录组数据预处理

文摘 2024-10-09 00:10 北京

前言：组织的转录组的数据主要包含芯片和RNAseq，转录组的芯片主要来源于Affymetrix和Illumina公司，RNAseq几乎是二代测序的数据；在做肿瘤微环境解析和后续的signature score计算时，首先需要对基因表达矩阵进行基因注释和标准化。

注释的目的：因为后续的TME解析工具或者评分计算都是主要识别gene symbol，而无法识别探针或者Ensembl或者entrez的编号；为了简化这个步骤，IOBR提供了anno_eset()；方便进行Affymetrix和illumina芯片的注释，以及针对RNAseq的Ensembl和entrez的注释，也就是将这些ID转化成我们比较熟悉的gene symbol。

标准化目的：实现更精确的样本间可比性；在RNA测序数据分析中，"counts per million"（CPM）和"transcripts per million"（TPM）是两种常用的标准化方法，用于将原始计数数据（reads counts）转换为相对丰度或表达量。这些方法旨在考虑样本之间的总读取量差异以及基因长度的影响。Counts per million (CPM): CPM是将原始计数数据（reads counts）除以每个样本的总读取数量，然后乘以一百万（million）得到的值。CPM反映了基因在样本中的相对表达水平，并考虑到了每个样本的总读取数量。Transcripts per million (TPM): TPM是一种更进一步的标准化方法，除了考虑每个基因的读取数量，还考虑了基因的长度。TPM除了将原始计数数据除以每个样本的总读取数量外，还将结果除以基因的长度（以千碱基数为单位），然后乘以一百万（million）。这种标准化方法更好地考虑了基因长度对表达量的影响，因此在比较不同基因之间的表达水平时更具有可比性。对于RNA测序数据，特别是用于表达量分析时，通常更常用TPM而不是CPM。因为TPM考虑了基因长度的因素，使得不同基因之间的表达量更具有可比性。所以，对于RNAseq数据标准化而言，使用TPM更为合适。

IOBR在此提供了count2tpm()用于RNAseq数据的count to TPM的转化，而芯片的数据大部分能获取到的都是已经标准化的，主要可通过affay和limma包进行处理。因为现在几乎没有人测芯片（不仅贵，而且信息相对匮乏），所以RNAseq的占比会越来越大。

1.加载相关R包

1library(IOBR)
2library(tidyverse)
3library(clusterProfiler)

2.使用R包下载基因表达数据

2.1.GEOquery包下载基因表达数据

下载胃癌数据集 GSE62254（芯片数据）。

1if (!requireNamespace("GEOquery", quietly = TRUE))  BiocManager::install("GEOquery")
2library("GEOquery")
3eset_geo<-getGEO(GEO     = "GSE62254", getGPL  = F, destdir = "./")
4eset    <-eset_geo[[1]]
5eset    <-exprs(eset)   # 获取表达矩阵

1eset[1:5,1:5]  # 查看表达矩阵的格式
2
3##           GSM1523727 GSM1523728 GSM1523729 GSM1523744 GSM1523745
4## 1007_s_at  3.2176645  3.0624323  3.0279131   2.921683  2.8456013
5## 1053_at    2.4050109  2.4394879  2.2442708   2.345916  2.4328582
6## 117_at     1.4933412  1.8067380  1.5959665   1.839822  1.8326058
7## 121_at     2.1965561  2.2812181  2.1865556   2.258599  2.1874363
8## 1255_g_at  0.8698382  0.9502466  0.8125414   1.012860  0.9441993

查看IOBR的不同注释文件，挑选与表达矩阵的基因探针相匹配的注释文件。IOBR内置了多种数据类型的注释文件，包括Affymetrix微阵列数据、Illumina微阵列和RNA-seq数据。

1#人RNA-seq转录组注释文件    
2head(anno_grch38)
3# id eff_length        gc symbol      mgi_id      gene_type     start       end transcript_id  ont
4# 1 ENSMUSG00000000001       3262 0.4350092  Gnai3   MGI:95773 protein_coding 108014596 108053462          <NA> <NA>
5#   2 ENSMUSG00000000003        902 0.3481153   Pbsn MGI:1860484 protein_coding  76881507  76897229          <NA> <NA>
6#   3 ENSMUSG00000000028       3506 0.4962921  Cdc45 MGI:1338073 protein_coding  18599197  18630737          <NA> <NA>
7#   4 ENSMUSG00000000031       2625 0.5588571    H19   MGI:95891         lncRNA 142129262 142131886          <NA> <NA>
8#   5 ENSMUSG00000000037       6397 0.4377052  Scml2 MGI:1340042 protein_coding 159865521 160041209          <NA> <NA>
9#   6 ENSMUSG00000000049       1594 0.5050188   Apoh   MGI:88058 protein_coding 108234180 108305222          <NA> <NA>

1#小鼠RNA-seq转录组注释文件
2head(anno_gc_vm32)
3# id eff_length        gc symbol      mgi_id      gene_type     start       end transcript_id  ont
4# 1 ENSMUSG00000000001       3262 0.4350092  Gnai3   MGI:95773 protein_coding 108014596 108053462          <NA> <NA>
5# 2 ENSMUSG00000000003        902 0.3481153   Pbsn MGI:1860484 protein_coding  76881507  76897229          <NA> <NA>
6# 3 ENSMUSG00000000028       3506 0.4962921  Cdc45 MGI:1338073 protein_coding  18599197  18630737          <NA> <NA>
7# 4 ENSMUSG00000000031       2625 0.5588571    H19   MGI:95891         lncRNA 142129262 142131886          <NA> <NA>
8# 5 ENSMUSG00000000037       6397 0.4377052  Scml2 MGI:1340042 protein_coding 159865521 160041209          <NA> <NA>
9# 6 ENSMUSG00000000049       1594 0.5050188   Apoh   MGI:88058 protein_coding 108234180 108305222          <NA> <NA>

使用anno_hug133plus2文件对基因表达矩阵进行注释，并删除重复基因

1eset <- anno_eset(eset       = eset,
2                  annotation = anno_hug133plus2,
3                  symbol     = "symbol",
4                  probe      = "probe_id",
5                  method     = "mean")

2.2.使用UCSCXenaTools下载RNAseq数据

 1if (!requireNamespace("UCSCXenaTools", quietly = TRUE))  
 2BiocManager::install("UCSCXenaTools")
 3library(UCSCXenaTools)
 4    
 5# 下载GDC的TCGA胃癌数据（count数据）
 6eset_stad<-XenaGenerate(subset = XenaCohorts =="GDC TCGA Stomach Cancer (STAD)") %>% 
 7XenaFilter(filterDatasets    = "TCGA-STAD.htseq_counts.tsv") %>% 
 8XenaQuery() %>%
 9XenaDownload() %>% 
10XenaPrepare()

1head(eset_stad)
2# A tibble: 6 × 408
3#   Ensembl_ID      `TCGA-D7-5577-01A` `TCGA-D7-6818-01A` `TCGA-BR-7958-01A` `TCGA-D7-8572-01A` `TCGA-VQ-A91Z-01A` `TCGA-HU-A4HD-01A` `TCGA-VQ-A8PQ-01A`
4#   <chr>                        <dbl>              <dbl>              <dbl>              <dbl>              <dbl>              <dbl>              <dbl>
5# 1 ENSG0000000000…              11.3               10.6               10.2               10.8               11.8               12.7               10.5 
6# 2 ENSG0000000000…               1                  0                  0                  2                  1                  2.81               2   
7# 3 ENSG0000000041…              11.2                9.96              10.2               11.6               12.3               12.7               11.1 
8# 4 ENSG0000000045…               9.13               8.77               8.73              10.1               10.4               11.1               11.1 
9# 5 ENSG0000000046…               8.42               8.89               8.46               9.87              10.5               10.7                9.29

2.3.标准化和注释

将基因表达矩阵转换为 TPM 格式，并进行后续注释。

1eset_stad$Ensembl_ID<-substring(eset_stad$Ensembl_ID, 1, 15)
2eset_stad<-column_to_rownames(eset_stad, var = "Ensembl_ID")
3eset_stad<-(2^eset_stad)+1
4# 与anno_eset()功能类似，但最终输出的是TPM数据
5eset_stad<-count2tpm(countMat = eset_stad, 
6                     idType   = "Ensembl", 
7                     org      =  "hsa", 
8                     source   = "local" )

1head(eset_stad)[1:6,1:6]
2           TCGA-D7-5577-01A TCGA-D7-6818-01A TCGA-BR-7958-01A TCGA-D7-8572-01A TCGA-VQ-A91Z-01A TCGA-HU-A4HD-01A
3TSPAN6         33.3939485      16.90505921       24.2112184       20.2282326      32.43373343        32.716332
4TNMD            0.1236946       0.06733912        0.1222764        0.1695712       0.08577831         0.118845
5DPM1           15.9689278       5.34679603       11.1486682       16.3134541      23.62654342        15.676932
6SCYL3           4.9602209       3.15520429        5.5719893        8.0435150       8.08128172         6.801386
7C1orf112        3.4965175       3.95405854        5.3206837        7.8398147      10.52599048         5.913285
8FGR            10.6887954      18.30915305       18.1974228       20.8695954      11.66757533         4.091034

3.识别异常样本

这里以ACRG微阵列数据作为示例

1# 异常样本识别的可视化
2res <- find_outlier_samples(eset = eset, 
3                            project = "ACRG", 
4                            show_plot = TRUE ）

异常样本显示

删除异常样本

1eset1 <- eset[, !colnames(eset)%in%res]

4.主成分分析

 1data("pdata_acrg", package = "IOBR")
 2res<- iobr_pca(data       = eset1,# PCA 的输入数据。它应该是一个矩阵或一个数据框。
 3               is.matrix   = TRUE,# 表达谱数据是否为矩阵数据
 4               scale       = TRUE,# 指定是否缩放输入数据。默认为 TRUE。
 5               is.log      = FALSE,# 指定是否对输入数据进行对数转换。默认为 FALSE。
 6               pdata       = pdata_acrg, # 与主成分相关的附加数据。包含样本分组数据。默认为 NULL。
 7               id_pdata    = "ID", #在 'pdata' 中代表与 'data' 匹配的 ID 的列名。默认为 "ID"。
 8               group       = "Subtype",#在 'pdata' 中代表要对点进行着色的组/类别的列名。默认为 NULL。
 9               geom.ind    = "point", #PCA 图中点的几何表示类型。默认为 "point"。
10               cols        = "normal", #用于组类别的颜色方案。默认为 "normal"。
11               palette     = "jama", # 用于组类别的颜色调色板。默认为 "jama"。
12               repel       = FALSE, #指定是否排斥数据点以避免重叠。默认为 FALSE。
13               ncp         = 5, #在 PCA 中保留的维数。默认为 5。
14               axes        = c(1, 2), #要绘制的维度/轴。默认为 c(1, 2)。
15               addEllipses = TRUE) # 指定是否在图中添加椭圆。默认为 TRUE。

表型相关PCA图

5.去除批次效应

去除批次首先可以通过比如PCA降维发现批次，用户可以通过以上的iobr_pca查看目标分类与主成分是否有关联，比如dataset不同、检测技术不同、时间上的批次、送检导致的批次、实验处理上的批次（非干预导致）；假如存在明显的非生物因素或者不合理因素导致的距离，那可能就要尝试去除批次来实现后续的分析和比较。remove_batcheffect（）这个函数可以针对array（芯片）和count（RNAseq）选择各自的校正算法进行校正，需要通过data_type进行指定。

5.1. 微阵列数据批次校正

下载胃癌数据集 GSE57303（芯片数据）

 1eset_geo<-getGEO(GEO     = "GSE57303", 
 2                 getGPL  = F, 
 3                 destdir = "./")
 4eset2    <-eset_geo[[1]]
 5eset2    <-exprs(eset2)
 6eset2<-anno_eset(eset       = eset2,
 7                 annotation = anno_hug133plus2,
 8                 symbol     = "symbol",
 9                 probe      = "probe_id",
10                 method     = "mean")

去除批次效应

 1eset_com <- remove_batcheffect( eset1       = eset1,  
 2                                eset2       = eset2,   
 3                                eset3       = NULL,
 4                                id_type     = "symbol",
 5                                data_type   = "array", 
 6                                cols        = "normal", 
 7                                palette     = "jama", 
 8                                log2        = TRUE, 
 9                                check_eset  = TRUE,
10                                adjust_eset = TRUE,
11                                repel       = FALSE,
12                                path        = "result")

校正前后结果对比

5.2. count数据批次校正

加载IOBR内置数据集"eset_stad"和"eset_blca"。这两个数据是来源于胃癌和膀胱癌的基因表达矩阵（count），可通过ComBat-seq在count层面上进行批次校正。校正后的count会通过count2tpm进行转化，最终实现样本间数据的可比性。

1data("eset_stad", package = "IOBR")
2head(eset_stad)
3# TCGA-BR-6455 TCGA-BR-7196 TCGA-BR-8371 TCGA-BR-8380 TCGA-BR-8592 TCGA-BR-8686 TCGA-BR-A4IV TCGA-BR-A4J4 TCGA-BR-A4J9 TCGA-FP-7916
4# ENSG00000000003         8006         2114          767         1556         2806         2923         1524         7208          711         2747
5# ENSG00000000005            1            0            5            5           60            1           22            2            0            3
6# ENSG00000000419         3831         2600         1729         1760         2273         1934         2838         4418         2426         2824
7# ENSG00000000457         1126          745         1040         1260         1814          707         1683         1335         1590         1672
8# ENSG00000000460          857          463          231          432          635          323          270          423          276          773
9# ENSG00000000938          758         1126          557          557          828          666          760          597          370          688

1data("eset_blca", package = "IOBR")
2head(eset_blca)
3# TCGA-2F-A9KO TCGA-2F-A9KP TCGA-2F-A9KQ TCGA-2F-A9KR TCGA-2F-A9KT
4# ENSG00000000003         6092        11652         5426         4383         3334
5# ENSG00000000005            0            4            1            1            0
6# ENSG00000000419         3072         2656         1983         2061         2930
7# ENSG00000000457         1302          984         1134         1092          496
8# ENSG00000000460          779          924          421          386          318
9# ENSG00000000938          436          116          312          590          362

去除批次效应：

1eset_com <- remove_batcheffect(eset_stad, 
2                               eset_blca, 
3                               id_type = "ensembl", 
4                               data_type = "count")

校正前后结果对比

此部分的英文教程来自：https://iobr.github.io/book/rna-data-preprocessing.html

研究成果正在投稿中，如果用户在发表的研究中使用了IOBR2，请引用以下预印版本：

Zeng D, Fang Y, …, Liao W (2024) IOBR2: Multidimensional Decoding of Tumor Microenvironment for Immuno-Oncology Research. bioRxiv, 2024.01. 13.575484

Wang et al., (2019). The UCSCXenaTools R package: a toolkit for accessing genomics data from UCSC Xena platform, from cancer multi-omics to single-cell RNA-seq. Journal of Open Source Software, 4(40), 1627, https://doi.org/10.21105/joss.01627

Zhang et al., ComBat-seq: batch effect adjustment for RNA-seq count data, NAR Genomics and Bioinformatics, Volume 2, Issue 3, September 2020, lqaa078, https://doi.org/10.1093/nargab/lqaa078

Leek, J. T., et al., (2012). The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics, 28(6), 882-883.

原文链接：https://www.biorxiv.org/content/10.1101/2024.01.13.575484v2.full.pdf

详细教程：https://iobr.github.io/book/

IOBR代码仓库：https://github.com/IOBR/IOBR

反馈问题：https://github.com/IOBR/IOBR/issues

如有任何问题请邮件联系：fyr_nate@163.com 或 interlaken@smu.edu.cn

灵活胖子的科研进步之路

医学博士，R语言及Python爱好者，科研方向为真实世界研究，生信分析与人工智能研究。