docs: def create_features(bam_in, loci_file, reference, out_dir):
      """
      Use feature extraction module from CoRaL
      """
      lenvec_plus = op.join(out_dir, 'genomic_lenvec.plus')
      lenvec_minus = op.join(out_dir, 'genomic_lenvec.minus')
      compute_genomic_cmd = ("compute_genomic_lenvectors "
                             "{bam_in} {lenvec_plus} "
                             "{lenvec_minus} "
                             "{min_len} "
                             "{max_len} ")
      index_genomic_cmd = ("index_genomic_lenvectors "
                           "{lenvec} ")
      genomic_lenvec = op.join(out_dir, 'genomic_lenvec')
      feat_len_file = op.join(out_dir, 'feat_lengths.txt')
      compute_locus_cmd = ("compute_locus_lenvectors "
                           "{loci_file} "
                           "{genomic_lenvec} "
                           "{min_len} "
                           "{max_len} "
                           "> {feat_len_file}")
      cov_S_file = op.join(out_dir, 'loci.cov_anti')
      coverage_anti_cmd = ("coverageBed -S -counts -b "
                           "{bam_in} -a {loci_file} "
                           "> {cov_S_file}")
      feat_posentropy = op.join(out_dir, 'feat_posentropy.txt')
      entropy_cmd = ("compute_locus_entropy.rb "
                     "{counts_reads} "
                     "> {feat_posentropy}")
      with utils.chdir(out_dir):
          run(compute_genomic_cmd.format(min_len=min_trimmed_read_len, max_len=max_trimmed_read_len, **locals()), "Run compute_genomic")
          run(index_genomic_cmd.format(lenvec=lenvec_plus), "Run index in plus")
          run(index_genomic_cmd.format(lenvec=lenvec_minus), "Run index in minus")
          run(compute_locus_cmd.format(min_len=min_trimmed_read_len, max_len=max_trimmed_read_len, **locals()), "Run compute locus")
          run(coverage_anti_cmd.format(**locals()), "Run coverage antisense")
          feat_antisense = _order_antisense_column(cov_S_file, min_trimmed_read_len)
  
          counts_reads = _reads_per_position(bam_in, loci_file, out_dir)
          run(entropy_cmd.format(**locals()), "Run entropy")
  
          rnafold = calculate_structure(loci_file, reference)