认识TCGA的XENA的转录组测序表达量矩阵

认识TCGA的XENA的转录组测序表达量矩阵 by 生信菜鸟团

转录组测序表达量矩阵（TCGA的XENA）解析 htseq_counts.tsv.gz，6万多行是因为gtf文件，ensembl的ID，最后五行并不是基因。

如果你直接搜索TCGA的XENA，只能是看到它的主页：https://xena.ucsc.edu/，它上面比较好用的是一个网页工具。。。

但是对我们生信工程师来说，最方便的其实是直接下载其各种ngs组学数据，所以前往：https://xenabrowser.net/datapages ，这里我们讲解的是 GDC TCGA 相关数据库文件，任意选择一个进入即可：

GDC TCGA Acute Myeloid Leukemia (LAML) (15 datasets)
GDC TCGA Adrenocortical Cancer (ACC) (14 datasets)
GDC TCGA Bile Duct Cancer (CHOL) (14 datasets)
GDC TCGA Bladder Cancer (BLCA) (14 datasets)
GDC TCGA Breast Cancer (BRCA) (20 datasets)

比如我们选择：cohort: GDC TCGA Breast Cancer (BRCA)

生物信息学上游流程

今天我们要介绍的就是gene expression RNAseq ，它是 HTSeq - Counts (n=1,217) GDC Hub

在 https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/ 可以看到其转录组测序表达量矩阵来源的生物信息学流程，主要是比对和定量。

其中比对的过程就顺便走了融合基因流程，流程图是：

感兴趣的可以去阅读转录组数据分析实战目录，如下所示：

这里需要注意的是定量过程使用的是 htseq-count 这个软件：

htseq-count \
-m intersection-nonempty \
-i gene_id \
-r pos \
-s no \
- gencode.v22.annotation.gtf

而且需要注意是 gencode.v22.annotation.gtf 版本的参考基因组注释文件哦。

表达量矩阵文件

可以看到这个矩阵文件（TCGA-BRCA.htseq_counts.tsv.gz）其实有点古老了：

cohort   GDC TCGA Breast Cancer (BRCA)
dataset ID   TCGA-BRCA.htseq_counts.tsv
download  https://gdc-hub.s3.us-east-1.amazonaws.com/download/TCGA-BRCA.htseq_counts.tsv.gz ; 
samples  1217
version  07-18-2019

这个链接很方便下载，https://gdc-hub.s3.us-east-1.amazonaws.com/download/TCGA-BRCA.htseq_counts.tsv.gz

其主页介绍的很清楚，这个矩阵列是基因名字，行是每个病人的id

矩阵

这个时候的基因名字是ensembl数据库的ID，病人的id是tcga自己命名系统。

ensembl数据库的ID

一般来说我们定量使用的参考基因组注释文件（比如 gencode.v22.annotation.gtf 版本），里面的gene_id也就是Ensembl ID，而gene_name就是我们的gene symbol ID，这些id都有自己的规则系统。

在Ensembl数据库中，每一个转录本（或是一条成熟的mRNA）有一个编码，即 transcript ID。

人类的转录本都由“ENST”加上后面的唯一、固定的11位数字组成，例如：“ENST00000369985”；
其它物种则会在开头字母ENST中插入其它三个字母，如：小鼠的转录本ID由“ENSMUST”开头；
部分 transcript ID的数字后还有小数点，这表示这个ID在数据库中做了几次变更。参考：

Ensemble 一共提供四个数据库服务器访问地址：

ensembldb.ensembl.org：欧洲服务器，只有该服务器可访问 GRCh37 数据集
useastdb.ensembl.org：美洲服务器
asiadb.ensembl.org：亚洲服务器

https://asia.ensembl.org/Help/View?id=151

tcga的病人id

TCGA官网有样本命名的规则，详见：https://docs.gdc.cancer.gov/Encyclopedia/pages/TCGA_Barcode/ ，看完你就明白了 -01A 和 -11A 的区别，虽然很复杂，其实绝大部分小伙伴只需懂Sample列信息即可；

这个矩阵数值并不是整数

细心地小伙伴可能发现了，它矩阵里面的表达量并不是整数，但是我们的转录组定量理论上是数每个基因上面的reads的数量，不应该是小数。所以需要仔细看官方文档，里面写的是：log2(x+1) transformed.，所以我们读取的时候需要注意，得转换：

## dataset: gene expression RNAseq - HTSeq - Counts 
## unit: log2(count+1)
a <- fread( 'TCGA-LAML.htseq_counts.tsv.gz' ,
           data.table = F) #全部是symbol
head(a[ ,1:4])
tail(a[ ,1:4])

mat <- a[1:60483,]
mat = a
rownames(mat) <- a$Ensembl_ID  # keep point
mat[1:4,1:4]
mat <- mat[,-1]
mat[1:4,1:4]

## tidy data
org_dat <- floor(2^mat - 1 ) #floor 返回小于输入值的最大整数值，即输出值将不大于输入值
org_dat[1:4,1:4] 

## filter data
keep_feature <- rowSums (org_dat > 2) > 10
table(keep_feature)

mat <- org_dat[keep_feature, ]
mat <- mat[, colSums(mat) > 1000000]  # !!note: can change to get more genes!!
mat[1:4,1:4]
dim(mat) 
head(colnames(mat))

而且是 HTSeq - Counts软件定量，这个表达量矩阵它最后的五行其实是需要删除的：

> head(a[ ,1:4])
          Ensembl_ID TCGA-AB-2949-03B TCGA-AB-2918-03A TCGA-AB-2943-03A
1 ENSG00000000003.13         5.129283         3.700440         5.209453
2  ENSG00000000005.5         1.000000         1.000000         0.000000
3 ENSG00000000419.11         9.972980         9.885696         9.868823
4 ENSG00000000457.12         9.980140        10.052568        10.965784
5 ENSG00000000460.15         9.517669         9.823367        10.388017
6 ENSG00000000938.11        11.398744        11.913263        11.738937
> tail(a[ ,1:4])
                  Ensembl_ID TCGA-AB-2949-03B TCGA-AB-2918-03A TCGA-AB-2943-03A
60483      ENSGR0000281849.1          0.00000          0.00000          0.00000
60484           __no_feature         24.67842         24.44323         24.49330
60485            __ambiguous         20.57133         20.50281         20.39089
60486        __too_low_aQual          0.00000          0.00000          0.00000
60487          __not_aligned          0.00000          0.00000          0.00000
60488 __alignment_not_unique         24.61085         24.30559         24.75108

转换成功之后是：

> mat[1:4,1:4]
                   TCGA-AB-2949-03B TCGA-AB-2918-03A TCGA-AB-2943-03A TCGA-AB-2851-03A
ENSG00000000003.13               34               12               36                3
ENSG00000000419.11             1004              944              934              630
ENSG00000000457.12             1009             1060             1999              705
ENSG00000000460.15              731              904             1339              317

人类基因通常是2万个但是gtf文件里面是6万多个基因

就需要解析前面的定量过程使用的参考基因组注释文件（比如 gencode.v22.annotation.gtf 版本），下载链接是：https://gdc-hub.s3.us-east-1.amazonaws.com/download/gencode.v22.annotation.gene.probeMap

b <- fread( 'gencode.v22.annotation.gene.probeMap' ,header = T, data.table = F)
head(b)  # ensg has point
> dim(b)  
[1] 60483     6

可以看到确实是6万多个基因，而且里面就有ensembl的数据库的基因的id转换成为基因的名字的信息：

> head(b)  # ensg has point
                 id         gene chrom chromStart chromEnd strand
1 ENSG00000223972.5      DDX11L1  chr1      11869    14409      +
2 ENSG00000227232.5       WASH7P  chr1      14404    29570      -
3 ENSG00000278267.1    MIR6859-3  chr1      17369    17436      -
4 ENSG00000243485.3 RP11-34P13.3  chr1      29554    31109      +
5 ENSG00000274890.1    MIR1302-9  chr1      30366    30503      +
6 ENSG00000237613.2      FAM138A  chr1      34554    36081      -

我们现在看看这些基因的属性。

library(AnnoProbe)
tmp= annoGene(b$gene,'SYMBOL')
tail(sort(table(tmp$biotypes)))

可以看到，因为我们使用annoGene函数的时候最新版gtf文件，转换这个v22版本会一半失败：

> tmp= annoGene(b$gene,'SYMBOL')
Warning message:
In annoGene(b$gene, "SYMBOL") :
  40.03% of input IDs are fail to annotate...

6万多个基因里面，留下来的绝大部分都是蛋白质编码基因，他们的基因名字很少变化：

> tail(sort(table(tmp$biotypes)))

unprocessed_pseudogene 
                  1352 
                 miRNA 
                  1674 
                 snRNA 
                  1789 
                lncRNA 
                  3818 
  processed_pseudogene 
                  4981 
        protein_coding 
                 18334

非编码基因他们的名字是一直在变化的，但是他们的ensembl数据库的id还是可以追踪的，这就是为什么大家并不喜欢使用常规的基因的名字来进行表达量矩阵文件保存。

ixxmu / mp_duty