docs: 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