rrwick / compare-annotations Goto Github PK

 match, identity, length_diff = compare_features(old_feature, new_feature, old_record, new_record, args.match_identity_threshold)
            if match:
                for j in range(old_offset):
                    print_in_old_not_new(old_features[old_i + j])
                for j in range(new_offset):
                    print_in_new_not_old(new_features[new_i + j])
                print_match(old_features[old_i + old_offset], new_features[new_i + new_offset], identity, length_diff)
                old_i += old_offset
                new_i += new_offset
                break
        else:
            sys.exit('\nERROR: Failed to find alignment')

        if old_feature is None and new_feature is None:
            break

        old_i += 1
        new_i += 1


def print_match(f1, f2, identity, length_diff):
    print('', flush=True)
    if identity == 1.0:
        print('Exact match')
    else:
        print('Inexact match (' + '%.2f' % (identity * 100.0) + '% ID, ', end='')
        if length_diff == 0:
            print('same length)')
        elif length_diff > 0:
            print('old seq longer)')
        elif length_diff < 0:
            print('new seq longer)')
    print('  old: ', end='')
    print_feature_one_line(f1)
    print('  new: ', end='')
    print_feature_one_line(f2)
    p1 = f1.qualifiers['product'][0].lower()
    p2 = f2.qualifiers['product'][0].lower()
    if 'hypothetical' in p1 and 'hypothetical' in p2:
        print('  still hypothetical')
    if 'hypothetical' in p1 and 'hypothetical' not in p2:
        print('  no longer hypothetical')
    if 'hypothetical' not in p1 and 'hypothetical' in p2:
        print('  became hypothetical')


def print_in_old_not_new(f):
    print('')
    print('In old but not in new:')
    print('  ', end='')
    print_feature_one_line(f)


def print_in_new_not_old(f):
    print('')
    print('In new but not in old:')
    print('  ', end='')
    print_feature_one_line(f)


def print_feature_one_line(f):
    f_str = f.qualifiers['product'][0]
    strand = '+' if f.location.strand == 1 else '-'
    f_str += ' (' + str(f.location.start) + '-' + str(f.location.end) + ' ' + strand + ', '
    f_str += str(f.location.end - f.location.start) + ' bp)'
    print(f_str)


def compare_features(f1, f2, r1, r2, match_identity_threshold):
    if f1 is None or f2 is None:
        return False
    s1 = f1.extract(r1).seq
    s2 = f2.extract(r2).seq
    score = pairwise2.align.globalms(s1, s2, 1, 0, 0, 0, score_only=True)
    identity = score / max(len(s1), len(s2))
    match = identity > match_identity_threshold
    length_diff = len(s1) - len(s2)
    return match, identity, length_diff


if __name__ == '__main__':
    main()