Extract a protein sequence from a PDB file.
    
    Arguments
    =========
    pdb: str
        The path to a PDB file or a gzipped PDB file.  Alternatively, a PDB id 
        to fetch and download from the internet.

    chain: str
        Which chain or chains to extract the sequence from, and how to combine 
        sequences from multiple chains.  Different values for this option can 
        trigger quite different behaviors:

        - If no chain is specified (e.g. chain=None, the default), behavior 
          depends on how many chains are in the file.  If there's only one, its 
          sequence will be returned.  Otherwise a ValueError is raised and a 
          chain must be specified.

        - If a single chain is specified (e.g. chain="A"), the sequence of that 
          chain will be returned.

        - If multiple chains are joined by "+" signs (e.g. chain="A+B"), the 
          sequences for those chains will be concatenated together in the given 
          order.  This effectively combines the chains into one.
          
        - If multiple chains are joined by "~" signs (e.g. chain="A~B"), the 
          indicated chains will be considered as multiple copies of the same 
          protein, and mutations from every chain will be merged into one 
          sequence.  This behavior requires that the `ref` parameter is given.

        Note that mixing "+" and "~" is not currently allowed.

    ref: str
        The 

    Returns
    =======
    seq: str
        The extracted protein sequence (as a string of one-letter codes).
          
    Note that this function determines sequence from the ATOM records, not the 
    SEQRES records.  This means missing residues will not be included in the sequence.

def sequence_from_pdb(pdb, chain=None, ref=None):
    """
    Extract a protein sequence from a PDB file.
    
    Arguments
    =========
    pdb: str
        The path to a PDB file or a gzipped PDB file.  Alternatively, a PDB id 
        to fetch and download from the internet.

    chain: str
        Which chain or chains to extract the sequence from, and how to combine 
        sequences from multiple chains.  Different values for this option can 
        trigger quite different behaviors:

        - If no chain is specified (e.g. chain=None, the default), behavior 
          depends on how many chains are in the file.  If there's only one, its 
          sequence will be returned.  Otherwise a ValueError is raised and a 
          chain must be specified.

        - If a single chain is specified (e.g. chain="A"), the sequence of that 
          chain will be returned.

        - If multiple chains are joined by "+" signs (e.g. chain="A+B"), the 
          sequences for those chains will be concatenated together in the given 
          order.  This effectively combines the chains into one.
          
        - If multiple chains are joined by "~" signs (e.g. chain="A~B"), the 
          indicated chains will be considered as multiple copies of the same 
          protein, and mutations from every chain will be merged into one 
          sequence.  This behavior requires that the `ref` parameter is given.

        Note that mixing "+" and "~" is not currently allowed.

    ref: str
        The 

    Returns
    =======
    seq: str
        The extracted protein sequence (as a string of one-letter codes).
          
    Note that this function determines sequence from the ATOM records, not the 
    SEQRES records.  This means missing residues will not be included in the sequence.
    """

    seqs = chains_from_pdb(pdb)
    
    if chain is None:
        if len(seqs) == 1:
            return next(iter(seqs.values()))
        else:
            raise scripting.UserError("'{}' contains multiple chains, but no chain was specified.".format(pdb))

    single = re.match('^([A-Z])$', chain)
    merge = re.match('^(?:([A-Z])[~])+([A-Z])$', chain)
    concat = re.match('^(?:([A-Z])[+])+([A-Z])$', chain)
    
    if single:
        chain = single.group(1)
        return seqs[chain]

    if merge:
        if ref is None:
            raise ValueError('template required to merge sequences')

        mutants = [seqs[k] for k in merge.groups()]
        return merge_sequences(ref, *mutants)

    if concat:
        return ''.join(seqs[k] for k in concat.groups())

    raise ValueError("chain '{}' not understood".format(chain))