ich arbeite in der klinischen Humangenetik und beschäftige mich zurzeit damit, aus Sequenzierdaten einen Zugewinn oder Verlust von genomischen Material zu detektieren. Dazu wird das gesamte Genom mit sehr niedriger Abdeckung sequenziert und die Anzahl der Sequenzen, welche pro Chromosom zugeordnet werden, gezählt. Der Hintergrund ist einfach: Ein Verlust eines Chromosoms 3 sollte dazu führen, das für dieses Chromosom nur noch ungefähr halb soviel Sequenzen gefunden werden können. Für Zugewinne entsprechend mehr.
Ich hoffe, dass ich für euch ein interessantes Beispiel mitgebracht habe. Zum einen interessiert mich, ob die verwendete Methodik korrekt ist und zum anderen wäre ich über eine Diskussion zum zweiten Punkt (am Ende des Beitrags) interessiert.
Die Auswertung erfolgt im Moment so:
Die gezählten Sequenzen pro Chromosom (raw_read_count) werden durch die Gesamtanzahl der Sequenzen im Sample (alle Chromosomen) dividiert (=relative_read_count). Anschließend wird der zuvor ermittelte Wert (relative_read_count) durch den Mittelwert der Sequenzen des entsprechenden Chromosoms aller Kontrollsamples (mean_relative_read_count) dividiert (=normalized_relative_read_count). Der ZScore wird wie folgt berechnet: (relative_read_count - mean_relative_read_count) / mean_sd;
Meine erste Frage ist, ob ich hier korrekt normalisiere und der ZScore korrekt berechnet wird. Es gibt Meinungsäußerungen, dass es hier besser wäre, den Median und die IQR zu verwenden. In einem anderen Versuch wurde zur Normalisierung der Median und für den ZScore die SD verwendet. Die Daten wurde von einer Kollegin mittels Chi-Quadrat Test positiv auf normalverteilt geprüft.
Folgend ein Auszug eines "auffälligen" Datensatzes.
- Code: Alles auswählen
chr_id raw_read_count relative_read_count normalized_relative_read_count zscore
1 124529 0.081706259 1.028254056 0.782
2 130181 0.085414663 1.003419617 0.091
3 98823 0.064839978 0.939839292 -1.140
4 91334 0.059926278 0.994322387 -0.135
5 94004 0.061678125 1.003380166 0.086
6 85497 0.056096492 0.993048896 -0.150
7 81898 0.053735108 1.014597275 0.269
8 75258 0.049378455 0.993395724 -0.119
9 88133 0.057826030 1.408522793 7.562
10 72208 0.047377282 0.976682743 -0.536
11 71385 0.046837293 0.980431495 -0.555
12 69388 0.045527017 0.982323534 -0.369
13 48633 0.031909198 0.958932496 -0.957
14 47351 0.031068049 1.020121829 0.290
15 28902 0.018963248 0.684355230 -4.666
16 43330 0.028429781 0.968643082 -0.511
17 39470 0.025897149 0.968191302 -0.418
18 40037 0.026269170 0.996580865 -0.049
19 32475 0.021307573 1.122080048 1.441
20 37427 0.024556691 0.965946772 -0.520
21 19089 0.012524719 0.989224360 -0.152
22 20312 0.013327157 1.046766760 0.500
23 83856 0.055019795 0.953661548 -1.053
24 217 0.000142379 1.059115389 0.385
25 369 0.000242109 0.764575662 -0.214
Mit der beschriebenen Methodik können wir auffällige von unaufälligen Proben unterscheiden.
Anbei drei Bilder zur Veranschaulichung:
Ein "auffälliger" Fall:
Die Farbe der Punkte korreliert mit dem ZScore. Die Verteilung auf der Y-Achse mit dem normalized_realtive_read_count (d.h. der biologischen Signifikanz).
Eines der "besseren" Kontrollsamples
Eines der "schlechteren" Kontrollsamples
Ich möchte abschließend zur Diskussion stellen, ob es vernünftig / statistisch machbar ist, eine Auswahl der Kontrollen anstatt aller Kontrollen zu verwenden.
Ich kenne es aus anderen Programmen, dass Kontrollen nach ihrer "Ähnlichkeit" zum Sample gereiht und davon nur die besten / ähnlichsten 30 genommen werden. Mindestens allerdings 15 egal wie gut oder schlecht sie passen.
Ich komme hier auf keine Lösung, wie ich ein Kontrollsample mit meinem Testsample in ihrer Gesamtheit auf Ähnlichkeit prüfen kann. Für jedes einzelne Chromosom ja aber insgesamt?
Warum interessiert mich das?
Da sehr viele Kontrollen wunderbar auf der 1er Linie liegen bekomme ich für einige Chromosomen relativ rasch hohe ZScores bei kleinen Abweichungen. Dies würde ich gerne zum Zwecke der Interpretierbarkeit umgehen.
Macht das für Euch Sinn?
Beste Grüße,
Thomas