search orfs against pfam
    # seqs[id] = [gene, model, [[i-gene_pos, i-model_pos, i-length, iseq, [orfs], [introns], orfs?, introns?, {orf: annotations}], ...]]

def annotate_orfs(fa, seqs, threads, threshold = 1e-6):
    """
    search orfs against pfam
    # seqs[id] = [gene, model, [[i-gene_pos, i-model_pos, i-length, iseq, [orfs], [introns], orfs?, introns?, {orf: annotations}], ...]]
    """
    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