summaryrefslogtreecommitdiff
path: root/Code/prepare_gene_file.py
blob: febcbef6b9601987b73c326a91e0470dc3ac4f2e (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
# Usage: python prepare_gene_file.py all-ath-gene-position.txt > gene_file.txt
#        all-ath-gene-position contains all gene IDs in the genome, and their positions, orientation, etc, in BED format.
#        See ../Data/information/all-ath-gene-position.txt
# Purpose: get gene name and gene annotation
# 2 JAN 2017 hui SLCU

import sys
import os

###################################################################################
GENE_DESCRIPTION  =  '../Data/information/gene_description_20140101.txt'
AGI_TO_GENE_NAMES =  '../Data/information/AGI-to-gene-names.txt'
###################################################################################


def get_description(x, d):
        result = ''
	if x in d:
	        result = '\t' + d[x]
	else:
		result = '\tNot Found'
        return result
        

def make_AGI_to_gene_name_dict(fname):
    d = {}
    f = open(fname)
    for line in f:
        line = line.strip()
        lst = line.split()
        agi = lst[0]
        name = lst[1]
        if not agi in d:
            d[agi] = name
        else:
            d[agi] += ';' + name
    f.close()
    return d


## main

# make a dictionary of gene description
f0 = open(GENE_DESCRIPTION)
d = {}
for line in f0:
	line = line.strip()
	lst = line.split('\t')
	id = lst[0]
        id = id[0:9] # AGI id, omit .1, .2, .3, etc
        s = '\t'.join(lst[1:])
	if not id in d:
		d[id] = s
        else:
                d[id] += '\t' + s

f0.close()

agi2genename_dict = make_AGI_to_gene_name_dict(AGI_TO_GENE_NAMES)

locus_file = sys.argv[1] # location of genes
f = open(locus_file) # see all-ath-gene-position.txt
for line in f:
        line = line.strip()
        lst = line.split('\t')
        x = lst[3]
        c = lst[0]
        ss  = lst[1]
        ee  = lst[2]
        strand = lst[5]
        result = [x]
        if x in agi2genename_dict and not x == agi2genename_dict[x]:
                result.append(agi2genename_dict[x])
        else:
                result.append(' ') # if no gene name, use a space
        result.extend([c, ss, ee, strand, get_description(x, d)])
        print('\t'.join(result))

f.close()