从GTF文件中提取TSS上下游1kb的区间，要多少行代码？

欢迎关注”生信修炼手册”!

在ATAC_seq数据分析中，需要绘制reads在TSS位点附近的分布图， 如下所示

左侧为NFR reads在TSS位点两侧的分布图，右侧为单个核小体边界reads在TSS位点两侧的分布图，可以看到，NFR reads在TSS位点两侧有明显的富集趋势。

在上述热图中，每一行代表一个转录本/基因，对于TSS附近区域，换个为等长的bin,比如上图中选取了TSS上下游1kb的区域，那么可以按照100bp划分为等长的窗口，统计每个窗口内的测序深度，然后进行可视化。

要绘制这样的热图，首选要根据基因结构注释文件（通常是GFF或者GTF格式）来获取TSS附近区域的染色体区间信息。TSS表示转录起始位点，本身这个概念是针对基因而言的，但是基因有多个转录本，对应的转录起始位点可能不同，所以在统计TSS时，以转录本为单位进行统计。

本文介绍一种方法python提取TSS区间信息的方法，通过以下两个模块来实现

1. gffutils

2. pybedtools

gffutils用于读取GFF/GTF文件，将所有的信息存在一个sqlite数据库中，方便检索和读取，官网文档链接如下

https://pythonhosted.org/gffutils/index.html

首先需要对gtf文件建立sqlite数据库，用法如下

第一个参数为gtf文件的名称，第二个参数为生成的sqlite db文件名称。创建数据库的过程是比较慢的，但是只需创建一次，以后直接读取这个db文件就可以了。

pybedtools用于区间操作，灵活简便, 官方文档链接如下

https://daler.github.io/pybedtools/

通过pybedtools可以轻松的从bed/gtf文件中提取感兴趣的染色体区间，对于TSS位点以及上下游1kb区间的提取方法如下

上述代码加起来不超过15行，python强大的生态使得我们可以只通过几行代码就实现一个TSS区间提取的功能。小到文件提取，格式转换，大到一个成熟的pipeline, python都可以完美驾驭，而且代码简洁高效，对于从事生信的人而言，python是技能列表中的必备项。

·end·

—如果喜欢，快分享给你的朋友们吧—

原创不易，欢迎收藏，点赞，转发！生信知识浩瀚如海，在生信学习的道路上，让我们一起并肩作战！

本公众号深耕耘生信领域多年，具有丰富的数据分析经验，致力于提供真正有价值的数据分析服务，擅长个性化分析，欢迎有需要的老师和同学前来咨询。

  更多精彩

• KEGG数据库，除了pathway你还知道哪些

• 全网最完整的circos中文教程

• DNA甲基化数据分析专题

• 突变检测数据分析专题

• mRNA数据分析专题

• lncRNA数据分析专题

• circRNA数据分析专题

• miRNA数据分析专题

• 单细胞转录组数据分析专题

• chip_seq数据分析专题

• Hi-C数据分析专题

• HLA数据分析专题

• TCGA肿瘤数据分析专题

• 基因组组装数据分析专题

• CNV数据分析专题

• GWAS数据分析专题
  

  写在最后

转发本文至朋友圈，后台私信截图即可加入生信交流群，和小伙伴一起学习交流。

扫描下方二维码，关注我们，解锁更多精彩内容！

一个只分享干货的

生信公众号