imgt-to-vdj.py 3.82 KB
Newer Older
1
'''
2 3 4
Takes either:
  - a .fa file, a _Summary.txt file as produced by IMGT/V-QUEST
  - or a results file produced by MiXCR
5 6
and creates a .vdj file to be checked by should-vdj-to-tap.py

7 8 9
python imgt-to-vdj.py data-curated/curated_IG.fa data-curated/curated_ig_Summary.txt > data-curated/imgt-IG.vdj
python imgt-to-vdj.py data-curated/curated_TR.fa data-curated/curated_tr_Summary.txt > data-curated/imgt-TR.vdj
python imgt-to-vdj.py data-curated/mixcr.results > data-curated/mixcr.vdj
10 11 12 13 14 15 16 17 18 19 20 21 22 23
'''

import sys


def parse_gene_and_allele_to_vdj(s):
    '''
    Parse lines such as:
    Homsap IGHV3-30*03 F, or Homsap IGHV3-30*18 F or Homsap IGHV3-30-5*01 F'
    and produce a .vdj line
    '''

    genes = []
    for term in s.replace(',', '').split():
24
        if term in ['Homsap', '[F]', '(F)', 'F', 'P', 'or', 'and', '(see', 'comment)', 'ORF', '[ORF]']:
25 26 27 28 29 30 31 32 33 34 35 36
            continue
        genes += [term]

    if not genes:
        return ''

    if len(genes) == 1:
        return genes[0]

    return '(%s)' % ','.join(genes)


37 38 39 40
def N_to_vdj(s):
    return '/%s/' % s


41 42
class Result():
    '''Stores a tabulated result'''
43 44 45 46

    def __init__(self, l):

        self.d = {}
47
        self.result = self.parse(l)
48 49

        for i, data in enumerate(l.split('\t')):
50
            self.d[self.labels[i]] = data
51 52 53 54

    def __getitem__(self, key):
        return self.d[key]

55 56 57 58
    def __str__(self):
        return str(self.d)


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
class MiXCR_Result(Result):

    def parse(self, l):
        self.labels = mixcr_labels
        return ('\t' in l.strip())

    def to_vdj(self):

        if not self.result:
            return 'no result'

        s = ''
        s += self['Best V hit']
        s += ' '

        if self['Best D hit']:
            s += N_to_vdj(self['N. Seq. VDJunction'])
            s += ' '
            s += self['Best D hit']
            s += ' '
            s += N_to_vdj(self['N. Seq. DJJunction'])
        else:
            s += N_to_vdj(self['N. Seq. VJJunction'])

        s += ' '
        s += self['Best J hit']

        return s

88 89 90 91 92 93 94 95

class IMGT_VQUEST_Result(Result):
    '''Stores a IMGT/V-QUEST result'''

    def parse(self, l):
        self.labels = vquest_labels
        return ('No result' not in l)

96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134
    def to_vdj(self):

        if not self.result:
            return 'no result'

        s = ''
        s += parse_gene_and_allele_to_vdj(self['V-GENE and allele'])
        s += ' '
        s += parse_gene_and_allele_to_vdj(self['D-GENE and allele'])
        s += ' '
        s += parse_gene_and_allele_to_vdj(self['J-GENE and allele'])
        return s



def header_vquest_results(ff_fasta, ff_vquest):
    f_fasta = open(ff_fasta).__iter__()
    f_vquest = open(ff_vquest).__iter__()

    vquest_first_line = f_vquest.next()
    globals()['vquest_labels'] = vquest_first_line.split('\t')

    # print vquest_labels

    while True:

        fasta = ''
        # Advance until header line
        while not '>' in fasta:
            fasta = f_fasta.next().strip()

        vquest = ''
        # Advance until non-empty line
        while not vquest:
            vquest = f_vquest.next().strip()

        yield (fasta, vquest)


135 136 137 138 139 140 141 142 143 144 145 146 147
def header_mixcr_results(ff_mixcr):

    f = open(ff_mixcr).__iter__()

    mixcr_first_line = f.next()
    globals()['mixcr_labels'] = mixcr_first_line.split('\t')

    while True:
        l = f.next()
        result = MiXCR_Result(l)
        yield result['Description R1'], result.to_vdj()


148 149
if __name__ == '__main__':

150 151 152 153 154 155 156 157
    if 'mixcr' in sys.argv[1]:
        for (header, result) in header_mixcr_results(sys.argv[1]):
            print "#%s" % header
            print ">%s" % result
            print

        sys.exit(0)

158 159 160 161 162 163 164
    for (header, result) in header_vquest_results(sys.argv[1], sys.argv[2]):
        # print "=========="
        print header.replace('>', '#')
        r = IMGT_VQUEST_Result(result)
        # print r
        print ">%s" % r.to_vdj()
        print