|
| 1 | +import json |
| 2 | +import logging |
| 3 | +from argparse import Namespace |
| 4 | +from typing import Any |
| 5 | + |
| 6 | +from Bio import SeqIO |
| 7 | +from Bio.SeqUtils.ProtParam import ProteinAnalysis |
| 8 | + |
| 9 | +from jsrc.core import DataFormatError |
| 10 | + |
| 11 | +logger = logging.getLogger(__name__) |
| 12 | + |
| 13 | +_AA_VOLUMES: dict[str, float] = { |
| 14 | + "A": 1.0, "R": 6.13, "N": 2.95, "D": 2.78, "C": 2.43, |
| 15 | + "Q": 3.95, "E": 3.78, "G": 0.0, "H": 4.66, "I": 4.0, |
| 16 | + "L": 4.0, "K": 4.77, "M": 4.43, "F": 5.89, "P": 2.72, |
| 17 | + "S": 1.6, "T": 2.6, "W": 8.08, "Y": 6.47, "V": 3.0, |
| 18 | +} |
| 19 | + |
| 20 | + |
| 21 | +def _aliphatic_index(seq: str) -> float: |
| 22 | + """Aliphatic index per Ikai (1980).""" |
| 23 | + aa = seq.upper() |
| 24 | + total = len(aa) |
| 25 | + if total == 0: |
| 26 | + return 0.0 |
| 27 | + a = aa.count("A") |
| 28 | + v = aa.count("V") |
| 29 | + ile = aa.count("I") |
| 30 | + leu = aa.count("L") |
| 31 | + return (a + 2.9 * v + 3.9 * (ile + leu)) / total * 100.0 |
| 32 | + |
| 33 | + |
| 34 | +def cmd(args: Namespace) -> None: |
| 35 | + records = list(SeqIO.parse(args.fa, "fasta")) |
| 36 | + if not records: |
| 37 | + raise DataFormatError("No protein sequences found in FASTA") |
| 38 | + |
| 39 | + results = [] |
| 40 | + for rec in records: |
| 41 | + seq = str(rec.seq).upper().replace("U", "T") |
| 42 | + clean = "".join(c for c in seq if c.isalpha()) |
| 43 | + if not clean: |
| 44 | + continue |
| 45 | + pa = ProteinAnalysis(clean) |
| 46 | + try: |
| 47 | + mw = pa.molecular_weight() |
| 48 | + except Exception: |
| 49 | + mw = 0.0 |
| 50 | + try: |
| 51 | + pi = pa.isoelectric_point() |
| 52 | + except Exception: |
| 53 | + pi = 0.0 |
| 54 | + try: |
| 55 | + ec = pa.molar_extinction_coefficient() |
| 56 | + except Exception: |
| 57 | + ec = (0, 0) |
| 58 | + try: |
| 59 | + ii = pa.instability_index() |
| 60 | + except Exception: |
| 61 | + ii = 0.0 |
| 62 | + try: |
| 63 | + ai = _aliphatic_index(clean) |
| 64 | + except Exception: |
| 65 | + ai = 0.0 |
| 66 | + try: |
| 67 | + gv = pa.gravy() |
| 68 | + except Exception: |
| 69 | + gv = 0.0 |
| 70 | + try: |
| 71 | + ar = pa.aromaticity() |
| 72 | + except Exception: |
| 73 | + ar = 0.0 |
| 74 | + try: |
| 75 | + ss = pa.secondary_structure_fraction() |
| 76 | + except Exception: |
| 77 | + ss = (0.0, 0.0, 0.0) |
| 78 | + charge = pa.charge_at_pH(args.ph) if args.ph is not None else None |
| 79 | + |
| 80 | + entry: dict[str, Any] = { |
| 81 | + "id": rec.id, |
| 82 | + "length": len(clean), |
| 83 | + "molecular_weight": round(mw, 2), |
| 84 | + "isoelectric_point": round(pi, 2), |
| 85 | + "extinction_coefficient_reduced": round(ec[0], 2), |
| 86 | + "extinction_coefficient_oxidized": round(ec[1], 2), |
| 87 | + "instability_index": round(ii, 2), |
| 88 | + "aliphatic_index": round(ai, 2), |
| 89 | + "gravy": round(gv, 4), |
| 90 | + "aromaticity": round(ar, 4), |
| 91 | + "helix_fraction": round(ss[0], 4), |
| 92 | + "turn_fraction": round(ss[1], 4), |
| 93 | + "sheet_fraction": round(ss[2], 4), |
| 94 | + } |
| 95 | + if charge is not None: |
| 96 | + entry["charge_at_pH"] = round(charge, 4) |
| 97 | + results.append(entry) |
| 98 | + |
| 99 | + if args.json: |
| 100 | + print(json.dumps(results, ensure_ascii=False, indent=2)) |
| 101 | + return |
| 102 | + header = ( |
| 103 | + "id\tlength\tmw\tpI\tEC_red\tEC_ox\tinstability\taliphatic\t" |
| 104 | + "gravy\taromaticity\thelix\tturn\tsheet" |
| 105 | + ) |
| 106 | + if args.ph is not None: |
| 107 | + header += "\tcharge" |
| 108 | + print(header) |
| 109 | + for r in results: |
| 110 | + row = ( |
| 111 | + f"{r['id']}\t{r['length']}\t{r['molecular_weight']}\t" |
| 112 | + f"{r['isoelectric_point']}\t{r['extinction_coefficient_reduced']}\t" |
| 113 | + f"{r['extinction_coefficient_oxidized']}\t{r['instability_index']}\t" |
| 114 | + f"{r['aliphatic_index']}\t{r['gravy']}\t{r['aromaticity']}\t" |
| 115 | + f"{r['helix_fraction']}\t{r['turn_fraction']}\t{r['sheet_fraction']}" |
| 116 | + ) |
| 117 | + if args.ph is not None: |
| 118 | + row += f"\t{r.get('charge_at_pH', '')}" |
| 119 | + print(row) |
| 120 | + logger.info("Analyzed %d protein sequences", len(results)) |
| 121 | + |
| 122 | + |
| 123 | +def register(subparsers: Any) -> None: |
| 124 | + p = subparsers.add_parser( |
| 125 | + "protparam", help="Protein physicochemical properties" |
| 126 | + ) |
| 127 | + p.add_argument("-fa", required=True, help="Protein FASTA file") |
| 128 | + p.add_argument("--json", action="store_true", help="Print JSON output") |
| 129 | + p.add_argument("--ph", type=float, default=None, help="pH for net charge") |
| 130 | + p.set_defaults(func=cmd) |
0 commit comments