Diagramy Venna i transkryptomika

tl;dr: Rok temu napisałem pracę, która chodziła za mną od lat, i która moim zdaniem jest bardzo ważna. Nikogo najwyraźniej jednak nie zainteresowała, dlatego teraz Wam o niej opowiem. Będą w notce setki błędnych prac, będzie jeden zwrot akcji, a nawet wystąpi piła łańcuchowa.

Już ładnych parę lat temu ukazała się praca (Niewenhuis et al. 2011), którą od tamtej pory zamęczam wszystkich wokół, czy chcą tego, czy nie. Praca opisuje bardzo powszechny błąd statystyczny – Niewenhuis et al. sprawdzili kilka setek prac z topowych czasopism w dziedzinie neurobiologii. Błąd był w połowie prac, w których można go było popełnić.

Wyobraźcie sobie, że badacie dwie linie komórkowe (albo dwa szczepy myszy, albo…). Jedna linia to zwykłe komórki HeLa, nazwiemy tę linię “WT” (wildtype). W drugiej dokonaliście modyfikacji, wyłączając pewien gen. Tę linię nazwiemy “KO” (knockout).

Pytanie brzmi, czy wyłączenie genu sprawiło, że komórki reagują inaczej na stymulację pewną cytokiną. Oba szczepy, WT i KO, poddajecie więc zabiegowi, i mierzycie jak wiele IL6 produkują komórki. Mamy więc już cztery grupy: kontrolę (K) i stymulację (S) zarowno w WT, jak i KO.

I teraz clou: w WT istnieje statystycznie istotna (p = 0.049) różnica między K a S. Natomiast w KO różnica nie jest statystycznie istotna (p = 0.051). Czy wobec tego uprawnione jest twierdzenie, że rzeczywiście istnieje różnica między WT a KO? Czy one rzeczywiście inaczej reagują na stymulację?

Prawidłowa odpowiedź brzmi: nie, takie stwierdzenie nie jest uprawnione. Jak to ładnie ujął Andrew Gelman, różnica między istotnością statystyczną a jej brakiem sama nie jest statystycznie istotna.

Dlaczego? Przyjrzyjmy się dwóm sytuacjom. Po pierwsze, taki trywialny przypadek:

Jak widać, w tym wypadku rzeczywiście jest tak, jak podejrzewamy: w WT produkcja IL6 idzie w górę, a w KO nawet trochę wędruje w dół. Różnica w KO jest niewielka, i pewnie właśnie dlatego różnica jest istotna w WT, a nie jest istotna w KO. Ale moment, bo możliwa jest też taka sytuacja:

Widać wyraźnie, że obie linie zachowują się tak samo. Co prawda i tutaj jest mała różnica: w WT produkcja IL6 rośnie trochę szybciej niż w KO. Może dlatego różnica w KO jest nieistotna statystycznie. Jednak jeśli policzymy jaka jest różnica między stymulacją a jej brakiem w WT, a potem zrobimy to samo w KO, i sprawdzimy, czy te różnice się różnią¹, to się okaże że nie. Taką “różnicę różnic” nazywamy w statystyce interakcją.

Interakcję da się normalnie przebadać odpowiednią statystyką, ale w tym celu musielibyśmy sformułować i przetestować właściwy model. O tym, czy jest różnica między KO a WT w ich reakcji na stymulację nie dowiemy się po prostu porównując dwie wartości p, bo nie jesteśmy w stanie rozróżnić tych dwóch przypadków powyżej.

Dowcip polega na tym, że te dwie sytuacje:

• “istnieją statystycznie istnotne różnice”, oraz
• “nie ma statystycznie istotnych różnic”

nie są symetryczne!

W rzeczywistości to drugie stwierdzenie powinno bowiem brzmieć inaczej: “nie udało nam się wykryć statystycznie istotnych różnic”. Ale to, że nam się nie udało, nie znaczy że różnic nie ma. “Brak dowodów to nie dowód braku” (absence of proof is not proof of absence).

W naszej sytuacji możemy więc wyprowadzić wniosek, że różnice istnieją w WT, ale nic nie możemy powiedzieć na temat KO. To, że mamy statystycznie istotną różnicę w WT, a nie mamy jej w KO, nie oznacza, że w KO nie ma różnicy - a jedynie, że nie udało nam się tej różnicy wykryć.

Jakie to wszystko ma znaczenie?

Kiedy uczę statystyki, czasami na zajęciach pokazuję studentkom i studentom piłę mechaniczną - jako metaforę. Statystyka daje nam potężne narzędzia, którymi można zrobić wiele ciekawych rzeczy. Niestety, można też się poważnie skaleczyć, jeśli nie zachowujemy podstawowych zasad bezpieczeństwa.

Błędne użycie statystyki prowadzi do błędnych wniosków. Nieuwenhuis et al. twierdzą, że w dwóch trzecich przypadków błędnej analizy, które wykryli, błąd mógł mieć “poważne skutki”, bo sytuacja miała się tak, jak na drugim z obrazków powyżej. Badacze publikują więc obserwacje na temat różnic, które najprawdopodobniej nie istnieją.

Co więcej, jeszcze jedna rzecz może odgrywać rolę. Znacznie łatwiej natrafić na sytuację, w której jedno porównanie daje istotność, a drugie nie – niż pokazać istotną interakcję. Innymi słowy, błędna analiza częściej dostarczy nam wyników na poparcie naszej hipotezy.

Nieuwenhuis et al. przyjrzeli się też pracom z dziedziny biologii komórki i biologii molekularnej. Tu nie udało im się znaleźć ani jednej prawidłowej analizy interakcji. Pokrywa się to, niestety, z moim doświadczeniem.

Za to pracując przez lata w dziedzinie analiz transkryptomowych natrafiłem na pewną bardzo szczególną mutację tego błędu. Szczególną, bo nie dość, że wyjątkowo popularną, to jeszcze dającą pozornie obłędne wyniki!

Dość dawno temu zauważyłem, że w wielu pracach w których pojawia się analiza transkryptomu można zobaczyć takie oto ilustracje:

To, co widzimy powyżej, to diagram Venna. Powyższa ilustracja pochodzi z pracy, w której zbadano dwie grupy pacjentów, G1 oraz G2 (zaraz wyjaśnię, co to za grupy: na razie to tajemnica). Każdemu pacjentowi czy pacjentce pobrano krew i przeprowadzono RNA-Seq. RNA-Seq mierzy względny poziom ekspresji kilkunastu tysięcy naszych genów. Bez wchodzenia w szczegóły, mamy więc dla każdej osoby kilkanaście tysięcy pomiarów, każdy odpowiadający jednemu genowi.

Ekspresja genu: liczba molekuł mRNA wytworzonych przez komórkę dla danego genu. Im większa, tym silniej gen ulega ekspresji, i tym samym tym więcej komórka może syntetyzować kodowanego przezeń białka. Wysoka ekspresja = dużo białka; niska ekspresja = mało białka. Ekspresja genów mówi nam zatem, jak komórka reaguje na pewne bodźce.

W każdej z obu grup byli chorzy na COVID19 oraz kontrole, które nie miały infekcji Sars-Cov-2. Pytanie zadane przez badaczy brzmiało: czy istnieją różnice między G1 i G2 w tym, jak chorzy reagują na Sars-Cov-2? Hipoteza była taka, że grupy G1 oraz G2 różnią się w odpowiedzi immunologicznej, co wychodzi gdy ktoś złapie Sars-Cov-2.

Jak widzicie, jest to typowe pytanie o interakcję: interakcję między infekcją Sars-Cov-2 (lub jej brakiem), a byciem w grupie G1 lub G2. Ale autorzy pracy zamiast poprawnie policzyć interakcję dla każdego genu, postąpili inaczej.

Autorzy osobno porównali infekcję Sars-Cov-2 z kontrolą najpierw w grupie G1, a potem w grupie G2. A potem podzielili geny na cztery grupy:

• geny istotnie zmieniające swoją ekspresję tylko w grupie G1
• geny istotnie zmieniające swoją ekspresję tylko w grupie G2
• geny istotnie zmieniające swoją ekspresję w obu grupach
• cała reszta

Wyniki ilustruje właśnie powyższy diagram Venna. Zgodnie z wnioskowaniem autorów, 431 genów było specyficznych dla grupy G1, i nie pojawiały się w grupie G2. A 278 było specyficznych dla G2, i nie pojawiały się w grupie G1. To jeszcze nic! Te różnice totalnie pasowały do pierwotnej hipotezy. Analiza bogacenia zbioru genów (gene set enrichment analysis) pokazała, że wzbogacone grupy genów to właśnie grupy genów odpowiadające za różne rodzaje odpowiedzi immunologicznej:

Powyżej widać wyniki analizy wzbogacenia dla grupy G1 (po lewej) i dla grupy G2 (po prawej). W G1 mamy wzbogacenie w geny związane np. ze ścieżką TLR4, podczas gdy w G2 widzimy odpowiedź na interleukinę 7. Najwyraźniej w G1 mamy do czynienia z innym rodzajem odpowiedzi immunologicznej - co potwierdza hipotezę autorów.

Analiza wzbogacenia grupy genów: mamy jakąś grupę genów, które nas interesują, np. geny związane z odpowiedzią na infekcję wirusową. Nazwijmy sobie tę grupę “X”. I teraz pytamy: czy pośród genów, których ekspresja ulega zmianie, genów z grupy “X” jest proporcjonalnie więcej niż wśród pozostałych genów?

A teraz opowiemy sobie, co to za grupy i kim byli autorzy.

Autorzy pracy to ja i dwóch moich kolegów (Weiner et al. 2022), a grupy G1 i G2 niczym się nie różniły: losowo wyciągnąłem je z homogennej populacji w której byli pacjenci COVID19 i kontrole. Powtórzę: w rzeczywistości nie było żadnych różnic między G1 a G2.

Zrobiłem to, żeby zademonstrować błędność takiego rozumowania.

Poprawnie przeprowadzona analiza interakcji nie wykazała żadnych różnic. Wszelkie różnice widoczne na obrazku powyżej to totalnie artefakty.

Dowcip polega więc właśnie na błędnej analizie interakcji, dokładnie jak w pracy Nieuwenhuis et al. To, że gen jest statystycznie istotny w grupie G1, a nie jest w grupie G2, to jeszcze nie znaczy, że jego ekspresja nie różni się w grupie G2. Po prostu tego nie wiemy. Ale na obrazku wygląda to znakomicie, nie?

No dobrze, ale czy faktycznie jest to taki częsty błąd?

Napisałem tę pracę właśnie dlatego, że ciągle się z takimi obrazkami i dokładnie takim rozumowaniem stykam. Widzę te diagramy Venna pracach naukowych z górnej półki i na posterach konferencyjnych, i bardzo często są właśnie połączone z błędną interpretacją. Geny osiągające statystyczną istotność w jednej tylko grupie określa się jako “specyficzne” dla tej grupy; a ich liczba ma świadczyć o różnicach między grupami.

Systematyczna analiza 282 artykułów z kilku czasopism pokazała, że jedna trzecia artykułów wymieniających “differential expression” i “venn diagram” popełnia ten sam błąd. Dotyczy to zarówno artykułów ze “Scientific Reports”, jak i “Nature Communications” a nawet “Science Immunology”. Podejrzewam, że moja metoda wyszukiwania problematycznych artykułów nawet zaniża te proporcje.

Czy to jednak diagramy Venna są winne? Same w sobie oczywiście są tylko narzędziem. Ale taka wizualizacja danych aż się prosi o błędną interpretację. Kiedy wysyłaliśmy artykuł, redaktor miał olbrzymi problem z oskarżeniem diagramów Venna, przyrównał to do zabraniania noży kuchennych dlatego że można ich użyć do popełnienia przestępstwa. Odpowiedziałem, że gdyby jedna trzecia przypadków użycia noża kuchennego kończyła się morderstwem, to wypada ostrzec wszystkich, którzy ich używają.

Ale to wciąż nie koniec całej historii; póki co pokazałem tylko występowanie tego samego błędu, o jakim pisali Niewenhuis i ska.

Wspomniałem przed chwilą o analizie bogacenia. To potężna i często stosowana metoda. W telegraficznym skrócie sprawa wygląda tak.

Analiza ekspresji genów dostarcza nam gigantyczne listy genów, których ekspresja ulega zmianie pod wpływem jakiegoś bodźca. Najczęściej takie listy trudno zinterpretować, bo co - będziemy po kolei opisywać tysiąc różnych genów, czy jak? Analiza bogacenia genów (gene set enrichment analysis) zaczyna od uprzednio zdefiniowanych zestawów genów. Na przykład, geny reagujące na interferon (interferon signalling genes, ISG), albo geny związane z pewnym rodzajem komórek. W analizie bogacenia pytamy, jak zachowują się te grupy, a nie pojedyncze geny. Czy wśród genów odpowiadających na pewien bodziec jest więcej genów reagująych na interferon niż oczekiwalibyśmy z losowego rozkładu, a więc czy ta grupa genów jest nadreprezentowana w wynikach. Jeśli tak, to możemy wyciągnąć wniosek, że mamy do czynienia z aktywacją ścieżki interferonowej.

Taką analizę błędne prace często łączą z diagramami Venna. Po prostu pyta się, czy wśród tych 431 genów istotnych są nadreprezentowane geny związane z pewnymi ścieżkami. W świetle tego, co Wam opowiedziałem przed chwilą, jest to oczywiście nieprawidłowe. Ale, niestety, bardzo przekonujące.

Problem w tym, że te 431 genów to nie są przypadkowe geny. To są geny, które faktycznie zmieniają ekspresję w grupie G1; a ponieważ między G1 a G2 nie ma istotnych różnic, to one prawie na pewno również zmieniają swoją ekspresję w grupie G2, tylko że nie udało nam się tej zmiany wykryć. Innymi słowy, to są geny powiązane ogólnie z odpowiedzią na Sars-Cov-2, a więc wzbogacone w grupy genów odpowiedzi immunologicznej na zakażenie wirusowe.

Nie tylko więc te 431 genów jest artefaktem, ale jest, co gorsza, artefaktem bardzo atrakcyjnym: bo autentycznie wskazującym na różne grupy genów związanych z odpowiedzią immunologiczną, a więc pasującym do naszej hipotezy!

Myślę, że właśnie dlatego te nieszczęsne diagramy Venna (i związane z nimi błędne rozumowanie) są tak rozpowszechnione. Wyobraźcie sobie eksperymentatorów przeprowadzających test kliniczny, testujących swoją hipotezę tak, jak to opisaliśmy. Pięknie wychodzą różnice, można wokół nich napisać całą historię, jak to w G1 mamy do czynienia z TLR4, a za to w G2 mamy aktywację IL7. A potem przychodzi bioinformatyczka i mówi im, że w ich danych nic nie ma, grupy się niczym nie różnią, żadne wzbogacenia nie wychodzą. Całe badanie psu w odbyt, a zamiast artykułu w Nature publikacja w Sci Rep. W zależności od pozycji bioinformatyczki, albo eksperymentatorzy nie ulegną, albo ulegną, ale następnym razem nie będą jej pytać o zdanie.

To nie jest wymyślona sytuacja. To mi się regularnie zdarza. I dlatego napisałem tę pracę.

• Nieuwenhuis S, Forstmann BU, Wagenmakers EJ. Erroneous analyses of interactions in neuroscience: a problem of significance. Nature Neuroscience. 2011 Sep;14(9):1105-7.
• Gelman A, Stern H. The difference between “significant” and “not significant” is not itself statistically significant. The American Statistician. 2006 Nov 1;60(4):328-31.
• Weiner 3rd J, Obermayer B, Beule D. Venn diagrams may indicate erroneous statistical reasoning in transcriptomics. Frontiers in Genetics. 2022 Apr 14;13:818683.