# PROGRAM USED TO CALCULATE PD DISTANCE
# Surendra Negi    
# March 10th 2014
print("\nHelp: python PD_Distance.py\n")
print("""Please cite following papers in your publications:
Ivanciuc, O., Schein, C.H. and Braun, W.  SDAP:  database and computational tools for allergenic proteins.  
Nucleic Acids Res.  31(1):359-362, 2003. 
Ivanciuc, O.I., Schein, C.H. and Braun, W.  Data mining of sequences and 3D structures of allergenic proteins.  
Bioinformatics  18(10):1358-1364, 2002.
Ivanciuc, O., Midoro-Horiuti, T., Schein, C.H., Xie, L., Hillman, G.R., Goldblum, R.M. and Braun, W.  
The property distance index PD predicts peptides that cross-react with IgE antibodies. 
Mol. Immunol., 46(5), 873-883, 2009.
Negi, S.S. and Braun, W.  Automated Detection of Conformational Epitopes using Phage Display Peptide Sequences.  
Bioinform. Biol. Insights,  3, 71-81, 2009.
Venkatarajan, M.S. and Braun, W.  New quantitative descriptors of amino acids based on 
multidimensional scaling of a large number of physical-chemical properties.  
J. Mol. Model.  7(12):445-453, 2001.\n\n\n""")

# Please feel free to modify your peptide sequences
# Input sequence1
Sequence1 = "VIPAARLFK"
# Input sequence2
Sequence2 = "VIPAARLWD"
# Gap penalty, Please feel free to modify this number. Default value is 20;
gap = 20.0

vector1 = {
    "A": 0.354, "R": 7.573, "N": 11.294, "D": 13.420, "C": -5.846,
    "Q": 6.599, "E": 9.788, "G": 9.655, "H": 1.019, "I": -15.634,
    "L": -11.825, "K": 10.762, "M": -10.585, "F": -14.571, "P": 7.662,
    "S": 8.813, "T": 3.012, "W": -13.110, "Y": -6.245, "V": -12.135, "-": 0.000
}

vector2 = {
    "A": 3.762, "R": -10.135, "N": 1.067, "D": -1.600, "C": 4.885,
    "Q": -5.166, "E": -7.861, "G": 15.778, "H": -4.969, "I": 1.993,
    "L": 0.505, "K": -9.517, "M": -3.959, "F": -0.646, "P": 8.029,
    "S": 6.682, "T": 4.127, "W": -5.222, "Y": -1.600, "V": 3.818, "-": 0.000
}

vector3 = {
    "A": -11.036, "R": 2.486, "N": 2.718, "D": -0.325, "C": 1.626,
    "Q": -0.697, "E": -7.318, "G": -0.558, "H": 0.953, "I": -2.045,
    "L": -6.157, "K": -1.022, "M": -3.601, "F": 1.673, "P": 9.456,
    "S": -0.348, "T": -0.348, "W": 9.038, "Y": 9.874, "V": -4.345, "-": 0.000
}

vector4 = {
    "A": -0.649, "R": -4.291, "N": 1.963, "D": 3.742, "C": 9.397,
    "Q": 0.582, "E": 2.611, "G": 0.299, "H": 4.657, "I": -3.243,
    "L": -4.557, "K": -5.405, "M": 5.339, "F": -0.033, "P": -3.576,
    "S": -1.131, "T": -2.195, "W": 1.380, "Y": -1.597, "V": -3.260, "-": 0.000
}

vector5 = {
    "A": 2.828, "R": -5.687, "N": -0.859, "D": 2.437, "C": -5.843,
    "Q": -1.750, "E": 4.734, "G": 1.656, "H": -0.328, "I": -1.672,
    "L": 3.219, "K": -0.422, "M": 1.203, "F": 3.250, "P": 6.000,
    "S": -3.062, "T": -4.281, "W": 4.640, "Y": -1.422, "V": -4.672, "-": 0.000
}

seq_lgth1 = len(Sequence1)
seq_lgth2 = len(Sequence2)

if seq_lgth1 != seq_lgth2:
    print(f"ERROR: Peptides are of different length. Length of Peptide1: {seq_lgth1} and Peptide2: {seq_lgth2}")
    exit()

PD_Distance = 0
for m in range(seq_lgth1):
    seq1 = Sequence1[m]
    seq2 = Sequence2[m]
    if seq1 == "X":
        seq1 = "-"
    if seq2 == "X":
        seq2 = "-"

    if seq1 == "-" or seq2 == "-":
        print("Gap or unknown residue found in peptide sequence")
        DifferencePD = (5.0 * gap * gap) ** 0.5
    else:
        value1 = vector1[seq1]
        value2 = vector2[seq1]
        value3 = vector3[seq1]
        value4 = vector4[seq1]
        value5 = vector5[seq1]
        xvalue1 = vector1[seq2]
        xvalue2 = vector2[seq2]
        xvalue3 = vector3[seq2]
        xvalue4 = vector4[seq2]
        xvalue5 = vector5[seq2]
        DifferencePD = ((value1 - xvalue1) ** 2 + (value2 - xvalue2) ** 2 + (value3 - xvalue3) ** 2 + 
                        (value4 - xvalue4) ** 2 + (value5 - xvalue5) ** 2) ** 0.5

    PD_Distance += DifferencePD

finalpd_distance = PD_Distance / seq_lgth1
print(f"Sequence 1: {Sequence1} :: Sequence 2: {Sequence2} :: PD_DISTANCE: {finalpd_distance}")