code: def annotate_orfs(fa, seqs, threads, threshold = 1e-6):
      
      db = '%s/pfam/Pfam-A.hmm' % os.environ['databases']
      base = fa.rsplit('.', 1)[0]
      out = '%s.pfam.hmmscan' % (base)
      tblout = '%s.pfam.hmmscan-tblout' % (base)
      orf2pfam = {} # orf2pfam[scaffold id][orf id] = pfam
      if os.path.exists(tblout) is False:
          p = subprocess.Popen('hmmscan --cpu %s --tblout %s -o %s %s %s'\
                  % (threads, tblout, out, db, fa), shell = True)
          p.communicate()
      for line in open(tblout):
          if line.startswith('#'):
              continue
          line = line.strip().split()
          hit, id, e, inc = line[0], line[2], float(line[4]), int(line[17])
          contig = id.split('_', 1)[0]
          if e > threshold:
              continue
          if contig not in orf2pfam:
              orf2pfam[contig] = {}
          if id not in orf2pfam[contig]:
              orf2pfam[contig][id] = hit
      for seq, info in list(seqs.items()):
          for insertion in info[2]:
              annotations = {}
              for orf in insertion[4]:
                  orf  = orf[0].split('>')[1].split()[0]
                  contig = orf.split('_', 1)[0]
                  if contig in orf2pfam and orf in orf2pfam[contig]:
                      annotations[orf] = orf2pfam[contig][orf]
                  else:
                      annotations[orf] = 'n/a'
              insertion.append(annotations)
      return seqs
docs: search orfs against pfam
      # seqs[id] = [gene, model, [[i-gene_pos, i-model_pos, i-length, iseq, [orfs], [introns], orfs?, introns?, {orf: annotations}], ...]]