python 学习之 提取 Ensembl,Gencode 和 Ucsc 基因 TSS 位点

此事古难全,但愿人长久

1引言

Chip-seq 和 ATAC-seq 等经常需要绘制离 TSS 位点处的结合强度:

绘图工具也很多,例如 deeptools 软件,此外我们需要提供基因的 TSS(转录起始位点)的坐标信息才能进行绘制, 我们可以下载注释文件去获取基因的 TSS 位点坐标信息,但是每次自己去提还是会稍微麻烦一点。

于是写了个 python 脚本用来提取 Ensembl,Gencode 和 Ucsc 的基因的 TSS 位点信息, 输出结果为 bed 格式,后面可以直接可以给 deeptools 的 computeMatrix 和 plotProfile 计算绘图即可。

2基本思路

其实 linux 里面用 shell 也可做这样的事,也尝试一下用 python。

Ensembl 和 Gencode 数据库的 GTF 注释文件含有 gene 一行信息,包括染色体,起始位置,终止位置,正负链等,根据这些信息提取就行了, 正链则取坐标起始位置第一个碱基,负链则取右边结束最后一个碱基位置。此外 Ucsc 的注释文件是 没有 gene 这行信息的, 因此选用 transcript 的坐标作为参考。

需要注意的是,gtf 文件坐标系统是 1-based,而 bed 文件是 0-based,所以得 位置减去 1。

输出的 bed 文件应该是一个碱基的位置,我参考了一下 macs2 callpeak 的结果 summit 文件,应该是两个连续的值, summit 表示的是 peak 最高的位置。

3使用方法

使用方法很简单,python 运行脚本,提供三个参数即可。

例如提取 gencode 的 TSS 位点:

$ python GetTss.py gencode gencode.v19.annotation.gtf gencode.bed
You GTF file have: 57820 genes.

Your task has down!

bed 结果为:

$ head -5 gencode.bed
chr1    11868   11869   ENSG00000223972.4       .       +
chr1    29806   29807   ENSG00000227232.4       .       -
chr1    29553   29554   ENSG00000243485.2       .       +
chr1    36081   36082   ENSG00000237613.2       .       -
chr1    52472   52473   ENSG00000268020.2       .       +

提取 ucsc 的 TSS 位点:

$ python GetTss.py ucsc hg19.ncbiRefSeq.gtf myucsc.bed
You GTF file have: 104178 transcripts.

Your task has down!

bed 结果为:

$ head -5 myucsc.bed
chrMT   16023   16024   TRNP    .       -
chrMT   15887   15888   TRNT    .       +
chrMT   14746   14747   CYTB    .       +
chrMT   14742   14743   TRNE    .       -
chrMT   14673   14674   ND6     .       -

后面可以拿这个结果去计算富集密度了:

$ computeMatrix reference-point -S normal.bw treat.bw 
                -R myTSS.bed 
                --referencePoint center 
                -a 3000 -b 3000 -p 25 
                -out matrix.tab.gz

然后绘图:

$ plotProfile -m matrix.tab.gz 
              -out profile.pdf 
              --perGroup 
              --plotTitle 'test profile'

4python 脚本

以下是脚本代码: