repseq_vdj.py 3.93 KB
Newer Older
1
'''
2
Parses output of various RepSeq programs.
3 4 5
Takes either:
  - a .fa file, a _Summary.txt file as produced by IMGT/V-QUEST
  - or a results file produced by MiXCR
6 7
and creates a .vdj file to be checked by should-vdj-to-tap.py

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

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():
25
        if term in ['Homsap', '[F]', '(F)', 'F', 'P', 'or', 'and', '(see', 'comment)', 'ORF', '[ORF]']:
26 27 28 29 30 31 32 33 34 35 36 37
            continue
        genes += [term]

    if not genes:
        return ''

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

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


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

41 42
def CDR3_to_vdj(s):
    return '{%s}' % s if s else ''
43

44 45
class Result():
    '''Stores a tabulated result'''
46 47 48 49

    def __init__(self, l):

        self.d = {}
50
        self.result = self.parse(l)
51 52

        for i, data in enumerate(l.split('\t')):
53
            self.d[self.labels[i]] = data
54 55 56 57

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

58 59 60 61
    def __str__(self):
        return str(self.d)


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
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']

89 90 91
        s += ' '
        s += CDR3_to_vdj(self['AA. Seq. CDR3'])

92 93
        return s

94 95 96 97 98 99 100 101

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)

102 103 104 105 106 107 108 109 110 111 112
    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'])
113 114 115 116

        s += ' '
        s += CDR3_to_vdj(self['AA JUNCTION'])

117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141
        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()

Mathieu Giraud's avatar
Mathieu Giraud committed
142
        r = IMGT_VQUEST_Result(vquest)
143
        yield (fasta.replace('>', ''), r.to_vdj())
144 145


146 147 148 149 150 151 152 153 154 155 156 157 158
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()


159 160
if __name__ == '__main__':

161
    if 'mixcr' in sys.argv[1]:
162 163 164 165 166 167 168 169
        gen = header_mixcr_results(sys.argv[1])
    else:
        gen = header_vquest_results(sys.argv[1], sys.argv[2])

    # output .vdj data
    for (header, result) in gen:
        print "#%s" % header
        print ">%s" % result
170
        print