30 lines
2.5 KiB
Python
30 lines
2.5 KiB
Python
#!/usr/bin/env python3
|
|
from __future__ import annotations
|
|
import argparse
|
|
from pathlib import Path
|
|
import pandas as pd
|
|
import yaml
|
|
from panfamily.io import read_fasta
|
|
|
|
|
|
def main():
|
|
ap=argparse.ArgumentParser(); ap.add_argument("--candidates",required=True,type=Path); ap.add_argument("--hmm-hits",required=True,type=Path); ap.add_argument("--homology-hits",required=True,type=Path); ap.add_argument("--proteins",required=True,type=Path); ap.add_argument("--rules",required=True,type=Path); ap.add_argument("--output",required=True,type=Path); args=ap.parse_args()
|
|
candidates=pd.read_csv(args.candidates,sep="\t",dtype={"gene_id":str}); hmm=pd.read_csv(args.hmm_hits,sep="\t",dtype={"sequence_id":str}); hom=pd.read_csv(args.homology_hits,sep="\t",names=["query_id","subject_id","identity","alignment_length","mismatch","gapopen","qstart","qend","sstart","send","evalue","bitscore","query_length","subject_length"],dtype={"query_id":str,"subject_id":str}); proteins=read_fasta(args.proteins); rules=yaml.safe_load(args.rules.read_text(encoding="utf-8"))
|
|
min_cov=float(rules.get("member_filter",{}).get("minimum_core_hmm_coverage",0.5)); frag=float(rules.get("member_filter",{}).get("fragment_hmm_coverage",0.4))
|
|
best_hmm=hmm.sort_values(["sequence_id","domain_ievalue","domain_score"],ascending=[True,True,False]).drop_duplicates("sequence_id").set_index("sequence_id") if not hmm.empty else pd.DataFrame()
|
|
best_hom=hom.sort_values(["subject_id","evalue","bitscore"],ascending=[True,True,False]).drop_duplicates("subject_id").set_index("subject_id") if not hom.empty else pd.DataFrame()
|
|
rows=[]
|
|
for row in candidates.itertuples(index=False):
|
|
gene=row.gene_id; hc=float(best_hmm.loc[gene,"hmm_coverage"]) if gene in best_hmm.index else 0.0; has_h=bool(row.has_hmm); has_d=bool(row.has_reference_homology)
|
|
if has_h and has_d and hc>=min_cov: level="A"
|
|
elif has_h and hc>=min_cov: level="B"
|
|
elif has_d and hc<min_cov: level="C"
|
|
elif has_h and hc>=frag: level="Fragment"
|
|
else: level="Rejected"
|
|
rows.append({"gene_id":gene,"accession":row.accession,"member_level":level,"has_core_domain":int(has_h),"has_reference_homology":int(has_d),"hmm_coverage":hc,"protein_length":len(proteins.get(gene,"")),"domain_status":"complete" if hc>=min_cov else "partial" if hc>0 else "missing"})
|
|
args.output.parent.mkdir(parents=True,exist_ok=True); pd.DataFrame(rows).to_csv(args.output,sep="\t",index=False)
|
|
|
|
|
|
if __name__ == "__main__": main()
|
|
|