做长链非编码RNA（lncRNA）的数据分析，有一个部分是比较mRNA和lncRNA在染色上的分布密度，做完Hisat2——stringtie流程能够分别拿到mRNA和lncRNA的gtf格式注释文件，那如何根据这两个文件按指定的步长计算基因密度呢？
经过搜索找到了非常方便的工具是tbtools
参考推文
image.png
1 是输入的gtf文件
2 是步长
3 Defined Feature Tag 这个是什么意思暂时没有搞明白 默认的是Guess,如果采用默认我这边会遇到报错，改成exon就可以了
输出文件的路径
最终部分结果
image.png
这里遇到一个问题是不是从小到大依次排列下来的，这个可以后续改
也可以先把自己的gtf文件里的顺序更改一下，使用到的工具是 Tbtools里的 GXF Fix
这里参考
image.png
还找到了一个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")
但是这个结果好像不对，暂时看不太懂这个代码
