Kallisto定量结果：如何将transcript-level表达矩阵转为gene-level

生信技能树

发布于 2025-02-12 15:34:43

5100

代码可运行

文章被收录于专栏：生信技能树

运行总次数：0

代码可运行

前面关于Kallisto写了两个笔记，主要是： Kallisto 软件的上、下游分析流程。但是曾老板说他没有看到他想让我做的，嗯：就是Kallisto的定量结果是转录本水平的定量，他想让我学习怎么转换为基因水平的定量。转录本水平的定量研究也是有很多的，我老以为他让我看转录本水平定量的count是小数怎么处理，为什么是小数以及怎么给转录本ID添加转录本symbol和基因symbol。就是我没有get到他的意思！这就来看看~

前面的两个Kallisto笔记：

转换数据准备

这里将transcript-level表达矩阵转为gene-level的矩阵，使用的是R包 tximport，学习链接：https://bioconductor.org/packages/release/bioc/vignettes/tximport/inst/doc/tximport.html。这个软件可以接受这些软件的转录本定量结果：

Salmon (Patro et al. 2017)
Alevin (Srivastava et al. 2019)
Sailfish (Patro, Mount, and Kingsford 2014)
kallisto (Bray et al. 2016)
RSEM (Li and Dewey 2011)
StringTie (Pertea et al. 2015)

还需要一个文件为：tx2gene数据框，包含两列，列名并不重要，但必须第一列是transcript ID，第二列是gene ID这种列的顺序。转录本 ID 必须与在abundance.tsv文件中使用的相同。

第一列：transcript ID
第二列：gene ID

有两种方法可以得到这个对应关系，第一种是定量的时候配套版本的gtf文件，见帖子：在什么情况下基因ID转换会100%失败？

### gtf 
library(tidyverse)
library(rtracklayer)

# 读取gtf文件
gtf <- rtracklayer::import('GSE270697/Mus_musculus.GRCm39.113.gtf.gz') 
# 转为数据框，简单探索
gtf_df <- as.data.frame(gtf)
colnames(gtf_df)

# 过滤transcript行
gtf_df <- filter(gtf_df, type=="transcript")
head(gtf_df)
colnames(gtf_df)

gtf_df_trans <- gtf_df[,c("transcript_id","gene_id")]
head(gtf_df_trans)

第二种是定量的时候使用的参考序列：/nas1/zhangj/database/genome/Mus_musculus/release113/Mus_musculus.GRCm39.cdna.all.fa

使用bash命令得到：

# 转录本id 和基因id的关系
grep '^>' /nas1/zhangj/database/genome/Mus_musculus/release113/Mus_musculus.GRCm39.cdna.all.fa  |cut -d' ' -f 1,4 |sed -e 's/>//g' -e 's/gene://g'  >GRCm39.tx2gene.txt

# gene ID 和name也保存一份
grep '^>' /nas1/zhangj/database/genome/Mus_musculus/release113/Mus_musculus.GRCm39.cdna.all.fa  |cut -d' ' -f 4,7 |sed -e 's/gene_symbol://g' -e 's/gene://g' >GRCm39.g2name.txt

GRCm39.tx2gene.txt：

ENSMUST00000103301.3 ENSMUSG00000076500.3
ENSMUST00000166255.2 ENSMUSG00000090395.2
ENSMUST00000095364.3 ENSMUSG00000090765.2
ENSMUST00000103304.3 ENSMUSG00000094491.3
ENSMUST00000103305.2 ENSMUSG00000096580.2
ENSMUST00000103306.2 ENSMUSG00000076505.2
ENSMUST00000200586.2 ENSMUSG00000076508.4
ENSMUST00000103309.3 ENSMUSG00000076508.4
ENSMUST00000103310.3 ENSMUSG00000094345.3
ENSMUST00000196768.2 ENSMUSG00000096632.3

GRCm39.g2name.txt：

ENSMUSG00000076500.3 Gm20730
ENSMUSG00000090395.2 Gm54608
ENSMUSG00000090765.2 Gm54637
ENSMUSG00000094491.3 Igkv1-133
ENSMUSG00000096580.2 Igkv1-132
ENSMUSG00000076505.2 Igkv1-131
ENSMUSG00000076508.4 Igkv17-127
ENSMUSG00000076508.4 Igkv17-127
ENSMUSG00000094345.3 Igkv14-126
ENSMUSG00000096632.3 Igkv9-124

万事俱备，转换

rm(list=ls())
library(tximport)
library(readr)
library(stringr)

# 定量结果
dir <- file.path('PRJNA1128028/rawdata/kallisto_quant/')
dir
files <- list.files(pattern="*tsv",dir,recursive=T)
files
files <- file.path(dir,files)
files
all(file.exists(files))

转换：

# tx2gene
tx2gene <- read.table(file = "GRCm39.tx2gene.txt",stringsAsFactors = F )
head(tx2gene)


# 导入转换
txi <- tximport(files, type = "kallisto", tx2gene = tx2gene)
names(txi) 
head(txi$length)
head(txi$counts)
str(txi)


# 添加样本名字
files
sapply(strsplit(files,'\\/'), function(x) x[5]) # 看清楚在哪个字段
temp <- sapply(strsplit(files,'\\/'), function(x) x[5]) 
temp <- gsub(".quant", "", temp)
temp
colnames(txi$counts) <- temp

# 提取基因水平count，小数转为整数
tmp <- txi$counts
head(tmp)
exprSet <- apply(tmp,2,as.integer)
rownames(exprSet) <- rownames(tmp) 
dim(exprSet)
exprSet[1:4,1:4]

顺便对基因id做一下symbol转换：

# 基因id转换
id2name <- read.table("GRCm39.g2name.txt",sep = " ", quote = "\"", fill = T)
id2name <- unique(id2name)
id2name <- id2name[id2name$V2!="", ]
head(id2name)

mat <- exprSet %>% 
  as.data.frame() %>% 
  rownames_to_column("ID")
head(mat)

mat_symbol <- merge(id2name, mat, by.x="V1", by.y="ID")
mat_symbol1 <- limma::avereps(mat_symbol[,-c(1,2)], ID=mat_symbol$V2)
dim(mat_symbol1)
mat_symbol1[1:5, ]
keep <- rowSums(mat_symbol1)>0;table(keep)
mat_symbol1 <- mat_symbol1[keep, ]
mat_symbol1[1:5, ]