最近有个任务需要从基因组中提取所有的内含子（intron）区域。一开始觉得很简单，用ensembl上的biomart或者ucsc genome browser一顿猛点应该就可以搞定。结果上去一顿操作下来发现不行啊，只能提取功能feature （CDS, exon, UTR，TSS之类的），没有intron啊。网上有大牛说是用biomart 或者ucsc可以，但是不想再花更多时间去摸了。之前写了一个python程序实现这个功能了，但是最近学校服务器坏了，修了半个月了还没好，也不能瞎等，于是再写了一遍，还是放到这里来比较方便，一来给有需要的朋友，二来，下次服务器坏了，也不至于要重写。
代码是基于python2.7的，需要的同学可以拷贝过去直接使用。
#!~/tools/anaconda2/bin/python
'''
Extract intron from given gtf file, and generate bed files
Usage: python extract_intron_from_gtf.py --gtf <xxx.gtf> --out <xxx.bed>
'''
import argparse
import re
# Function --------------------------------
def extract_exon(gtf_file):
keep_chr = []#需要保留的染色体
for i in range(22):
keep_chr.append('chr'+str(i+1))
keep_chr = keep_chr + ['chrX','chrY']
with open(gtf_file,'r') as fh:
transcript_exon = {}
geneid,transcript_id = '',''
for line in fh.readlines():
line = line.strip()
w = line.split("\t")
chrn = w[0]
strand = w[6]
entry_type = w[2]
geneid_patt = re.search(r'gene_id\s"(\w+)";',w[8])
transcript_patt = re.search(r'transcript_id\s"(\w+)";',w[8])
if geneid_patt:
gene_id = geneid_patt.group(1)
else:
gene_id = "NA"
if transcript_patt:
transcript_id = transcript_patt.group(1)
transcript_id = chrn+":"+transcript_id+":"+gene_id+":"+strand
if entry_type == "exon" and chrn in keep_chr:
exon = (int(w[3]),int(w[4]))
if transcript_id not in transcript_exon:
transcript_exon[transcript_id] = [exon]
else:
transcript_exon[transcript_id].append(exon)
else:
continue
else:
continue
return transcript_exon
def extract_intron(exon_list):
intron_list = []
for i in range(len(exon_list)-1):
intron = (exon_list[i][1],exon_list[i+1][0])
intron_list.append(intron)
return intron_list
# Main -----------------------------------
parser = argparse.ArgumentParser(description='')
parser.add_argument('--gtf',help="specify a gtf file")
parser.add_argument('--out',help="specify a filename for output")
args = parser.parse_args()
exon_dict = extract_exon(args.gtf) #调用上面定义的函数读取gtf
fho = open(args.out,'w')
for t in exon_dict:
chrname,tid,gid,strand = t.split(":")
if len(exon_dict[t]) > 1:
introns = extract_intron(exon_dict[t])
for intron in introns:
print >>fho,"%s\t%d\t%d\t%s\t%d\t%s" % (chrname,intron[0]-1,intron[1],gid+":"+tid,intron[1]-intron[0]+1,strand)
fho.close()
print "Generate intron bed in file:%s" % (args.out)
print "Done!"
gtf文件里面可能会存在alternative contig, 我这里会自动去除。只保留了染色体1-22，X和Y上面的内含子。另外，有些gtf文件里面，ensembl geneid 是带有版本号的（比如ENSG00000469289.1）,这种情况，代码里面需要稍作调整（re 模块正则表达式里面需要添加版本匹配）。嗯，暂时就这样吧，有其他问题的朋友可以留言。
更新2021-05-12：我放了一个更新的版本在github上，功能更加完善，能同时产生exon和intron两个bed文件。
仓库链接：https://github.com/Xu-Dong/Exon_Intron_Extractor
