guohanghui/SPM

1#!/usr/bin/env python32# -*- coding: utf-8 -*-3 4# This python script for Sequence Pattern Matching (SPM) alignment was written by WuLab & YanLab, School of Life Science, Westlake University.5#6# It can help you search target protein sequence with your input peptides sequence from your given sequence database. 7# Input variants are a peptide sequence and a protein sequence database. And output is a ranking protein list, each with a score and matching postion. 8# The most possible searching candidate is at the top of the output files with the lowest score. 9#10# Variant:11# query_seq: String type. It is your given one-letter sequence, which can be as short as 10-20 residues, the longer, the better. The script will search proteins based on it.12# db_fasta: String type. The protein sequence database file's path and name.13# output_dir: String type. It is the output file's path and name.14 15# All rights reserved. Please cite us if you use it.16 17import numpy as np18import os19import argparse20 21 22volume = {'A':15, 'C':47, 'D':59, 'E':73, 'F':91, 'G':1, 'H':81, 'I':57, 'K':72, 'L':57, 23          'M':75, 'N':58, 'P':41, 'Q':72, 'R':100, 'S':31, 'T':45, 'V':43, 'W':130, 'Y':107, 'X':0}24 25def volumeScoring(query_seq_volume, uniprot_info, db_seq):26    """27    Arguments:28        - query_seq_volume: list of residue volumes of query sequence.29        - uniprot_info: information read from .fasta database.30        - db_seq: corresponding uniprot sequence.31    """32    query_len = len(query_seq_volume)33    db_seq_volume = np.array([volume[i] for i in db_seq])34    score = 999935    position = 136    for i in range(len(db_seq_volume)-query_len):37            s = np.sum(np.abs(db_seq_volume[i:i+query_len]-query_seq_volume))38            score = min(s, score)39            if s <= score:40                position = i+141    return [uniprot_info, score, position-1]42 43def loadUniprotDB(db_fasta):44    """45    Argument:46        - db_fasta: directory of .fasta file containing a set of sequences against which you want to query.47    """48    uniprot_info, seq, database = None, None, dict()49    for i in open(db_fasta).readlines():50        if i[0] == '>':51            database[uniprot_info] = seq52            uniprot_info, seq = i.replace(' ', '_').strip(), ''53        else:54            seq += i.split()[0]55    database[uniprot_info] = seq56    database.pop(None)57    return database58 59def peptideSearching(db_fasta, query_seq, output_file):60    """61    Arguments:62        - db_fasta: directory of .fasta file containing a set of sequences against which you want to query.63        - query_seq: string of protein sequence in one-letter code without gapping, e.g. 'DKLSPIRRAAVVN'.64        - output_file: file of output containing information of scoring and best matching positions, e.g. '/ssd/output.txt'65    """66    database = loadUniprotDB(db_fasta)67    query_seq_volume = np.array([volume[i] for i in query_seq])68    score = [volumeScoring(query_seq_volume, i, j) for i,j in database.items()]69    np.savetxt(output_file, np.array(score), fmt='%s', header='#UniProt_INFO #Score #Position')70    os.system('bash scripts/SPM_ranking.sh %s' % (output_file))71 72 73if __name__ == "__main__":74    parser = argparse.ArgumentParser("Peptide Searching by WuLab & YanLab in Westlake University")75    parser.add_argument('-q', '--query_seq', type=str, help='Input query sequence in one-letter code without gapping, e.g. DKLSPIRRAAVVN')76    parser.add_argument('-d', '--db_fasta', type=str, help='Input database fasta file, e.g. uniprot.fasta')77    parser.add_argument('-o', '--output_file', type=str, help='Output file path, e.g. output.txt')78    args = parser.parse_args()79    80    81    #query_seq = 'RLMHARFIAWKII'82    query_seq = args.query_seq83    db_fasta = args.db_fasta84    output_file = args.output_file85    86    peptideSearching(db_fasta=db_fasta, query_seq=query_seq, output_file=output_file)87    88