def search_and_extract_orfs_matching_protein_database(self,
                                                      unpack,
                                                      search_method,
                                                      maximum_range,
                                                      threads,
                                                      evalue,
                                                      min_orf_length,
                                                      restrict_read_length,
                                                      diamond_database,
                                                      output_search_file,
                                                      hit_reads_fasta,
                                                      hit_reads_orfs_fasta):
        '''As per aa_db_search() except slightly lower level. Search an
        input read set (unpack) and then extract the proteins that hit together
        with their containing nucleotide sequences.

        Parameters
        ----------
        output_search_file: str
            path to hmmsearch output table or diamond basename
        hit_reads_fasta: str
            path to nucleotide sequences containing hit proteins
        hit_reads_orfs_fasta: str
            path to hit proteins, unaligned
            
        Returns
        -------
        direction_information: dict
                    
            {read_1: False
             ...
             read n: True}
            
            where True = Forward direction
            and False = Reverse direction
        
        result: DBSearchResult object containing file locations and hit 
        information
        
        '''

        # Define method of opening sequence files to stdout
        if unpack.is_zcattable():
            clazz = ZcatOrfM
        else:
            clazz = OrfM
        orfm = clazz(min_orf_length=min_orf_length,
                     restrict_read_length=restrict_read_length)
        extracting_orfm = OrfM(min_orf_length=min_orf_length,
                      restrict_read_length=restrict_read_length)

        if search_method == 'hmmsearch':
            # run hmmsearch
            search_result = self.hmmsearch(
                                           output_search_file,
                                           unpack.read_file,
                                           unpack,
                                           unpack.sequence_type(),
                                           threads,
                                           evalue,
                                           orfm
                                           )

        elif search_method == 'diamond':
            # run diamond
            search_result = Diamond(
                                     database=diamond_database,
                                     threads=threads,
                                     evalue=evalue,
                                     ).run(
                                           unpack.read_file,
                                           unpack.sequence_type(),
                                           daa_file_basename=output_search_file
                                           )
            search_result = [search_result]

        else:  # if the search_method isn't recognised
            raise Exception("Programming error: unexpected search_method %s" % search_method)

        orfm_regex = OrfM.regular_expression()
        
        hits = self._get_sequence_directions(search_result)
    
        orf_hit_readnames = hits.keys() # Orf read hit names
        if unpack.sequence_type() == 'nucleotide':
            hits={(orfm_regex.match(key).groups(0)[0] if orfm_regex.match(key) else key): item for key, item in hits.iteritems()}
            hit_readnames = hits.keys() # Store raw read hit names   
        else:
            hit_readnames=orf_hit_readnames
        
        hit_reads_fasta, direction_information = self._extract_from_raw_reads(
                                                       hit_reads_fasta,
                                                       hit_readnames,
                                                       unpack.read_file,
                                                       unpack.format(),
                                                       hits
                                                       )

        
        if not hit_readnames:
            hit_read_counts = [0, len(hit_readnames)]
            result = DBSearchResult(None,
                                    search_result,
                                    hit_read_counts,
                                    None)

            return result, direction_information


        if unpack.sequence_type() == 'nucleotide':
            # Extract the orfs of these reads that hit the original search
            self._extract_orfs(
                               hit_reads_fasta,
                               extracting_orfm,
                               orf_hit_readnames,
                               hit_reads_orfs_fasta,
                               search_method,
                               list(search_result[0].each([SequenceSearchResult.QUERY_ID_FIELD,
                                                           SequenceSearchResult.ALIGNMENT_DIRECTION,
                                                           SequenceSearchResult.QUERY_FROM_FIELD,
                                                           SequenceSearchResult.QUERY_TO_FIELD])
                                    )
                               )
            
            hit_reads_fasta = hit_reads_orfs_fasta
        slash_endings=self._check_for_slash_endings(hit_readnames)
        result = DBSearchResult(hit_reads_fasta,
                                search_result,
                                [0, len([itertools.chain(*hits.values())])],  # array of hits [euk hits, true hits]. Euk hits alway 0 unless searching from 16S
                                slash_endings)  # Any reads that end in /1 or /2     

        if maximum_range:
            n_hits = sum([len(x["strand"]) for x in hits.values()])
        else:
            n_hits = len(hits.keys())
        logging.info("%s read(s) detected" % n_hits)
                
        return result, direction_information