从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是技能列表中的必备项。