使用Tbtools根据gtf文件统计基因密度

用户7010445

发布于 2021-07-30 10:41:27

1.4K00

代码可运行

运行总次数：0

代码可运行

做长链非编码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")

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

也计划自己写脚本统计了，但是这个逻辑自己还没想明白

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2021-07-10，如有侵权请联系 cloudcommunity@tencent.com 删除

https