用python根据指定染色体及坐标得到位置信息(第十题)

2020-03-05  本文已影响0人  多啦A梦詹

基因的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
上一篇 下一篇

猜你喜欢

热点阅读