基因的chr,start,end都是已知的 (这个文件需要下载,可以是UCSC,或者NCBI),也就是TSS文件。
测试文件见我的github.
import sys#引入sys模块
import os
os.chdir('D:/python/question10')
args = sys.argv#调用命令行参数
class Genome_info:#创建类Genome_info
def __init__(self):
self.chr = ""
self.start = 0
self.end = 0#初始化属性
class Gene(Genome_info):#创建父类Genome_info之下的类Gene
def __init__(self):
Genome_info.__init__(self)
self.orientation = ""
self.id = "" #初始化属性
list_chr = {} #定义染色体列表
with open('TSS.bed') as fp_gene: #导入参数1,即TSS.txt
for line in fp_gene:
if line.startswith("#"): #如果某行以#开头则越过
continue
lines = line.strip("\n").split("\t") #每行去除换行,以制表符分割
id = lines[0] #第一栏为基因id
chr = lines[1] #第二栏为染色体号
start = int(lines[2]) #第三栏转为整数型
end = int(lines[3]) #第四栏转为整数型
orientation = lines[4] #第五栏为基因方向
if not chr in list_chr: #如果染色体号在列表里不存在就初始化一下
list_chr[chr] = {}
gene = Gene() #初始化基因
gene.chr= chr #初始化染色体
gene.start = start #初始化基因起始位点
gene.end = end #初始化基因结束位点
gene.id = id #初始化ID
gene.orientation = orientation #初始化基因方向
list_chr[chr][id] = gene #将基因键、值存入list_chr字典
with open('pos.txt') as fp_pos: #导入参数2,即pos.txt
for line in fp_pos:
gene_list = [] #初始化gene_list
lines = line.strip('\n').split('\t') #每行去除换行,用制表符分割
(chr, start, end) = (lines[0], int(lines[1]), int(lines[2])) #取出染色体号,起始坐标和结束坐标,后两者均转为整数型
for gene_id, gene in list_chr[chr].items(): #判断pos.txt中基因位置与TSS.txt中是否有重叠
if gene.start <= start <= gene.end or gene.start <= end <= gene.end or start <= gene.start <= end or start <= gene.end <= end:
gene_list.append(gene.id) #如有则将基因ID添加至列表gene_list
print(gene_list) #输出gene_list