def cleanup_methyl_protons(labels, Gs, envs, shifts, errs, ws, crysts, inds, atoms, bonds):
    """
    Gather methyl proton shifts into only one per methyl group
    
    Inputs:     - labels        List of labels for each distribution
                - Gs            Graph of each distribution
                - envs          Environment of each graph
                - shifts        Predicted shifts in each distribution
                - errs          Prediction errors in each distribution
                - ws            Maximum depth for each distribution
                - crysts        Crystals in each distribution
                - inds          Indices of the atoms in each distribution
                - atoms         List of atoms in the molecule
                - bonds         Bonded atoms for each atom in the molecule (by index)
    
    Outputs:    - new_labels    Cleaned list of labels for each distribution
                - new_Gs        Cleaned list of graphs
                - new_shifts    Cleaned list of predicted shifts in each distribution
                - new_errs      Cleaned list of prediction errors in each distribution
                - new_ws        Cleaned list of maximum depth for each distribution
                - new_crysts    Cleaned array of crystals in each distribution
                - new_inds      Cleaned array of indices of the atoms in each distribution
    """
    
    # Initialize new arrays
    new_labels = []
    new_Gs = []
    new_envs = []
    new_shifts = []
    new_errs = []
    new_ws = []
    new_crysts = []
    new_inds = []
    
    # Array to store already identified methyl groups
    methyls = []
    
    # Loop over all graphs
    for l, G, env, sh, er, w, cryst, ind in zip(labels, Gs, envs, shifts, errs, ws, crysts, inds):
        
        # Get the index of the central node
        i = G.nodes[0]["ind"]
        # Check that there is only one neighbour (H should be linked to only one atom
        if len(bonds[i]) == 1:
            j = bonds[i][0]
            
            # Check if there are at least three protons linked to the neighbour
            nH = [atoms[k] for k in bonds[j]].count("H")
            
            if nH >= 3:
                # If we find a new methyl proton, add it,
                #   otherwise skip it if it is part of an already known methyl
                if j not in methyls:
                    methyls.append(j)
                    new_labels.append(l)
                    new_Gs.append(G)
                    new_envs.append(env)
                    new_shifts.append(sh)
                    new_errs.append(er)
                    new_ws.append(w)
                    new_crysts.append(cryst)
                    new_inds.append(ind)
            
            # If this is not a methyl proton, add it
            else:
                new_labels.append(l)
                new_Gs.append(G)
                new_envs.append(env)
                new_shifts.append(sh)
                new_errs.append(er)
                new_ws.append(w)
                new_crysts.append(cryst)
                new_inds.append(ind)
        
        # If the proton is bonded to more than one atom, add it
        else:
            new_labels.append(l)
            new_Gs.append(G)
            new_envs.append(env)
            new_shifts.append(sh)
            new_errs.append(er)
            new_ws.append(w)
            new_crysts.append(cryst)
            new_inds.append(ind)
    
    return new_labels, new_Gs, new_envs, new_shifts, new_errs, new_ws, new_crysts, new_inds