def identify_native_substructures(native_file, d_cutoff, min_seq_separation, contact_sep_thresh, min_clustersize, atom='CA', labelsize=30, fontsize=30, max_res=None, plot=True, ax=None, native_contacts=[], verbose=False):
    '\n    Identify substructures within native file contact map\n    Using the following strategy\n    \n    We produce a contact map which is a bunch of dots\n    Contacts correspond to pairs of residues that are less than d_cutoff apart\n    6 Angstroms is generally a good value, but may want a larger value for helical proteins where residues interact primarily\n    via sidechains, and thus the alpha carbons are further apart\n    \n    We only count contacts if the residues are separated by min_seq_separation along the primary sequence\n    We set min_seq_separation relatively high because we don\'t care to include intra-helix contacts within our contact map\n    \n    Ce can calculate the "Manhattan" distnace between every pair of dots on that contact map \n    and build a graph of contacts in which two contacts are connected if they are less than some\n    threshold distance, contact_sep_thresh, apart in the contact map \n    \n    Then, we find all connected components of this graph, each of which is a substructure\n    But we only keep substructures whose number of residues is at least min_clustersize, to avoid sparse contacts here and there that we dont\' care about\n    \n    Gives you option to input native contacts a priori, but by defualt you don\'t do this (value is set to None)\n    \n    You set Max_res to something other than None if you wish to plot only up to a certain residue number (ex. to depict what substructures can be formed when the first 100 AAs are synthesized)\n    '
    if (len(native_contacts) == 0):
        (coords, resis) = read_PDB(native_file, atom)
        native_distances = compute_contacts_matrix(coords, mode='distances', min_seq_separation=min_seq_separation)
        native_contacts = np.zeros(np.shape(native_distances))
        native_contacts[np.where(((native_distances < d_cutoff) & (native_distances != 0)))] = 1
    positions = np.where((native_contacts == 1))
    positions = np.transpose(positions)
    M = metrics.pairwise.pairwise_distances(positions, metric='manhattan')
    (clusters, pairs_in_substructures, mean_intercluster, mean_intracluster) = loopCluster(contact_sep_thresh, positions, M, sort_orphans=False, min_clustersize=min_clustersize, verbose=verbose)
    pairs_in_substructures = [(np.array(C)[:, 0], np.array(C)[:, 1]) for C in pairs_in_substructures]
    nsubstructures = len(pairs_in_substructures)
    substructures = np.zeros((np.shape(native_contacts)[0], np.shape(native_contacts)[1], nsubstructures))
    for n in range(nsubstructures):
        SS = np.zeros(np.shape(native_contacts))
        SS[pairs_in_substructures[n]] = 1
        substructures[:, :, n] = SS
    if plot:
        visualize_substructures(native_contacts, substructures, max_res=max_res, ax=ax, labelsize=labelsize, fontsize=fontsize)
    return (native_contacts, substructures)
Identify substructures within native file contact map
Using the following strategy

We produce a contact map which is a bunch of dots
Contacts correspond to pairs of residues that are less than d_cutoff apart
6 Angstroms is generally a good value, but may want a larger value for helical proteins where residues interact primarily
via sidechains, and thus the alpha carbons are further apart

We only count contacts if the residues are separated by min_seq_separation along the primary sequence
We set min_seq_separation relatively high because we don't care to include intra-helix contacts within our contact map

Ce can calculate the "Manhattan" distnace between every pair of dots on that contact map 
and build a graph of contacts in which two contacts are connected if they are less than some
threshold distance, contact_sep_thresh, apart in the contact map 

Then, we find all connected components of this graph, each of which is a substructure
But we only keep substructures whose number of residues is at least min_clustersize, to avoid sparse contacts here and there that we dont' care about

Gives you option to input native contacts a priori, but by defualt you don't do this (value is set to None)

You set Max_res to something other than None if you wish to plot only up to a certain residue number (ex. to depict what substructures can be formed when the first 100 AAs are synthesized)
def identify_native_substructures(native_file, d_cutoff, min_seq_separation, contact_sep_thresh, min_clustersize, atom='CA', labelsize=30, fontsize=30, max_res=None, plot=True, ax=None, native_contacts=[], verbose=False):
    '\n    Identify substructures within native file contact map\n    Using the following strategy\n    \n    We produce a contact map which is a bunch of dots\n    Contacts correspond to pairs of residues that are less than d_cutoff apart\n    6 Angstroms is generally a good value, but may want a larger value for helical proteins where residues interact primarily\n    via sidechains, and thus the alpha carbons are further apart\n    \n    We only count contacts if the residues are separated by min_seq_separation along the primary sequence\n    We set min_seq_separation relatively high because we don\'t care to include intra-helix contacts within our contact map\n    \n    Ce can calculate the "Manhattan" distnace between every pair of dots on that contact map \n    and build a graph of contacts in which two contacts are connected if they are less than some\n    threshold distance, contact_sep_thresh, apart in the contact map \n    \n    Then, we find all connected components of this graph, each of which is a substructure\n    But we only keep substructures whose number of residues is at least min_clustersize, to avoid sparse contacts here and there that we dont\' care about\n    \n    Gives you option to input native contacts a priori, but by defualt you don\'t do this (value is set to None)\n    \n    You set Max_res to something other than None if you wish to plot only up to a certain residue number (ex. to depict what substructures can be formed when the first 100 AAs are synthesized)\n    '
    if (len(native_contacts) == 0):
        (coords, resis) = read_PDB(native_file, atom)
        native_distances = compute_contacts_matrix(coords, mode='distances', min_seq_separation=min_seq_separation)
        native_contacts = np.zeros(np.shape(native_distances))
        native_contacts[np.where(((native_distances < d_cutoff) & (native_distances != 0)))] = 1
    positions = np.where((native_contacts == 1))
    positions = np.transpose(positions)
    M = metrics.pairwise.pairwise_distances(positions, metric='manhattan')
    (clusters, pairs_in_substructures, mean_intercluster, mean_intracluster) = loopCluster(contact_sep_thresh, positions, M, sort_orphans=False, min_clustersize=min_clustersize, verbose=verbose)
    pairs_in_substructures = [(np.array(C)[:, 0], np.array(C)[:, 1]) for C in pairs_in_substructures]
    nsubstructures = len(pairs_in_substructures)
    substructures = np.zeros((np.shape(native_contacts)[0], np.shape(native_contacts)[1], nsubstructures))
    for n in range(nsubstructures):
        SS = np.zeros(np.shape(native_contacts))
        SS[pairs_in_substructures[n]] = 1
        substructures[:, :, n] = SS
    if plot:
        visualize_substructures(native_contacts, substructures, max_res=max_res, ax=ax, labelsize=labelsize, fontsize=fontsize)
    return (native_contacts, substructures)