做長鏈非編碼RNA（lncRNA）的數(shù)據(jù)分析，有一個部分是比較mRNA和lncRNA在染色上的分布密度，做完Hisat2——stringtie流程能夠分別拿到mRNA和lncRNA的gtf格式注釋文件，那如何根據(jù)這兩個文件按指定的步長計算基因密度呢？

經(jīng)過搜索找到了非常方便的工具是tbtools

參考推文

image.png

• 1 是輸入的gtf文件

• 2 是步長

• 3 Defined Feature Tag 這個是什么意思暫時沒有搞明白 默認(rèn)的是Guess,如果采用默認(rèn)我這邊會遇到報錯，改成exon就可以了

• 輸出文件的路徑

最終部分結(jié)果

image.png

這里遇到一個問題是不是從小到大依次排列下來的，這個可以后續(xù)改

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

這里參考

image.png

還找到了一個R語言的代碼可以統(tǒng)計基因密度

參考鏈接 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")

但是這個結(jié)果好像不對，暫時看不太懂這個代碼