def get_genes(exp_file, samples, threshold, max_only):
    """
    Reads in and parses the .bed expression file.
    File format expected to be:
        Whose format is tab seperated columns with header line:
        CHR  START  STOP  GENE  <sample 1>  <sample 2>  ...  <sample n>

    Args:
        exp_file (str): Name of expression file.
        samples (list): Names of the samples in the vcf file.
        threshold (float): Expression threshold to filter lowly/unexpressed genes.
        max_only (bool): if true, gene_dict value is 1 value = max expression
            if false gene_dict value is list of expression values
            YYY: WARNING: if want list to have meaning
                then values needs to be tied to header sample names

    Returns:
        gene_dict (dict): {gene_name: [expression_vals]}.
            Only include values for samples in the vcf.
    """
    data_cols = []
    gene_dict = {}

    print('start read exp_file:' + format(exp_file))

    if max_only:
        # read and only return max exp value in gene_dict
        with open(exp_file) as f:
            header = f.readline().strip().split('\t')
            for samp in header[4:]:
                if samp in samples:
                    data_idx = header.index(samp)
                    data_cols.append(data_idx)

            # Read in expression levels for each gene.
            for line in f:
                line = line.strip().split('\t')
                gene_name = line[3].upper()
                exp_val = -1e1000
                for idx in data_cols:
                    if float(line[idx]) > exp_val:
                        exp_val = float(line[idx])

                gene_dict[gene_name] = exp_val

    else:
        # read and return exp value list in gene_dict
        with open(exp_file) as f:
            header = f.readline().strip().split('\t')
            for samp in header[4:]:
                if samp in samples:
                    data_idx = header.index(samp)
                    data_cols.append(data_idx)

            # Read in expression levels for each gene.
            for line in f:
                line = line.strip().split('\t')
                gene_name = line[3].upper()
                exp_vals = []
                for idx in data_cols:
                    exp_vals.append(line[idx])

                gene_dict[gene_name] = exp_vals

    return gene_dict