获得转录本id和基因id的对应表格

使用kallisto软件获得转录本的表达量，利用tximport这个R包把所有样本的表达量合并到一起，可以获得转录本的表达量，如果提供转录本id和基因id的对应表格，也可以获得基因的表达量

我是用stringtie这个软件获得gtf文件，然后利用这个gtf文件获得转录本的fasta文件，然后用这个fasta文件去计算表达量

有了stringtie这个软件获得的gtf文件很容易就可以获得转录本和基因id的对应表格

首先使用gffread软件对gtf文件进行操作

gffread -E --keep-genes input.gtf -o input/output.gtf

这个output.gtf文件基本的内容

接下来在R语言里操作

library(readr)
library(stringr)
library(tidyverse)
df <- read_tsv("output.gtf",col_names = FALSE,comment = "#")
f%>%filter(X3=='transcript')%>%select(X9)%>%mutate(X9=str_replace_all(X9,"ID=|Parent=",""))%>%rename("TXNAME;GENEID"="X9")%>%write_csv
    (file = "tx2gene.csv")

然后用tximport读取kallisto的结果

help(package="tximport")
list.files("D:/Bioinformatics_Intro/pomeRTD/kallisto.output",
           recursive = TRUE,
           full.names = TRUE,
           pattern = "*.h5") -> files
library(stringr)
library(tximport)
library(readr)
library(tidyverse)

names(files)<- str_extract(files,pattern = "PRJ[A-z0-9]+/SRR[0-9]+") %>% 
  str_replace("/","_")

tx2gene<-read_delim("D:/Bioinformatics_Intro/pomeRTD/kallisto.output/pomeRTD_tx2gene.csv",
                    delim = ";")
head(tx2gene)

txi.ka.gene<-tximport(files,
                       type = "kallisto",
                       tx2gene = tx2gene)
txi.ka.tx<-tximport(files,
                      type = "kallisto",
                      txOut = TRUE)
txi.ka.gene$abundance %>% dim()
txi.ka.tx$abundance %>% dim()