regtools/upstreamGetter at master · combogenomics/regtools · GitHub

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
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
#!/usr/bin/python

############################
#            Imports             #
############################

from Bio import SeqIO
import argparse
from regtools import *

def get_replicons(strain):
    # curl is required... if not present, please
    # sudo apt-get install curl
    os.system('bash get_replicons.sh %s > %s' %(strain,'tmp'))
    replicons=[line.strip()[:-4] for line in open('tmp')]
    os.remove('tmp')
    return replicons

def export_upstreams(upstreams,nome_out,infos):
    out=open(nome_out,'w')
    organism,molecule=infos
    for up in upstreams:
        name,from_,size,start,end,strand,b_up,b_down=up.infos
        # upstream bases
        upstream = -(size-b_down)
        if upstream > 0:
            upstream=0
        to_write='>%s upstream from %s to %s; size: %s; feature type: cds; location:%s:%s:%s:%s:%s\n' %(name,upstream,b_down,size,organism,molecule,start,end,strand)
        seq=nice_write(up.seq)
        out.write(to_write)
        out.write(seq)

############################
#             Inputs             #
############################

def getOptions():
    description = "Get up and downstream sequences from genomes"
    parser = argparse.ArgumentParser(description = description)
    parser.add_argument('orgID', action='store',
                            help='Organism ID, as in the NCBI database (or genbank file if -f)')
    parser.add_argument('-u', metavar='upstream', action='store',
                        dest='b_up',
                        type=int,
                        default=400,
                        help='Upstream bases')
    parser.add_argument('-d', metavar='downstream', action='store',
                        dest='b_down',
                        type=int,
                        default=100,
                        help='Downstream bases')
    parser.add_argument('-f', '--file', action="store_true",
                            default=False,
                        help='Use a genbank file instead of downloading from NCBI (useful for draft unannotated genomes)')
    parser.add_argument('-l', '--linear', action="store_true",
                            default=False,
                        help='Linear or draft genome')
    parser.add_argument('-o', '--one-file', action="store_true",
                            dest='one_file',
                            default=False,
                        help='Save everything in one file')
    parser.add_argument('-a', '--all', action="store_true",
                            dest='all',
                            default=False,
                        help='Save ALL upstream regions (not considering genes overlaps)')

    return parser.parse_args()

if __name__=='__main__':
    options = getOptions()
    organism = options.orgID
    is_file = options.file
    b_up = options.b_up
    b_down = options.b_down
    is_circ = not options.linear
    one_file = options.one_file
    all_up = options.all

############################
#              Main              #
############################

if __name__ == '__main__':
    # do stuff
    all_upstream = []
    if is_file:
        sequences = [s for s in SeqIO.parse(open(organism), 'genbank')]
        replicons = [s.id for s in sequences]
    else:
        replicons = get_replicons(organism)

    print 'working on these molecules...'
    for replicon in replicons:
        print replicon
        if is_file:
            record = filter(lambda x: x.id == replicon, sequences)[0]
        else:
            record = get_from_genbank(replicon)
        genes=get_genes(record)
        if len(genes)==0:
            print 'no genes in', replicon
            continue
        whole_seq=get_whole_genome_seq(record)
        if not all_up:
            upstreams=get_all_upstreams(genes,whole_seq, b_up, b_down, is_circ)
        else:
            upstreams=get_all_upstreams_all(genes,whole_seq, b_up, b_down, is_circ)
        if not upstreams:continue
        file_name='%s_%s_upstreams.fasta' %(organism,replicon)
        infos=(organism,replicon)
        if not one_file:
            export_upstreams(upstreams,file_name,infos)
        else:
            for up in upstreams:
                all_upstream.append(up)
    if one_file:
        export_upstreams(all_upstream, '%s_upstreams.fasta' %(organism), (organism,'all'))