tximport 将 Salmon 定量结果导入 DESeq2

在《Salmon 进行转录本定量》文章介绍了 Salmon 定量。定量后将数据导入 DESeq2 官网有简单的 tximport 代码示例，不过里面缺少如何构建 tx2gene 步骤。这篇文章就写完整的 Salmon 结果导入 DESeq2 过程，新手照着做包会包爽。
转录本与基因注释信息存储在 GTF 注释文件，需要用 GenomicFeatures 的 makeTxDbFromGFF 读取构造成 TxDb 对象，TxDb 对象存储 UTR,CDS 等转录本信息。还有 makeTxDbFromUCSC, makeTxDbFromBiomart 等同类函数从这些数据库下载数据然后构建 TxDb 对象，因为网络容易出问题其实不是很推荐，看每个人条件自行尝试。

library(GenomicFeatures, quietly = TRUE)
library(tximport, quietly = TRUE)
library(DESeq2, quietly = TRUE)

> txd <- makeTxDbFromGFF(file = "/Example/GRCh38_hg38/gencode.v33.annotation.gtf", format = "gtf", dataSource = "gencode.v33.annotation.gtf", organism = "Homo sapiens")
Import genomic features from the file as a GRanges object ... OK
Prepare the 'metadata' data frame ... OK
Make the TxDb object ... OK
Warning message:
In .get_cds_IDX(mcols0$type, mcols0$phase) :
  The "phase" metadata column contains non-NA values for features of type
  stop_codon. This information was ignored.
  
> str(txd)
Reference class 'TxDb' [package "GenomicFeatures"] with 5 fields
 $ conn           :Formal class 'SQLiteConnection' [package "RSQLite"] with 7 slots
  .. ..@ ptr                :<externalptr> 
  .. ..@ dbname             : chr ""
  .. ..@ loadable.extensions: logi TRUE
  .. ..@ flags              : int 70
  .. ..@ vfs                : chr ""
  .. ..@ ref                :<environment: 0x7fd5439f9688> 
  .. ..@ bigint             : chr "integer64"
 $ packageName    : chr(0) 
 $ user_seqlevels : chr [1:25] "chr1" "chr2" "chr3" "chr4" ...
 $ user2seqlevels0: int [1:25] 1 2 3 4 5 6 7 8 9 10 ...
 $ isActiveSeq    : logi [1:25] TRUE TRUE TRUE TRUE TRUE TRUE ...
 and 16 methods, of which 2 are  possibly relevant:
   finalize, initialize

# 保存 TxDb 对象到本地
saveDb(txd, file="/Example/GRCh38_hg38/Gencode_V33_TxDb.sqlite")

其中 dataSource 记载数据来源，可以详尽些也可以不填。
TxDb 对象继承自 AnnotationDb 对象，可以用 AnnotationDbi 方法取得数据。AnnotationDbi 教程见《AnnotationDbi教程(以org.Hs.eg.db为例)》

> keytypes(txd)
[1] "CDSID"    "CDSNAME"  "EXONID"   "EXONNAME" "GENEID"   "TXID"     "TXNAME"
> txTogene <- AnnotationDbi::select(txd, keys=keys(txd, "TXNAME"), keytype="TXNAME", columns=c("TXNAME", "GENEID"))
'select()' returned 1:1 mapping between keys and columns
> head(txTogene)
             TXNAME            GENEID
1 ENST00000456328.2 ENSG00000223972.5
2 ENST00000450305.2 ENSG00000223972.5
3 ENST00000473358.1 ENSG00000243485.5
4 ENST00000469289.1 ENSG00000243485.5
5 ENST00000607096.1 ENSG00000284332.1
6 ENST00000606857.1 ENSG00000268020.3

用这个转录本与基因映射表和 Salmon 结果构建 tximport 对象。

sampleList <- c("UMUC-0-1", "UMUC-0-2", "UMUC-0-3", "UMUC-0-4", "UMUC-1-2", "UMUC-1-3", "UMUC-10-1", "UMUC-10-2", "UMUC-10-3", "UMUC-100-1", "UMUC-100-2", "UMUC-100-3")
fileList <- file.path("/datapool/pengguoyu/RNA/20200714DNJ_UMUC3/Salmon", sampleList, "quant.sf")
names(fileList) <- sampleList
txi <- tximport(fileList, type = "salmon", tx2gene = txTogene)

用于 tximport 的 tx2gene 必须第一列是转录本 ID 第二列是基因 ID, 对列名无要求。
再读入样品信息就可以用 DESeqDataSetFromTximport 函数导入到 DESeq2 了。

sampleGroup <- read.csv("/Example/SampleGroup.txt", header = TRUE, row.names = 1)
dds <- DESeqDataSetFromTximport(txi, colData = sampleGroup, design = ~ Group)