def readDatFile(dat_file,take_only_top=20):
    """
    Read a codeml dat file.

    take_only_top is an integer that forces the program to only sample the top
    N states.  By default it is set to 20, so all states will be sampled.
    """

    f = open(dat_file,'r')
    lines = f.readlines()
    f.close()

    # Remove blank lines and commented (#) lines
    lines = [l for l in lines if l[0] != "#" and l.strip() != ""]

    # Read the amino acids and posterior probabilities for each site into a
    # list.
    aa_list = []
    pp_list = []
    for l in lines:
        aa_list.append([])
        pp_list.append([])
        columns = l.split()[1:]
        for i in range(0,len(columns),2):

            if take_only_top and i > (2*take_only_top):
                break
                
            aa_list[-1].append(columns[i])
            pp_list[-1].append(float(columns[i+1]))
            
        # Normalize the values in pp to be between 0 and 1
        pp_list[-1] = [p/sum(pp_list[-1]) for p in pp_list[-1]]


    # Combine amino acid and posterior probability data into a single list
    anc_data = [aa_list,pp_list]

    return anc_data