W niniejszym badaniu dokonano ponownej analizy zanonimizowanych, ogólnodostępnych danych, co nie wiązało się z nowym naborem uczestników, interwencją ani pobieraniem próbek. W związku z powyższym, dla niniejszej analizy wtórnej nie było wymagane uzyskanie zgody komisji etyki instytucjonalnej ani nowej świadomej zgody uczestników. Zgoda komisji etycznej oraz świadoma zgoda dla badań pierwotnych zostały zgłoszone przez odpowiednich autorów danych.
Projekt badania i logika przepływu pracy
Przeprowadziliśmy retrospektywne bioinformatyczne badanie oparte na publicznych danych, wykorzystując zbiory danych zdeponowane w bazie Gene Expression Omnibus20,21. Nie generowano nowych próbek pacjentów, eksperymentów na hodowlach komórkowych, modeli zwierzęcych ani walidacji laboratoryjnej (wet-lab). Przepływ pracy oparto na projekcie z priorytetem dla kotwic (anchor-first design). Zbiór GSE117827 wykorzystano do szacowania wielkości efektu, konstrukcji modułów, oceny zgodności między przedziałami oraz analizy gradientu objawów. Niezależne zestawy danych wprowadzono dopiero po ustaleniu przynależności do modułów. Przepływ pracy obejmował kotwiczenie w parach przedziałów, mapowanie HGNC i filtrowanie białek kodujących, oddzielną konstrukcję modułów dla nosa i krwi, walidację z dopasowaniem tkankowym i podłużną, benchmarking kliniczny oraz analizy odpornościowe. Analizy potwierdzające obejmowały zablokowane testy z dopasowaniem tkankowym i testy podłużne. Analizy genów nakładających się, gradientu objawów, programów markerowych oraz analizy wariancji miały charakter eksploracyjny lub opisowy.
Kohorta zakotwiczenia i główny kontrast
Zbiór GSE117827 zawiera profile ekspresji z wymazów z nosa i krwi pełnej od 26 dzieci: 9 przypadków objawowego zakażenia pikornawirusami, 5 przypadków bezobjawowego wykrycia pikornawirusów, 6 przypadków objawowego zakażenia RSV oraz 6 bezobjawowych kontroli u których nie wykryto wirusów17. W dwóch przypadkach RSV brakowało próbek krwi. Zatem pełny projekt zakotwiczenia obejmował 50 próbek, a główny kontrast między grupą zakażoną a kontrolną zawierał 40 próbek po wykluczeniu bezobjawowych detekcji pikornawirusów z konstrukcji modułów. W głównym kontraście próbki z objawowym zakażeniem pikornawirusami oraz objawowym zakażeniem RSV oznaczono jako zakażone, natomiast bezobjawowe próbki u których nie wykryto wirusów oznaczono jako kontrole. Bezobjawowych detekcji pikornawirusów nie traktowano jako kontroli, ponieważ wykrycie wirusa bez objawów jest biologicznie odmienne od stanu zdrowia bez obecności wirusa. Próbki te pominięto podczas konstrukcji modułów, a następnie wykorzystano do analizy gradientu objawów.
Mapowanie sond i przetwarzanie wstępne
Zestaw GSE117827 został sprofilowany na platformie GPL23126. Identyfikatory klastrów transkrypcyjnych zmapowano na symbole genów, korzystając z pliku adnotacji Clariom D Human na36, hg38 udostępnionego przez platformę GEO GPL24539. Zachowano wyłącznie symbole zatwierdzone przez HGNC22. Usunięto sondy kontrolne, sondy ERCC, niezmapowane klastry transkrypcyjne oraz symbole niezatwierdzone. W przypadku wielu klastrów transkrypcyjnych mapujących na ten sam zatwierdzony symbol, dane scalono, obliczając medianę ekspresji. Po filtrowaniu i scalaniu do analizy kotwic dostępnych było 27 685 genów. Transformację log₂(x + 1) zastosowano tylko wtedy, gdy 95. percentyl macierzy ekspresji przekraczał 50, co wskazywało na niezalogarytmizowaną skalę intensywności. Geny z więcej niż 20% brakującymi wartościami zostały usunięte; pozostałe brakujące wartości zastąpiono medianą w obrębie danego genu. Następnie każdy gen został ustandaryzowany we wszystkich próbkach w danym zbiorze danych według wzoru z = (x − średnia)/odchylenie standardowe próbki (ddof = 1). Genom o zerowej wariancji przypisano ustandaryzowaną wartość 0. Konstrukcja modułów została ograniczona do wpisów sklasyfikowanych jako geny kodujące białka w kompletnym zbiorze HGNC.
Analiza wielkości efektu i konstrukcja modułów
Przedziały nosowe i krwi przeanalizowano oddzielnie. Próbki wirusowe objawowe porównano z kontrolami wirusowo ujemnymi, stosując standaryzowane różnice średnich Hedges g oraz dwustronne testy Welcha23. Wartość Hedges g obliczono jako różnicę średnich (grupa zakażona minus grupa kontrolna) podzieloną przez połączone odchylenie standardowe i pomnożoną przez poprawkę dla małych próbek 1 − 3/(4df − 1), gdzie df = nzakażone + nkontrolne − 2. Testy Welcha przeprowadzono przy założeniu nierówności wariancji. Współczynniki fałszywych odkryć Benjamini–Hochberga obliczono dla wszystkich badanych genów w obrębie każdego przedziału24. Ponieważ w małej kohorcie kotwic żadne geny kodujące białka nie spełniły wcześniej określonego rygorystycznego kryterium kombinowanego (FDR < 0,05 oraz Hedges g > 0,8), geny o dodatnim efekcie uszeregowano najpierw według rosnącej dwustronnej wartości P testu Welcha, a następnie według malejącej wartości Hedges g, stosując symbol genu jako ostateczny deterministyczny czynnik rozstrzygający przy identycznych wartościach. 50 najwyżej sklasyfikowanych genów utworzyło każdy główny moduł. Liczbę pięćdziesięciu genów wybrano a priori jako umiarkowaną wielkość modułu, która zachowuje szerokość biologiczną przy jednoczesnym ograniczeniu utraty danych z powodu brakujących genów pomiędzy platformami. Wielkość modułu nie była dostrajana do zewnętrznych wyników AUROC. Analizy czułości z wykorzystaniem 10, 25, 50, 100 i 200 genów przyniosły ten sam jakościowy wynik separacji kotwic. Celem była odtwarzalność na poziomie modułu, a nie identyfikacja pojedynczych genów.
Zgodność międzyprzedziałowa
Szacunki Hedgesa g dla nosa i krwi dopasowano według genów i porównano za pomocą korelacji Pearsona i Spearmana. Nakładanie się 50 najwyższych modułów dla nosa i 50 najwyższych dla krwi podsumowano za pomocą liczby wspólnych genów oraz indeksu Jaccarda. Geny nakładające się zinterpretowano jako eksploracyjne kandydaty do nakładania się rang międzyprzedziałowych, a nie jako zwalidowane, konserwatywne biomarkery.
Zewnętrzna walidacja z dopasowaniem do tkanki
Walidacja wyników modułów została przeprowadzona z wykorzystaniem niezależnych publicznych zbiorów danych z interpretowalnymi grupami źródłowymi. Walidacja dla górnych dróg oddechowych obejmowała próbki z płukania nosa RSV z zestawu GSE41374 oraz zbiory danych SARS-CoV-2 GSE152075 i GSE156063. Walidacja krwi obejmowała zbiór GSE171110 oraz dwie oddzielnie znormalizowane jednostki GSE38900: GPL10558 (36 próbek; 28 RSV i 8 kontrolnych) oraz GPL6884 (138 próbek; 107 RSV i 31 kontrolnych). Dla każdego zewnętrznego zbioru danych sondy zostały przypisane do symboli HGNC, a zdublowane symbole zostały scalone poprzez obliczenie średniej ekspresji. W obrębie każdego zbioru danych zastosowano te same reguły transformacji, filtrowania brakujących wartości, imputacji medianą oraz standaryzacji z> dla poszczególnych genów, które wykorzystano dla zbioru kotwiczącego. Wynik modułu obliczono jako nieważoną średnią zestandaryzowanych wartości ekspresji dla genów wchodzących w skład danego modułu. Wartości AUROC, średnia precyzja (average precision) oraz FDR z testu Welcha zostały podane jako mierniki przenośności, a nie jako szacunki klinicznej wydajności diagnostycznej.
Duże kliniczne zestawy referencyjne i zestawy danych z perturbacji ex vivo
Dwa zestawy danych z krwi pełnej zostały wykorzystane jako punkty odniesienia (benchmarks) kliniczne, a nie jako kohorty odkrywcze. W GSE63990 próbki przypisano na podstawie zdeponowanych metadanych stanu zakażenia do grupy wirusowych infekcji dróg oddechowych (n = 117), bakteryjnych infekcji dróg oddechowych (n = 73) lub chorób niezakaźnych (n = 90); próbki bez jednej z tych jednoznacznych etykiet zostały wykluczone18. W GSE40012 zdeponowane pola diagnostyczne zmapowano na zapalenie płuc wywołane wirusem grypy A (n = 39), bakteryjne zapalenie płuc bez grypy (n = 61), Zespół Odpowiedzi Systemowej na Zapalenie (SIRS) bez zapalenia płuc (n = 40), zdrowe kontrole (n = 36) lub mieszane bakteryjne/grypowe zapalenie płuc (n = 14)19. Próbki z mieszanym zapaleniem płuc bakteryjnym/grypowym zostały opisane, lecz wykluczono je z binarnych porównań referencyjnych.
Zbiór GSE53543 zawiera 196 profilów PBMC ex vivo pochodzących od 98 osób. Każda osoba dostarczyła jedną próbkę kontrolną (tylko medium) oraz jedną próbkę eksponowaną na rhinovirus 16 przez 24 h; identyfikatory badanych potwierdziły istnienie 98 pełnych par. Główny punkt odniesienia raportował AUROC dla 98 próbek stymulowanych i 98 próbek niestymulowanych, ponieważ AUROC jest miarą separacji rang. Parowanie pacjentów zachowano w metadanych i wykorzystano w sparowanej analizie czułości różnic w wynikach. Ten eksperyment był analizowany oddzielnie od naturalnych infekcji klinicznych i nie był wykorzystywany do wspierania roszczeń diagnostyki klinicznej.
Pobrano surowe macierze serii GEO dla GSE63990, GSE40012 oraz GSE53543. Dla GSE63990 zastosowano GPL571, dla GSE40012 GPL6947, a dla GSE53543 GPL10558. Sondy przypisano do symboli zatwierdzonych przez HGNC, a powtarzające się przypisania zredukowano poprzez obliczenie średniej ekspresji. Wykorzystano zdeponowane macierze znormalizowane. Ogólna zasada 95. percentyla spowodowała zastosowanie transformacji log₂(x + 1) wyłącznie dla macierzy w skali nielogarytmicznej. Zbiór GSE53543 został już poddany transformacji log₂, normalizacji niezależnej od rangi oraz korekcie ze względu na dzień przetwarzania przez pierwotnych badaczy, w związku z czym nie stosowano dodatkowej transformacji logarytmicznej. Po usunięciu genów z więcej niż 20% brakującymi wartościami i uzupełnieniu pozostałych braków wartością mediany, każdy gen został poddany standaryzacji z we wszystkich próbkach w obrębie danego zbioru danych. Warstwa referencyjna zawierała 470 klinicznych próbek krwi pełnej oraz 196 próbek PBMC ex vivo.
Walidacja sygnatury referencyjnej
Kotwiczne moduły nosowe i krwi zostały poddane walidacji w odniesieniu do trzech nieważonych zestawów referencyjnych. Oficjalny zestaw 33 mRNA według Pandyi został opracowany na podstawie uzupełniającej tabeli 1 z oryginalnego raportu klasyfikatora5. Komparator Andres-Terre obejmował 33 geny ukierunkowane na interferon, wyselekcjonowane z raportowanej sygnatury wielowirusowej3. Komparator Hallmark obejmował rdzenny podzbiór 33 genów interferonu alfa powiązany ze zbiorem MSigDB HALLMARK_INTERFERON_ALPHA_RESPONSE25,26. Pełne listy genów dla wszystkich trzech zestawów referencyjnych znajdują się w Supplementary Data 1. Wyniki referencyjne te nie odtwarzają oryginalnych klasyfikatorów ważonych. Dla każdego zbioru danych obliczono wynik jako nieważoną średnią dostępnych wyników z dla poszczególnych genów; wymagano obecności co najmniej trzech reprezentatywnych genów, a liczba reprezentowanych genów została odnotowana. Kontrasty walidacyjne obejmowały: wirusowe vs bakteryjne, wirusowe vs niezakaźne, wirusowe vs bakteryjne-lub-niezakaźne, wirusowe vs zdrowe/kontrolne oraz bakteryjne vs zdrowe/kontrolne, w zależności od dostępności wymaganych grup. Jako mierniki walidacyjne zinterpretowano AUROC oraz średnią precyzję. Wartości P z testu Welcha zostały skorygowane za pomocą procedury Benjamini–Hochberg dla wszystkich poprawnych kombinacji zbioru danych, kontrastu i modułu w tabeli walidacyjnej.
Niezależna walidacja podłużna
Sparowane zestawy danych z fazy ostrej i wypisu GSE97741 oraz GSE97742 zostały wykorzystane wyłącznie do walidacji podłużnej27. Próbki oznaczone w zdeponowanych metadanych jako pojedyncza infekcja RSV (RSVsi) lub rinowirus (hRV) stanowiły podstawę analizy; koinfekcje RSV (RSVco) zachowano jedynie w wynikach uzupełniających. Etykiety fazy ostrej i wypisu wyekstrahowano z nazw próbek. Próbki dopasowano według zestawu danych, przedziału, grupy wirusów i identyfikatora pacjenta, zachowując jedynie pacjentów z oboma punktami czasowymi. Główna grupa połączona składała się z 38 par RSVsi i 30 par hRV (68 par na przedział).
Do selekcji genów, doprecyzowania modułów ani dostrajania progów nie wykorzystano zbiorów GSE97741 ani GSE97742. Zbiory te zarezerwowano do walidacji zewnętrznej. Identyfikatory sond przypisano do zatwierdzonych symboli HGNC, a duplikaty symboli scalono poprzez obliczenie mediany ekspresji. Przed obliczeniem stałych wyników modułowych dla GSE117827 zastosowano tę samą zasadę transformacji logarytmicznej dla 95. percentyla, filtr wartości brakujących, imputację medianą oraz procedury standaryzacji z dla każdego genu w obrębie zbioru danych.
Dla każdej tkanki i modułu wyniki w fazie ostrej i w momencie wypisu zostały sparowane według obiektu i grupy wirusowej. Dla każdej pary obliczono różnicę między wynikiem w fazie ostrej a wynikiem w momencie wypisu. Wartość Cohen dz obliczono jako średnią różnicę sparowaną podzieloną przez odchylenie standardowe próby. Zgłoszono wyniki dwustronnych sparowanych testów t oraz dwustronnych testów znaków Wilcoxona, wraz z AUROC dla separacji wyników w fazie ostrej względem momentu wypisu. Wartości FDR według metody Benjamini–Hochberg obliczono dla wszystkich 20 prawidłowych sparowanych testów t w pełnych wynikach podłużnych (dwa zestawy danych, dwa źródła modułów i pięć z góry określonych podsumowań grup wirusowych).
Analiza gradientu objawów i analiza post-hoc
Po ustaleniu przynależności do modułów, wykrycia picornavirusu u osobników bezobjawowych w zbiorze GSE117827, które zostały wyłączone z analizy, wykorzystano wyłącznie do analizy gradientu objawów. Porządkowa ocena kliniczna wynosiła 0 dla kontroli ujemnych pod kątem wirusa, 1 dla bezobjawowych wykryć picornavirusu oraz 2 dla infekcji objawowej. Korelację Spearmana zastosowano do oceny trendu porządkowego, natomiast test Kruskala–Wallisa posłużył do oceny ogólnych różnic pomiędzy trzema grupami. Dwustronne testy post-hoc Mann–Whitneya U porównywały trzy pary grup. Korekcję Benjamini–Hochberga zastosowano oddzielnie w każdym przedziale dla rodziny trzech porównań parzystych.
Analiza wzbogacenia funkcjonalnego
Geny z modułów nosowych i krwi przeanalizowano za pomocą narzędzia Enrichr poprzez gseapy, korzystając z bibliotek MSigDB Hallmark 2020, Reactome 2022 oraz GO Biological Process 202325,26,28,29,30,31. Enrichr wykorzystał standardowy model nadreprezentacji oparty na dokładnym teście Fishera. Za istotne uznano zgłaszane przez biblioteki wartości P skorygowane metodą Benjamini–Hochberg < 0.05. Osiem najistotniejszych terminów dla każdego modułu z trzech analizowanych bibliotek uszeregowano według skorygowanej wartości P. Wyniki wzbogacenia wykorzystano wyłącznie do interpretacji.
Analizy odporności, programu markerów i wariancji
Analizy odporności miały na celu sprawdzenie, czy wyniki zależały od wielkości modułu lub przypadkowego doboru. Oceńiono moduły o wielkości 10, 25, 50, 100 i 200 genów. W analizie zerowej dla genów losowych uniwersum kwalifikowalnych genów obejmowało geny kodujące białka HGNC reprezentowane w przetworzonej macierzy GSE117827. Pobrano pięćset 50-genowych zestawów bez zwracania, używając ziarna losowości NumPy 20260622. Ponowny dobór metodą bootstrap ograniczono do 1000 genów kodujących białka o dodatnim efekcie, które zajęły najwyższe miejsca w oryginalnej analizie kotwic dla każdego przedziału. W każdej próbie geny kodujące białka o dodatnim efekcie uszeregowano według rosnącej dwustronnej wartości P testu Welcha, a następnie według malejącej różnicy średnich. Wybrano 50 najwyżej sklasyfikowanych genów. Stabilność genów obliczono jako liczbę wyborów podzieloną przez 100. Przedstawiono 20 genów o najwyższej częstotliwości wyboru w każdym przedziale.
Wyniki programów markerowych wykorzystano jako pomoce opisowe, a nie jako szacunki frakcji komórkowych. Przeanalizowano sześć wyselekcjonowanych programów: nabłonkowy (EPCAM, KRT8, KRT18, KRT19, MUC1, KRT5, KRT14, SCGB1A1, FOXJ1), monocytarny/makrofagiczny (LYZ, LST1, S100A8, S100A9, FCGR3A, FCGR1A, CD14, MS4A7, CTSS), neutrofilowy (S100A8, S100A9, MPO, ELANE, CEACAM8, FCGR3B, OLFM4, MMP8), T/NK (CD3D, CD3E, TRAC, NKG7, GNLY, PRF1, GZMB, KLRD1, IL7R), B/plazmatyczny (MS4A1, CD79A, CD79B, MZB1, JCHAIN, IGHG1, SDC1) oraz mieloidalny interferonowy (SIGLEC1, IFI27, IFI44L, ISG15, MX1, OAS1, RSAD2, IFIT3). Wynik każdego programu obliczono jako średnią dostępnych wartości z dla poszczególnych genów, przy wymaganiu obecności co najmniej trzech reprezentatywnych genów. Wartości P testu Spearmana skorygowano za pomocą procedury Benjamini–Hochberg dla pełnej rodziny korelacji zbiór danych–moduł–program. Jednoczynnikowe eta-kwadrat obliczono jako sumę kwadratów międzygrupowych podzieloną przez całkowitą sumę kwadratów dla każdej pary moduł–czynnik. Wiersze, w których brakowało odpowiedniego czynnika, wykluczono z tych obliczeń. Analiza ta miała charakter opisowy i nie uwzględniała korekty dla czynników wzajemnie skorelowanych.
Powtarzalność
Wszystkie analizy przeprowadzono przy użyciu skryptowych przepływów pracy w języku Python 3.12.13, z wykorzystaniem pakietów oprogramowania i wersji podanych w Tabeli materiałów. W procedurach randomizowanych zastosowano stałe ziarno (seed) o wartości 20260622. Publiczne macierze ekspresji z bazy GEO przetwarzano przy użyciu spójnej struktury projektu, która rozdzielała surowe dane wejściowe, dane przetworzone, wyniki statystyczne, tabele i ryciny. Definicje modułów, rekordy kurateli próbek, szacunki wielkości efektu, statystyki walidacyjne, wyniki wzbogacenia, analizy odporności oraz dane źródłowe rycin zostały zachowane w celu umożliwienia niezależnej weryfikacji. Kod analityczny, specyfikacje zależności oraz wspierające pochodne dane są dostępne zgodnie z opisem w oświadczeniu o dostępności danych.