stringtie merge與gene_id

1.不進行merge的pipeline

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


DE_pipeline_refonly.png

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


DE_pipeline.png

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

最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成,瀏覽時請結(jié)合常識與多方信息審慎甄別。
平臺聲明:文章內(nèi)容(如有圖片或視頻亦包括在內(nèi))由作者上傳并發(fā)布,文章內(nèi)容僅代表作者本人觀點,簡書系信息發(fā)布平臺,僅提供信息存儲服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容