original_docstring: Generate datasets for all features in hdf5 file.
      Args:
          h5          - h5 filename
          gene_annot  - Gene annotation bed filename
          reg_bound   - Tuple for the boundary of regulatory region, i.e.
                          (upstream distance, downstream distance)
          feat_dict   - Dictionary of feature types and feature names
      Returns:
          NULL
      
docstring: Generate datasets for all features in hdf5 file.
  Args:
  h5          - h5 filename
  gene_annot  - Gene annotation bed filename
  reg_bound   - Tuple for the boundary of regulatory region, i.e.
  (upstream distance, downstream distance)
  feat_dict   - Dictionary of feature types and feature names
  Returns:
  NULL
docstring_tokens: Generate datasets for all features in hdf5 file . Args : h5 - h5 filename gene_annot - Gene annotation bed filename reg_bound - Tuple for the boundary of regulatory region i . e . ( upstream distance downstream distance ) feat_dict - Dictionary of feature types and feature names Returns : NULL
code: def generate_features(h5, gene_annot, reg_bound, feat_dict):
      """Generate datasets for all features in hdf5 file.
      Args:
          h5          - h5 filename
          gene_annot  - Gene annotation bed filename
          reg_bound   - Tuple for the boundary of regulatory region, i.e.
                          (upstream distance, downstream distance)
          feat_dict   - Dictionary of feature types and feature names
      Returns:
          NULL
      """
      FEAT_COL = ['gene_idx', 'mtx_start', 'mtx_end', 'peak_score']
  
      ## Parse gene list and regulator DNA
      gene_bed = BedTool(gene_annot)
      genes = sorted(pd.unique(gene_bed.to_dataframe()['name']))
      regdna_bed = create_regdna(gene_bed, feat_dict['dna_sequence'], reg_bound)
  
      ## Write grouped datasets into hdf5
      with h5py.File(h5, 'w') as f:
          ## Gene list 
          f.create_dataset(
              'genes', 
              data=np.array(genes, dtype='S'), 
              compression='gzip')
  
          for k1, v1 in feat_dict.items(): 
              ## Regulatory DNA sequence
              if k1 == 'dna_sequence': 
                  logger.info('Working on {} {}'.format(k1, v1))
                  ## Store one-hot encode sequence in A, C, G, T channels
                  # logger.debug('... one-hot DNA')
                  regdna_df = get_onehot_dna_sequence(regdna_bed, v1, reg_bound, genes)
                  g = f.require_group(k1)
                  for i, alphabet in enumerate(['A', 'C', 'G', 'T']):
                      g.create_dataset(
                          alphabet, 
                          data=regdna_df.loc[regdna_df['alphabet'] == i, FEAT_COL].values.astype(int), 
                          compression='gzip')
                  # Store single and di-nucleotide frequencies
                  # logger.debug('... di-nucleotide freqs')
                  nt_freq_dict = get_nt_frequency(regdna_bed, v1, genes)
                  g = f.require_group(k1 + '_nt_freq')
                  for nt, freq_arr in nt_freq_dict.items():
                      g.create_dataset(
                          nt,
                          data=np.array(freq_arr, dtype='float16'),
                          compression='gzip')
  
              ## Gene expression level features
              elif k1 == 'gene_expression' or k1 == 'gene_variation': 
                  for k2, v2 in v1.items():
                      logger.info('Working on {} {}'.format(k1, k2))
                      g = f.require_group('gene_expression')
                      ## Load gene expression matrix and sort genes in the 
                      ## same dimension as other features
                      expr_df = pd.read_csv(v2, index_col=0)
                      expr_df = expr_df.loc[genes]
                      if k1 == 'gene_variation':
                          g.create_dataset(
                              'variation',
                              data=expr_df.values[:, 0].astype(float),
                              compression='gzip')
                      else:
                          ## Store individual gene expression profile
                          for sample in sorted(expr_df.columns):
                              g.create_dataset(
                                  sample,
                                  data=expr_df[sample].values.astype(float),
                                  compression='gzip')
  
              ## Epigenetic features, e.g. TF binding data
              else: 
                  g = f.require_group(k1)
                  for k2, v2 in v1.items():  ## Feature name, e.g. TF name
                      logger.info('Working on {} {}'.format(k1, k2))
                      v2_bed = BedTool(v2)
                      ## Intersect peak features with regulatory region
                      feat_df = intersect_peak_regdna(v2_bed, regdna_bed, gene_bed.to_dataframe())
                      feat_df = calculate_matrix_position(feat_df, reg_bound[0])
                      feat_df['gene_idx'] = [genes.index(x) for x in feat_df['gene']]
                      g.create_dataset(
                          k2, 
                          data=feat_df[FEAT_COL].values.astype(float), 
                          compression='gzip')
code_tokens: def generate_features ( h5 , gene_annot , reg_bound , feat_dict ) : FEAT_COL = [ 'gene_idx' , 'mtx_start' , 'mtx_end' , 'peak_score' ] gene_bed = BedTool ( gene_annot ) genes = sorted ( pd . unique ( gene_bed . to_dataframe ( ) [ 'name' ] ) ) regdna_bed = create_regdna ( gene_bed , feat_dict [ 'dna_sequence' ] , reg_bound ) with h5py . File ( h5 , 'w' ) as f : f . create_dataset ( 'genes' , data = np . array ( genes , dtype = 'S' ) , compression = 'gzip' ) for k1 , v1 in feat_dict . items ( ) : if k1 == 'dna_sequence' : logger . info ( 'Working on {} {}' . format ( k1 , v1 ) ) regdna_df = get_onehot_dna_sequence ( regdna_bed , v1 , reg_bound , genes ) g = f . require_group ( k1 ) for i , alphabet in enumerate ( [ 'A' , 'C' , 'G' , 'T' ] ) : g . create_dataset ( alphabet , data = regdna_df . loc [ regdna_df [ 'alphabet' ] == i , FEAT_COL ] . values . astype ( int ) , compression = 'gzip' ) nt_freq_dict = get_nt_frequency ( regdna_bed , v1 , genes ) g = f . require_group ( k1 + '_nt_freq' ) for nt , freq_arr in nt_freq_dict . items ( ) : g . create_dataset ( nt , data = np . array ( freq_arr , dtype = 'float16' ) , compression = 'gzip' ) elif k1 == 'gene_expression' or k1 == 'gene_variation' : for k2 , v2 in v1 . items ( ) : logger . info ( 'Working on {} {}' . format ( k1 , k2 ) ) g = f . require_group ( 'gene_expression' ) expr_df = pd . read_csv ( v2 , index_col = 0 ) expr_df = expr_df . loc [ genes ] if k1 == 'gene_variation' : g . create_dataset ( 'variation' , data = expr_df . values [ : , 0 ] . astype ( float ) , compression = 'gzip' ) else : for sample in sorted ( expr_df . columns ) : g . create_dataset ( sample , data = expr_df [ sample ] . values . astype ( float ) , compression = 'gzip' ) else : g = f . require_group ( k1 ) for k2 , v2 in v1 . items ( ) : logger . info ( 'Working on {} {}' . format ( k1 , k2 ) ) v2_bed = BedTool ( v2 ) feat_df = intersect_peak_regdna ( v2_bed , regdna_bed , gene_bed . to_dataframe ( ) ) feat_df = calculate_matrix_position ( feat_df , reg_bound [ 0 ] ) feat_df [ 'gene_idx' ] = [ genes . index ( x ) for x in feat_df [ 'gene' ] ] g . create_dataset ( k2 , data = feat_df [ FEAT_COL ] . values . astype ( float ) , compression = 'gzip' )
short_docstring: Generate datasets for all features in hdf5 file.
short_docstring_tokens: Generate datasets for all features in hdf5 file .
comment: """Generate datasets for all features in hdf5 file.
      Args:
          h5          - h5 filename
          gene_annot  - Gene annotation bed filename
          reg_bound   - Tuple for the boundary of regulatory region, i.e.
                          (upstream distance, downstream distance)
          feat_dict   - Dictionary of feature types and feature names
      Returns:
          NULL
      """$SEP$## Parse gene list and regulator DNA$SEP$## Write grouped datasets into hdf5$SEP$## Gene list $SEP$## Regulatory DNA sequence$SEP$## Store one-hot encode sequence in A, C, G, T channels$SEP$# logger.debug('... one-hot DNA')$SEP$# Store single and di-nucleotide frequencies$SEP$# logger.debug('... di-nucleotide freqs')$SEP$## Gene expression level features$SEP$## Load gene expression matrix and sort genes in the $SEP$## same dimension as other features$SEP$## Store individual gene expression profile$SEP$## Epigenetic features, e.g. TF binding data$SEP$## Feature name, e.g. TF name$SEP$## Intersect peak features with regulatory region