1.不進行merge的pipeline
stringtie官方說明文檔中說明
如果對新的異構(gòu)體不感興趣,可以使用下面簡單的pipeline進行定量,而不進行merge

也就是說如果只對已知的轉(zhuǎn)錄本感興趣就可以不進行merge option
如果要merge則利用下面的pipeline

2.進行merge 之后的geneid問題
merge之后得到stringtie_merge.gtf,再用它作為參考注釋文件進行轉(zhuǎn)錄本重新組裝,為了后續(xù)方便輸入DEseq2,利用官方的prepDE.py提取矩陣。得到兩個矩陣。gene_count_matrix.csv,trancript_count_matrix.csv。這兩個矩陣的gene_id和transcript_id有很多是stringtie內(nèi)部根據(jù)新異構(gòu)體命名如STRG.XXXX和MSTRG.XXXX,單純只看stringtie的結(jié)果文件不太能很好區(qū)分這些異構(gòu)體來自哪個基因,則需要利用gffcomepare對merge后的stringtie_merge.gtf進行分析,來確定已注釋的轉(zhuǎn)錄本和有多少新的轉(zhuǎn)錄本(以及這些轉(zhuǎn)錄本所對應(yīng)的基因)
另外,gene_count_matrix中一些已注釋基因查不到,是因為在count的計算中:如果一個基因有多個轉(zhuǎn)錄本,且該轉(zhuǎn)錄本中含有新的異構(gòu)體,則那個基因的名字變?yōu)閟tringtie內(nèi)部的編號,并且counts數(shù)是所有轉(zhuǎn)錄本相加(不管是不是新的)。(別問我怎么知道的,都是淚)
所以要利用gffcompare中的annotated.gtf文件,對基因進行注釋,tracking文件也行。
以下以gene_count_matrix為例,trancript_count_matrix也類似
library(dplyr)
library(rtracklayer)
#-------------------------------------------------讀入轉(zhuǎn)錄本組裝后的matrix
expr.gene <- read.csv("stringtie_gene_count_matrix.csv")
#-------------------------------------------------讀入annotated.gtf
ensembl_anno <- rtracklayer::import('stringtie_merge.annotated.gtf')
ensembl_anno <- as.data.frame(ensembl_anno)
anno_result <- dplyr::left_join(expr.gene,ensembl_anno[,(11:14)],by ="gene_id")
anno_result <- anno_result[!duplicated(anno_result$gene_id),]
即可得到大部分基因的gene_name(即SYMBOL),之后再用biomaRt等r包進行轉(zhuǎn)換,就可以得到其他ID。
參考:
stringtie
gffcompare