Pipeline 2.0

[mbruhns]

Moin,

wie schon angedroht sollten wir nochmal die gesamte Pipeline evaluieren bzw. eventuell anpassen, gab ja einigen interessanten Input bei deiner Verteidigung. Bin jetzt mit meinem MPI-Bewerbungsshizzle durch und daher wieder etwas besser verfügbar.

Wie passt dir Montag 16:00? Können das persönlich oder sonst über Zoom machen.

[mbruhns]

Bonus: wir werden versuchen das Paper beim 20th Annual Cancer Immunotherapy (CIMT) Meeting 2022 einzureichen. Die Deadline für Abstracts ist der 04 März 2022.

[mbruhns]

For reference:

[donEnno]

[mbruhns]

Enthält die Matrix für BLOSUM62 negative Distanzen? 🤨

[donEnno]

Moin!

Kein Stress, nach Prüfungen muss auch mal Ruhe sein!

Schön, das von dir zu hören, ich stress mich da manchmal selber etwas.

Jetzt zur Pipeline: Ich bin heute die Links durchgegangen und habe (denke ich) das Wichtigste verstanden. Jetzt stehe ich vor der Frage, wo ich am schlausten anfange.

Die Pipeline kann:

1. Train/Test Split der DM
2. UMAP
3. Cluster-Vector -> Feature (absolute, relative, frequency, relative frequency als Zählmaße)
4. kNN-Graph aus der trainDM bauen
5. Via kNN Classification (hier bisher ein "Mehrheitsentscheid", also der mode der kNNs) die testDM klassifizieren, daraus einen Feature-Vector bauen
6. Mit logistischer Regression klassifizieren

Was soll rein?

7. low rank sparse composition um die DM zu refinen
8. spectral clustering neben Louvain bzw. Leiden
9. ARI vs Quantile
   10. Selektion der Kohorten (BL/FU/HD)
   11. Patienten-Profile zu Compositional Data

Probleme: Zu 4. Der naive Ansatz

tree = KDTree(train_dm)
dd, ii = tree.query(test_dm, k=k)

funktioniert nicht, für die große DM wegen der Laufzeit. Ich habe auch gelesen, dass ein KDTree für Daten mit hoher Dimensinalität nicht gut sein soll (O(ND), mit D = dimensionality).

Eine Alternative wäre ein BallTree, der für hohe D gut funktionierten soll.

Oder wird dieses Problem schon durch (7) gelöst? Also, dass auf der refinten DM der KDTree, dann wieder gut funktioniert?

Zu 9. Habe ich nicht ganz verstanden, was du meintest mit "auf der x-achse das ausgewählte Quantil und auf y-Achse der ARI". Braucht der ARI nicht eine ground truth, um ein Clustering zu bewerten?

Zum allgemeinen Vorgehen: Vor meinem inneren Auge ist die Pipeline schon fertig. Den genauen Weg dahin, kann ich mir gerade aber nicht ganz vorstellen. Meine Idee ist, erstmal mit 7. anzufangen und darauf aufbauen Louvain für Leiden auszutauschen und das Spectral Clustering zu integrieren. Dann über Verlgeich der Clusterings/gridSearch für gamma bei Louvain bzw. Leiden/über das selftuning SpectralClustering eine sinnvolle Anzahl der Cluster finden. Sobald das steht, möchte ich eine main schreiben, die die Klassifikation für gewünschte Parameter/Clusterings/... performt und die Daten sammelt. Zu guter letzt, muss das ganze natürlich noch visualisiert werden. Da wird Kilian mit Sicherheit dann mitreden wollen.

Was denkst du zu dem Fahrplan? Anmerkungen und Kritik sind wie wimmer gern gesehen. Ich habe auch einen neuen Branch erstellt, der den aktuellsten Code, sowie ein dummy.fasta file beinhaltet.

PS: Das mit den negativen Distanzen in der BLOSUM62 DM überprüfe ich nochmal.

[donEnno]

Moinsen,

die Pipeline steht jetzt soweit, ein paar Schliffe braucht es aber noch.

Das Clustering mit spectral, louvain und leiden funktioniert. Beim spectral clustering habe ich eine eigengap_heuristic Methode geschrieben, um die Zahl der Cluster abzuschätzen. Hier stelle ich auch noch ein Bild zu rein. Die Zahl der cluster kann man aber auch gut "zu Fuß" über die Klassifikations-Performance testen.

Die kNN Klassifizierung funktioniert gut und relativ zackig mit:

nn = NearestNeighbors(n_neighbors=k, metric='precomputed', algorithm='auto')
nn.fit(train_dm)
dd, ii = nn.kneighbors(test_dm, n_neighbors=11)

Hier habe ich mich gefragt, ob man als test_dm die Affinitymatrix oder die Distancematrix nimmt. Ähnlich wie bei UMAP, sollte das ja nur mit der DM, aber das Clustering basiert ja auf der AM.

Die logistische Regression mit 5-fold CV auf den Patienten und (vorerst) ohne den Wald Test mit diversen Performance-Maßen: balanced accuracy, f1, sensitivity, specificity und precision. Das sollte reichen, oder?

Als nächstes werde ich eine Methode schreiben, die das für die gewünschten Settings laufen lassen kann. Dabei kann man dann die ARI scores und die Klassifikations-Performance sammeln.

for sm in substitution_matrices
    for cut_off in k%
        for g in gamma // for n in number_of cluster
            compute_ari_score
            make_classification

Hier kann man ggf. auch mit Parallelisierung arbeiten, damit das fixer geht. Wie das funktioniern kann, habe ich mir aber noch nicht überlegt. Vielleicht können wir hier nochmal besprechen, ob da nicht noch andere Maße, Verteilungen etc. erhoben und was je Iteration visualisiert werden sollte.

Ende der Statusmeldung

[mbruhns]

Sehr stark! Wollen wir nachher kurz via Zoom sprechen? Dann kanst du vielleicht auch Ergebnisse zeigen, falls du welche hast.

[donEnno]

Scheinbar kann man Konsensus-Sequenzen über Biopython ermitteln.

import sys
from Bio import AlignIO
from Bio.Align import AlignInfo

alignment = AlignIO.read(sys.argv[1], 'fasta')
summary_align = AlignInfo.SummaryInfo(alignment)
summary_align.dumb_consensus(float(sys.argv[2]))

Wikiepedia zu Konsensus-Sequenzen:

In molecular biology and bioinformatics, the consensus sequence (or canonical sequence) is the calculated order of most frequent residues, either nucleotide or amino acid, found at each position in a sequence alignment.

Das kann man auch leicht selber machen und ist, vermute ich, was biopython macht. Den UMAP labels steht somit nichts im Wege 😆

[mbruhns]

Perfekt, das wäre in Sachen usability ein nettes Feature. Kannst du bei Gelegenheit noch die ersten Resultate aus der vollen Pipeline für Leiden+Spectral hochladen?

[donEnno]

Unten siehst du die Resultate für die relativen Feature-Vektoren zu allen drei SMs. Für diesen Run habe ich die FUs mal ausgeschlossen, damit alles etwas mehr balanced ist. Das hatte zur Folge, das Leiden auf einmal nur sehr wenige Cluster findet und die eigengap (so wie ich sie implementiert habe) gar nicht taugt. Ich passe noch 1 2 Sachen an, teste ein gutes gamma für die drei SMs, sodass mehr Cluster gefunden werden und lasse den Spaß dann nochmal laufen.

Zur Erklärung: self in der sm Spalte bedeuetet, hier wurde auf den Trainingsdaten klassifiziert.

[mbruhns]

Uff, das sind ja herausragend schlechte Ergebnisse 🤦🏼‍♂️😅 Was ist bei der Analyse der Cluster herausgekommen? Sieht ja so aus, als wären die nicht wirklich informativ.

[donEnno]

Uff

Cluster

Die Cluster sind ähnlich wie schon in der BA. Kein Cluster besteht nur aus einer Klasse, sondern alle haben einen Mix aus BL und HD Anteilen, der sich verändert. Ich würde erwarten, ähnlich wie in der BA, einen Cut-Off im HD:BL Ratio zu sehen, wenn man das auflistet.

Ich lasse gerade mal laufen, wie sich die Zahl der Cluster auf die Klassifikation auswirkt. Unten die Ergebnisse für die drei SMs mit gamma = 1 bis 1.2 in 0.01er Schritten. Oben immer die Klassifikation auf den Trainingsdaten und unten auf den Testdaten. Der kNN Classifier macht hier alles mit k=11, ich hab aber spaßeshalber mal k=111 und k=1111 laufen lassen, und das sah sehr ähnlich aus. Schwierig zu sagen, was genau hier passiert.. Am anfang sieht es nach einem harten "underfit" aus und am Ende nach einem harten "overfit", oder was meinst du?

Dasselbe mache ich heute noch für Spectral Clustering klar, mit ähnlichen Werten für die Anzahl der Cluster.

PAM70

BLOSUM45

BLOSUM62

Hier sieht man ganz gut, wie "unberechenbar" das gamma ist. Während bei B45 fast 2000 cluster gefunden werden, gibts bei B62 nicht mal 500.

[donEnno]

Hier die versprochenen Plots für das spectral-clustering:

PAM70

BLOSUM45

BLOSUM62

[donEnno]

Update

Info zu den Plots

Reihenfolge der Bilder is in allen drei Fällen:

1. Leiden: Klassifikation auf Trainingsdaten
2. Leiden: Klassifikation auf Testdaten
3. Spectral: Klassifikation auf Trainingsdaten
4. Spectral: Klassifikation auf Testdaten

Insgesamt ist "Gut" immer noch was anderes, aber ich bin schonmal erleichtert, dass ein bisschen mehr Dynamik in den Ergebnissen ist. Auch finde ich, dass bei den Spectral-Clusterings eine Tendenz zu besseren Ergebnissen zu erkennen ist.

Die nächste Analyse

Im nächsten Schritt würde ich gerne zwei Dinge ausprobieren:

• nochmal die Auswirkung des k%-Pruning anschauen
• mal den Run, der jetzt gelaufen ist, mit "absolute" und "freq" anschauen

Zum zweiten Punkt: Hier nochmal die Tabelle die ich im Januar erzeugt habe. Das war noch die alte Pipeline, aber die Ergebnisse für das "absolute" Zählmaß waren, mMn, so gut, dass man das mal ausprobieren sollte.

ARI: Ground Truth vs. Train-Test-Combined

Und noch eine Info am Rande, auch wenn ich hier keine Plots erzeugt habe. Ich habe gestern für den gesamten Datensatz ein Clustering erzeugt. Dann gesplittet, auf den Trainingsdaten geclustert und anschließend

• einen Dummy-Array aus Nullen appended
• die kNN Klassifikation das Clustering "auffüllen" lassen und anschließend jeweils den ARI berechnet. Bei dem dummy-Array kam ein Wert um 0.4 raus, während der kNN-Ansatz einen ARI von 0.87 geschafft hat. D.h. dass der kNN Ansatz zumindest nicht untauglich ist.

PAM70

BLOSUM45

BLOSUM62

[mbruhns]

Okay, also immerhin passiert jetzt was im Training. Besteht weiterhin das Problem, dass im Test die Performance quasi bei Null liegt 🤦🏼‍♂️ Bzw. scheint ja nur BLOSUM62+Spectral irgendwas zu lernen. In Sachen Implementierung bist du sicher, dass jetzt alles richtig läuft?

• die kNN Klassifikation das Clustering "auffüllen" lassen und anschließend jeweils den ARI berechnet. Bei dem dummy-Array kam ein Wert um 0.4 raus, während der kNN-Ansatz einen ARI von 0.87 geschafft hat. D.h. dass der kNN Ansatz zumindest nicht untauglich ist.

Gute Idee das so einmal zu validieren! Hast du die Performance für verschiedene k verglichen?

Zu #28 : Bevor wir andere Klassifizierer ins Spiel bringen, sollte die durchschnittliche Bal_Acc noch deutlich höher sein.

Mich verwundert es auch wenig, dass die Performance im Training scheinbar mit größerer Anzahl an Clustern monoton zunimmt, ich hätte da ein klareres Optimum erwartet. Das könnte in Kombination mit den anderen Resultaten aber auch ein Indiz für Overfitting sein.

Den aktuellen Stand verstehe ich so, dass das Modell prinzipiell etwas lernt (--> Trainings-Performance), auf neuen Daten bis auf eine Ausnahme (BLOSUM62+Spectral) aber nur rät. Daraus leite ich zwei Möglichkeiten ab: Overfitting und/oder Fehler in der Test-Pipeline.

---

Am besten machen wir uns beide nochmal prinzipiell Gedanken über das Ganze und quatschen morgen/übermorgen via Zoom.

[donEnno]

Moin!

In Sachen Implementierung bist du sicher, dass jetzt alles richtig läuft?

Ja, ziemlich sicher. Der kNN-Classifier sind nur wenige Zeilen und laut dem Initialen-ARI-Test, scheint das auch was zu taugen. Das Bauen der Feature habe ich mehr oder weniger 1:1 aus der BA übernommen und dort hatte ich mir das gut überlegt. Ich bleibe aber wachsam 😄

Gute Idee das so einmal zu validieren! Hast du die Performance für verschiedene k verglichen?

Nein. Tagesziel heute ist das Skript für den k-run zu schreiben, so dass er gleichzeitig ARIs berechnet.

Update zu dem absolute Durchlauf

Was mich verwundert ist, dass der monotone Anstieg im Training, im Test gar nicht zu sehen ist. Und wieder ist bei B62+spectral das beste Model dabei.

PAM70

BLOSUM45

BLOSUM62

[donEnno]

Moin!

Hier die versprochenen Plots für die verbleibenden k% gepruned:

Ich lasse es gerade nochmal laufen mit Timer und einem Bugfix. Das ist jetzt nach wie vor ohne CV. Du meintest du hattest da eine Idee? Ansonsten würde ich nach dem aktuellen Run, eine klassische 5-Fold CV auf den Patienten laufen lassen und dann nochmal die Plots teilen. Außerdem werde ich am WE versuchen die Transformation und die Cross-Features zu implementieren.

LG und bis Montag!

[donEnno]

5-Fold CV Ergebnisse

Laufzeit: 70h

(Ich heule grade in bisschen :D)

[mbruhns]

(Ich heule grade in bisschen :D)

Das wird schon 😅 Wie hast du den plot erstellt? Mir fehlen da irgendwie die Fehlerbalken 😄

[donEnno]

Hier noch das Zöliakie Paper. Die haben wie gesagt basically dasselbe gemacht wie wir, nur eben mit Zöliakie-Patienten und dem k-mer Ansatz. Wir benutzen aber nur die delta-Ketten Info, während die auch die gamma-Kette miteinbeziehen. LG ✌🏼