<jupyter_start><jupyter_text>Decomposition Example First step, set the paths to where to find the motif dictionary and associated code. Note that these are available at http://github.com/sdrogers/motifdb<jupyter_code>motifdbcodepath = '/Users/simon/git/motifdb/code/utilities/'
motifdbpath = '/Users/simon/git/motifdb/motifs/'<jupyter_output><empty_output><jupyter_text>Load some motifs from the motif database, I'm just loading the massbank ones<jupyter_code>import sys
sys.path.append(motifdbcodepath)
from motifdb_loader import load_db
db_list = ['gnps_binned_005']
spectra,motif_metadata = load_db(db_list,motifdbpath)<jupyter_output>Looking in /Users/simon/git/motifdb/motifs/gnps_binned_005/*.m2m
	 Found 92
Found total of 92 motif files<jupyter_text>spectra is a dictionary with key the filename of the motif and values another dictionary with feature / probability pairs<jupyter_code>import numpy as np<jupyter_output><empty_output><jupyter_text>make an index of the unique features -- we can only decompose onto features that are loaded from the database<jupyter_code>word_index = {}
word_pos = 0
for motif,word_probs in spectra.items():
    for word,probability in word_probs.items():
        if not word in word_index:
            word_index[word] = word_pos
            word_pos += 1<jupyter_output><empty_output><jupyter_text>Put the words and motifs into a matrix ($\beta$)<jupyter_code># Create a beta matrix
motif_index = {}
motif_pos = 0
n_motifs = len(spectra)
n_words = len(word_index)
beta = np.zeros((n_motifs,n_words),np.double)
for motif,word_probs in spectra.items():
    motif_index[motif] = motif_pos
    for word,probability in word_probs.items():
        beta[motif_pos,word_index[word]] = probability
    motif_pos += 1
import pylab as plt
%matplotlib inline
plt.imshow(beta,aspect='auto')<jupyter_output><empty_output><jupyter_text>Renomalise beta as when we save motifs we threshold on the probabilities and they end up not summing to 1*Addition*: made sure there were no zeros<jupyter_code># find the minimum value of beta that isn't zero
    pos = np.where(beta > 0)
    min_val = beta[pos].min()
    zpos = np.where(beta == 0)
    beta[zpos] = min_val/100
beta /= beta.sum(axis=1)[:,None]<jupyter_output><empty_output><jupyter_text>Load some data..I'm loading massbank data, but other datasets are available. Also assuming the motifs were made using binned features<jupyter_code>ldacodepath = '/Users/simon/git/lda/code/'
sys.path.append(ldacodepath)
from ms2lda_feature_extraction import LoadGNPS,MakeBinnedFeatures

import glob
massbank_data = glob.glob('/Users/simon/Dropbox/BioResearch/Meta_clustering/MS2LDA/fingerid-104-traindata/spectra_massbank/*.ms')
sub_data = massbank_data[:100]
l = LoadGNPS()
ms1,ms2,spectral_metadata = l.load_spectra(sub_data)
m = MakeBinnedFeatures()
corpus,word_mz_range = m.make_features(ms2)
corpus = corpus[corpus.keys()[0]]<jupyter_output>Processed 100 spectra
100 documents
After removing empty words, 3695 words left<jupyter_text>New code to load from the MGF<jupyter_code>ldacodepath = '/Users/simon/git/lda/code/'
sys.path.append(ldacodepath)
from ms2lda_feature_extraction import LoadMGF,MakeBinnedFeatures

mgf_file = '/Users/simon/git/lda/notebooks/clusters.mgf'
l = LoadMGF()
ms1,ms2,spectral_metadata = l.load_spectra([mgf_file])
m = MakeBinnedFeatures()
corpus,word_mz_range = m.make_features(ms2)
corpus = corpus[corpus.keys()[0]]<jupyter_output>2331 documents
After removing empty words, 75773 words left<jupyter_text>Do the decompositionSome points to note here: - Features in the spectra that are not in the motifs don't add anything to they get skipped. I compute (`proportion_in`) the percentage of intensity that is usable. - This should be taken into account when interpreting the spectra-motif probabilities (`theta`). These values can be interpreted as the proportion of the usable part of the spectra, and not the total. - I also compute overlap score, as it's useful to see how much of the motif is in the spectrum This is a handy function to compute the overlap score<jupyter_code>def compute_overlap(phi_matrix,motif_pos,beta_row,word_index):
    overlap_score = 0.0
    for word in phi_matrix:
        word_pos = word_index[word]
        overlap_score += phi_matrix[word][motif_pos]*beta_row[word_pos]
    return overlap_score<jupyter_output><empty_output><jupyter_text>- This is the actual decomposition.  - The main loop goes over 100 times, this is probably overkill but it jumps out if the total absolute change in `gamma` is less than `1e-6`. - This code doesn't compute alpha, but I can add that if necessary - I get a warning about a `log(0)`. This is due to some entries in beta being 0, and it doesn't cause any problem, but could be removed by setting all the zero values to somethign very small and then re-normalising the beta matrix. - I do some plotting at the end to show how stuff can be exposed - If you wanted to see the decomposition at the individual feature level then you would need to keep hold of `phi_matrix` which gives (for each word), the probabilities over each motif.<jupyter_code>from scipy.special import psi as psi
theta = {}
K = n_motifs
alpha = 1 # will have some effect, but there's no way to set it. Making it low means we'll get sparse solutions
overlap_scores = {}

doc_scans = []
for doc in corpus:
    doc_scans.append((doc,int(spectral_metadata[doc]['scans'])))
    
doc_scans = sorted(doc_scans,key = lambda x: x[1])
# doc_scans = filter(lambda x: x[1]<=130,doc_scans)
doc_list,scan_list = zip(*doc_scans)
for doc in doc_list:
#     print doc,spectral_metadata[doc]['scans']
    # Compute the proportion of this docs intensity that is represented in the motifs
    total_in = 0.0
    total = 0.0
    doc_dict = corpus[doc]
    for word,intensity in doc_dict.items():
        total += intensity
        if word in word_index:
            total_in += intensity
    proportion_in = (1.0*total_in)/total
#     print '\t',proportion_in
    if proportion_in > 0: # Be careful: if there is no overlap between the features in the spectrum and those across the motifs, bad things might happen
        phi_matrix = {}
        for word in doc_dict:
            if word in word_index:
                phi_matrix[word] = None
        gamma = np.ones(K)
        for i in range(100):
            temp_gamma = np.zeros(K) + alpha
            for word,intensity in doc_dict.items():
                if word in word_index:
                    word_pos = word_index[word]
                    if beta[:,word_pos].sum()>0:
                        log_phi_matrix = np.log(beta[:,word_pos]) + psi(gamma)
                        log_phi_matrix = np.exp(log_phi_matrix - log_phi_matrix.max())
                        phi_matrix[word] = log_phi_matrix/log_phi_matrix.sum()
                        temp_gamma += phi_matrix[word]*intensity
            gamma_change = np.sum(np.abs(gamma - temp_gamma))
            gamma = temp_gamma.copy()
            if gamma_change < 1e-6:
                break
        temp_theta = (gamma/gamma.sum()).flatten()
        theta[doc] = {}
        overlap_scores[doc] = {}
        for motif,motif_pos in motif_index.items():
            theta[doc][motif] = temp_theta[motif_pos]
            overlap_scores[doc][motif] = compute_overlap(phi_matrix,motif_pos,beta[motif_pos,:],word_index)

        # print some things!
        tm = zip(theta[doc].keys(),theta[doc].values())
        tm = sorted(tm,key = lambda x: x[1],reverse = True)
        for mo,th in tm[:3]:
            if overlap_scores[doc][mo]>=0.3:
                print doc,spectral_metadata[doc]['scans'],mo,th,overlap_scores[doc][mo],motif_metadata[mo]['ANNOTATION'][:40]
print scan_list<jupyter_output>(1, 3, 7, 8, 9, 11, 12, 18, 19, 20, 22, 40, 44, 45, 46, 49, 55, 84, 88, 92, 93, 94, 96, 102, 103, 104, 105, 106, 107, 110, 113, 118, 119, 120)
{'fraction_non_alphanumeric': 0.09667558272831486, 'fraction_numerical': 0.0211438033371545, 'mean_word_length': 3.9352608422375863, 'pattern_counts': {'":': 0, '<': 39, '<?xml version=': 0, '>': 41, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '249765', 'license': 'MIT', 'n_tokens_mistral': 2704, 'n_tokens_neox': 2500, 'n_words': 756, 'path': 'lda/notebooks/Decomposition notebook.ipynb', 'repo_name': 'RP0001/ms2ldaviz', 'size': '7851'}