stringtie merge与gene_id

1.不进行merge的pipeline

stringtie官方说明文档中说明
如果对新的异构体不感兴趣,可以使用下面简单的pipeline进行定量,而不进行merge


DE_pipeline_refonly.png

也就是说如果只对已知的转录本感兴趣就可以不进行merge option
如果要merge则利用下面的pipeline


DE_pipeline.png

2.进行merge 之后的geneid问题

merge之后得到stringtie_merge.gtf,再用它作为参考注释文件进行转录本重新组装,为了后续方便输入DEseq2,利用官方的prepDE.py提取矩阵。得到两个矩阵。gene_count_matrix.csv,trancript_count_matrix.csv。这两个矩阵的gene_id和transcript_id有很多是stringtie内部根据新异构体命名如STRG.XXXX和MSTRG.XXXX,单纯只看stringtie的结果文件不太能很好区分这些异构体来自哪个基因,则需要利用gffcomepare对merge后的stringtie_merge.gtf进行分析,来确定已注释的转录本和有多少新的转录本(以及这些转录本所对应的基因)

另外,gene_count_matrix中一些已注释基因查不到,是因为在count的计算中:如果一个基因有多个转录本,且该转录本中含有新的异构体,则那个基因的名字变为stringtie内部的编号,并且counts数是所有转录本相加(不管是不是新的)。(别问我怎么知道的,都是泪)
所以要利用gffcompare中的annotated.gtf文件,对基因进行注释,tracking文件也行。
以下以gene_count_matrix为例,trancript_count_matrix也类似


library(dplyr)
library(rtracklayer)
#-------------------------------------------------读入转录本组装后的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包进行转换,就可以得到其他ID。

参考:
stringtie
gffcompare

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

推荐阅读更多精彩内容