做长链非编码RNA（lncRNA）的数据分析，有一个部分是比较mRNA和lncRNA在染色上的分布密度，做完Hisat2——stringtie流程能够分别拿到mRNA和lncRNA的gtf格式注释文件，那如何根据这两个文件按指定的步长计算基因密度呢？

经过搜索找到了非常方便的工具是tbtools

参考推文 TBtools | 全基因组 - 基因密度统计，充实你的图片

• 1 是输入的gtf文件

• 2 是步长

• 3 Defined Feature Tag 这个是什么意思暂时没有搞明白 默认的是Guess,如果采用默认我这边会遇到报错，改成exon就可以了

• 输出文件的路径

最终部分结果

这里遇到一个问题是不是从小到大依次排列下来的，这个可以后续改

也可以先把自己的gtf文件里的顺序更改一下，使用到的工具是 Tbtools里的 GXF Fix

这里参考 完美 | GXF Fix 修复 / 优化基因结构注释信息文件 - GTF/GFF3

还找到了一个R语言的代码可以统计基因密度

参考链接 https://davetang.org/muse/2017/08/04/read-gtf-file-r/ https://www.biostars.org/p/169171/

代码

library(dplyr)
library(rtracklayer)
my_obj<-import("gene_density/i.gtf")
my_obj
my_obj@seqinfo@seqlengths<-read.table("gene_density/chr_len.txt",header=F,sep="\t")$V2

my_obj@seqnames@values
my_obj@seqnames@lengths


seqinfo(my_obj)
pome.windows = tileGenome(seqinfo(my_obj), tilewidth=100000, cut.last.tile.in.chrom=T)
pome.windows$totalgenes<-countOverlaps(pome.windows,my_obj)
data.frame(pome.windows) %>% 
  write.table(file = "gene_density/1.tsv",quote=F,row.names = F,sep="\t")

pome.windows = tileGenome(seqinfo(reduce(my_obj)), tilewidth=100000, cut.last.tile.in.chrom=T)
pome.windows$totalgenes<-countOverlaps(pome.windows,reduce(my_obj))
data.frame(pome.windows) %>% 
  write.table(file = "gene_density/1-1.tsv",quote=F,row.names = F,sep="\t")

但是这个结果好像不对，暂时看不太懂这个代码

欢迎大家关注我的公众号

小明的数据分析笔记本

小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！