od jednoho matematika z Hamburku, který by prozatím rád zůstal v anonymitě

Shrnutí

De novo metatranskriptomické sekvenování nebo kompletní sekvenování genomu jsou ve virologii uznávané metody detekce deklarovaných patogenních virů. Při tomto postupu nejsou detekovány žádné virové částice (viriony), a nejsou ve smyslu slova izolace ani izolovány, ani biochemicky charakterizovány. V případě SARS-CoV-2 se celková RNA často extrahuje ze vzorků pacientů (např. z bronchoalveolární laváže (BAL) nebo z výtěrů z krku a nosu) a následně se sekvenuje. Neexistuje žádný důkaz toho, že fragmenty RNA použité k výpočtu sekvencí virového genomu jsou virového původu.

Z toho důvodu jsme zkoumali publikaci „A new coronavirus associated with human respiratory disease in China“ [1] a s ní i příslušná sekvenační data biologického projektu s číslem ID PRJNA603194 z 27.01.2020 publikovaná k původnímu návrhu sekvence genu pro SARS-CoV-2 (GenBank: MN908947.3). Opakované sestavení de novo pomocí programu Megahit (verze 1.2.9) ukázalo, že publikované výsledky nelze reprodukovat. V rozporu s tím, co bylo uvedeno ve studii výše [1], jsme pravděpodobně detekovali (ribozomální) ribonukleové kyseliny lidského původu. Další analýzy poskytly důkazy o možné nespecifické amplifikaci načtených sekvencí během PCR konfirmace a určování genomových sekvencí, které nejsou asociovány se SARS-CoV-2 (MN908947.3).

Nakonec jsme provedli několik sestavení za použití jiných referenčních genových sekvencí, jako jsou SARS-CoV, virus HIV, hepatitida delta, spalničky, Zika, Ebola nebo Marburg a prozkoumali strukturní podobnost dostupných sekvenčních dat s příslušnými sekvencemi. Získali jsme první indicie, že lze některé sekvence virového genomu, které jsme v této práci sledovali, získat z RNA nesuspektních (neinfikovaných) lidských vzorků.

Klíčová slova

SARS-CoV-2, COVID-19, virus, de novo sestavování, celogenomové sekvenování, bioinformatika, PCR, SARS-CoV, netopýří SARS-CoV, HIV, virus hepatitidy delta, virus spalniček, virus Zika, virus Ebola, virus Marburg.

Úvod

Pro konstrukci sekvencí virového genomu jsou z různých zdrojů nukleových kyselin, jako např. z tekutin z bronchoalveolární laváže (BAL) [1, 2], z výtěrů z nosohltanu [3, 4, 5, 6, 12, 13], z komponentů buněčných kultur nebo supernatantů buněčných kultur [2, 11, 12, 13, 14, 16], stejně jako z lidských [8, 9, 10, 16] a zvířecích vzorků [7, 15] izolovány a sekvenovány nukleové kyseliny (RNA nebo DNA). Přitom nepocházejí získané nukleové kyseliny výhradně z již izolovaných (virových) částic, tj. oddělených od všeho ostatního, ale často z celého vzorku. To znamená, že původ fragmentů nukleových kyselin použitých k výpočtu sekvencí genomu je a priori nejasný.

V případě ribonukleové kyseliny (RNA) je tato nejprve pomocí RNA-dependentní DNA polymerázy přepsána do cDNA. DNA nebo cDNA je poté fragmentována pomocí enzymů a amplifikována pomocí polymerázové řetězové reakce (PCR), a sice před vlastním sekvenováním, tj. určením nukleotidové sekvence krátkých fragmentů DNA nebo cDNA. Při amplifikaci se kromě náhodných sekvencí primerů (náhodných hexamerů) používají v závislosti na referenčním nebo cílovém genomu také vysoce specifické sekvence primerů [např.: 1, 3, 4, 5, 6, 7, 8, 17, 18]. Takto získaná data ze sekvenace se nakonec zpracují pomocí bioinformatických algoritmů.

Pro stanovení sekvencí virového genomu de novo se běžně používají dvě metody: metatranskriptomické sestavování [1, 12] a kompletní sekvenování genomu [3, 4, 5, 6, 17, 18]. Zatímco při de novo metatranskriptomickém sestavování se většinou nepoužívají žádné referenční sekvence, nebo se používají pouze navazující referenční sekvence, při kompletním sekvenování genomu se používá velké množství specifických primerových sekvencí, z nichž některé dohromady samy o sobě pokrývají 4 až 17 % cílového genomu [1, 17]. Pro amplifikaci cDNA se často používá 35 až 45 cyklů [1, 6, 17].

V případě viru SARS-CoV-2 (GenBank: MN908947.3) [1] byla navržená sekvence virového genomu vypočtena pomocí de novo metatranskriptomického sestavení RNA z BAL pacienta z čínského Wuchanu. K sestavení kontigů (rekonstruovaných spojitých úseků) byly použity assemblery Megahit (verze 1.1.3) a Trinity (verze 2.5.1). Megahit vygeneroval celkem 384 096 (200 nt – 30 474 nt) a Trinity 1 329 960 (201 nt – 11 760 nt) kontigů. Pozoruhodné jsou velké rozdíly mezi oběma soubory. Podle [1] vykazoval nejdelší kontig, sestavený pomocí programu Megahit, vysoký překryv (89.1 %) s genomem netopýra SL-CoVZC45 (genová banka: MG772933) a byl použit k designu primerů pro PCR konfirmaci a určení genomových sekvencí.

Organizace virového genomu byla přizpůsobena seřazením sekvencí podle dvou reprezentativníh druhů rodu betacoronavirus, a sice koronavirus asociovaný člověku (SARS-CoV Tor 2, genová banka: AY274119) a koronavirus asociovaný netopýrům (netopýr SL-CoVZC45, genová banka: MG772933).

Ze vzorku pacienta přitom nebyla identifikována a biochemicky charakterizována žádná patogenní virová částice, která by byla jednoznačně spojena se sekvencí MN908947.3. Spíše došlo k tomu, že z BAL pacienta byla extrahována a zpracována veškerá RNA. Chybí tedy důkaz, že k vytvoření deklarovaného virového genomu SARS-CoV-2 byly použity výhradně nukleové kyseliny virového původu. Navíc nebyly v souvislosti s konstrukcí deklarovaného virového genomu zveřejněny žádné výsledky případných kontrolních experimentů. To platí i pro všechny ostatní referenční sekvence, kterým se tato práce věnuje. V případě viru SARS-CoV-2 by byla nasnadě kontrola spočívající ve skutečnosti, že deklarovaný virový genom nelze sestavit z neinfikovaných zdrojů lidské RNA nebo RNA jiného původu.

V této publikaci jsme zkoumali reprodukovatelnost de novo sestavení s využitím původních sekvenačních dat publikovaných k původní práci o koronaviru SARS-CoV-2 [1]. Dále jsme zkoumali strukturní podobnost sekvenačních dat, které jsou aktuálně k dispozici, s dalšími veřejně dostupnými referenčními virovými sekvencemi pro (netopýří) SARS-CoV [1, 7, 13, 14], virus lidské imunodeficience [8], hepatitidu delta [9], spalničky [11, 12], Ziku [10], Ebolu [15], Marburg [16] (tabulky a obrázky: tabulka 3). Za tímto účelem zde představujeme jednoduchý bioinformatický protokol. Při ověření našich výsledků jsme vzali v úvahu i náhodně vygenerované a fiktivní genomové sekvence, abychom vyloučili čirou náhodnost našich výsledků.

Hlavní část

Opakované de novo sestavení již publikovaných sekvenačních dat

Pro opakované sestavení de novo jsme 30. 11. 2021 za pomocí nástrojů SRA stáhli z internetu původní sekvenační data (SRR10971381) pocházející ze dne 27. 1. 2020 [19]. K přípravě párových sekvenačních čtení pro vlastní sestavování pomocí softwaru Megahit (verze 1.2.9) [20] jsme použili preprocesor FASTQ fastp (verze 0.23.1) [21]. Po filtraci párových čtení zůstalo 26 108 482 z původního celkového počtu 56 565 928 čtení o délce přibližně 150 bp. Velkou část sekvencí, pravděpodobně většinou lidského původu, autoři přepsali písmenem „N“ jako neznámá, a proto je fastp odfiltroval. To je třeba ve smyslu vědeckosti považovat za problematické, protože ne všechny kroky lze dohledat nebo reprodukovat. Pro komplexní generování kontigů ze zbývajících krátkých sekvencí jsme použili program Megahit (verze 1.2.9) ve standardním nastavení.

S počtem 28 459 (200 nt – 29 802 nt) kontigů jsme získali výrazně méně, než je popsáno v publikaci [1]. Jinak, než je zobrazeno v publikaci [1], obsahoval nejdelší námi sestavený kontig pouze 29 802 nt, a byl tedy o 672 nt kratší než nejdelší kontig s 30 474 nt, který podle [1] zahrnoval téměř celý virový genom. Náš nejdelší kontig vykazoval dokonalou shodu se sekvencí MN908947.3 o délce 29 801 nt (tabulky a obrázky: tabulky 1, 2). V rámci vědeckého ověření jsme tedy nejdelší kontig o délce 30 474 nt nedokázali reprodukovat. Z toho plyne, že u zveřejněných sekvenčních dat se nemůže jednat o původní čtení, použitá pro sestavování.

Po sestavení kontigů jsme mapováním krátkých sekvencí na 28 459 určených kontigů pomocí Bowtie2 (v.2.4.4) stanovili příslušné pokrytí [22]. Následně jsme 5.12.2021, resp. 20.12.2021, porovnali 50 kontigů s nejvyšším pokrytím a 50 nejdelších kontigů s nukleotidovou databází (Blastn). Podrobné výsledky vyhledávání naleznete v tabulkách a obrázcích: Tabulky 1, 2.

Srovnání našich výsledků (tabulky a obrázky: tabulka 1) s výsledky z [1, doplňující tabulka 1] vykazuje pozoruhodné rozdíly. V následujícím textu je před všechna Contig-Id z [1] uvedena „1_“, abychom tyto lépe odlišili od našich Contig-Id. Obecně lze konstatovat, že výsledky našeho vyhledávání týkající se přístupových kódů neodpovídají přesně výsledkům (neshodují se s nimi) uvedeným v [1]. Co se týká popisů subjektů, zde jsme zaznamenali většinou shodu. Kromě toho, s výjimkou nejdelšího kontigu (1_k141_275316), vykazovaly naše kontigy větší délku a taky větší bohatost pokrytí. Jednoznačný je případ srovnání kontigu 1_k141_179411 s kontigem k141_12253. První z nich má délku 2 733 nt, zatímco druhý 5 414 nt. To je první možný náznak toho, že při PCR konfirmaci s primery zkonstruovanými pro MN908947.3 z 1_k141_275316 došlo k nespecifické amplifikaci sekvenovaných čtení, která nejsou asociována se SARS-CoV-2.

Na tomto místě je třeba podrobněji probrat kontig s identifikací k141_27232, s nímž je spojeno 1 407 705 sekvencí, a tím asi 5 % ze zbývajících 26 108 482 sekvencí. Porovnání s nukleotidovou databází dne 5.12.2021 ukázalo vysokou shodu (98,85 %) s „Homo sapiens RNA, 45S pre-ribosomal N4 (RNA45SN4), ribosomal RNA“ (GenBank: NR_146117.1, ze dne 4.7.2020). Toto pozorování je v rozporu s tvrzením uvedeným v [1], že byla provedena deplece ribozomální RNA a že čtení lidských sekvencí byla filtrována pomocí lidského referenčního genomu (lidská verze 32, GRCh38.p13). Za zmínku stojí skutečnost, že sekvence NR_146117.1 byla publikována až po zveřejnění zde posuzované knihovny sekvencí SRR10971381.

Toto pozorování zdůrazňuje, že je obtížné určit a priori přesný původ jednotlivých fragmentů nukleových kyselin použitých ke konstrukci deklarovaných sekvencí virového genomu.

Analýza sekvenční struktury na základě referenčních zdrojů

V podstatě jsme mapovali párové čtení (2×151 bp) pomocí BBMap [23] k referenčním sekvencím, které zde hodnotíme (tabulky a obrázky: Tabulka 3), s použitím relativně nespecifických nastavení. Poté jsme pomocí reformat.sh měnili minimální délku (M1) a minimální (nukleotidovou) identitu (M2), abychom získali odpovídající podmnožiny již dříve zmapovaných sekvencí s odpovídající kvalitou. Zvýšení minimální délky M1 nebo minimální identity nukleotidů M2 zvyšuje (statistickou – pozn. překlad.) významnost příslušného mapování. Následně jsme s příslušnými podmnožinami vybrané kvality vytvořili konsenzuální sekvence s ohledem na příslušnou uvažovanou referenci. Všechny báze s kvalitou nižší než 20 jsme nastavili na „N“ (neznámé). Kvalita 20 znamená chybovost 1 % na nukleotid, což lze pro účely našich analýz považovat za dostatečné. Nakonec byla posouzena shoda referenčních a konsenzuálních sekvencí pomocí nástrojů BWA [24], Samtools [25] a Tablet [26]. Uspořádaná dvojice (M1; M2) = (37; 0,6) byla zvolena tak, aby byla chybovost F1 a F2 pro referenční LC312715.1 menší než 10 %. Výsledky všech výpočtů jsou uvedeny v tabulkách a obrázcích: tabulka 4.

Výpočty ukazují nejvyšší významnost pro volbu uspořádané dvojice (37; 0,6), což je patrné z hodnot nejvyšší chybovosti. Srovnatelnou významnost mají uspořádané dvojice (47; 0,50) a (25; 0,62). Zatímco sekvence genomu spojované s koronaviry vykazují chybovost vyšší než 10 % u všech uvažovaných uspořádaných dvojic (M1; M2), chybovost dvou sekvencí LC312715.1 (HIV) a NC_001653.2 (hepatitis delta) je nižší než 10 % a dále klesá u uspořádaných dvojic (32; 0,60) a (30; 0,60). Sekvence MG772933_short se skládá převážně z části, která nemůže být pokryta čtením spojeným se SARS-CoV-2 (viz tabulky a obrázky: obrázek 3). Ani tady nebylo možné snížením hodnot M1 a M2 dosáhnout žádného zlepšení. Chybovost sekvencí NC_039345.1 (ebola), NC_024781.1 (Marburg), AF266291.1 a KJ410048.1 (spalničky) je výrazně vyšší než u LC312715.1 a NC_001653.2. Zatímco sekvence nukleových kyselin pro výpočet prvních genomů byly namnoženy ve Vero buňkách, sekvence nukleových kyselin použité pro LC312715.1 a NC_001653.2 pocházejí přímo ze vzorků lidského původu (tabulky a obrázky: tabulka 3). Vyvstává proto otázka, zda je tento výsledek způsoben strukturními rozdíly v příslušných zdrojích nukleových kyselin nebo použitými sekvenačními protokoly. Reverzní transkriptáza použitá k převodu RNA na cDNA nebo sekvence primerů použitých k amplifikaci, jakož i amplifikační cykly by mohly vést k rozdílům v získaných knihovnách sekvencí.

Nejvyšší chybovost F1 a F2 vykazují náhodně vygenerované fiktivní sekvence genomu rnd_uniform, rnd_wuhan, rnd_wh_mk_1 a rnd_wh_mk_2, takže tyto zjištěné výsledky nejsou čistě náhodné.

Grafická analýza rozložení pokrytí a délky čtení

Poté, co jsme pozorovali, že je možné vytvořit konsenzuální sekvence s vysokou kvalitou vzhledem k některým referenčním sekvencím, jsme analyzovali rozložení pokrytí asociovaných krátkých sekvenačních čtení (tabulky a obrázky: obrázky 1–22) a rozložení délek čtení (tabulky a obrázky: obrázky 23–25). Za tímto účelem jsme nejprve mapovali čtení krátkých sekvencí na příslušné referenční sekvence pomocí BBMap ((M1; M2) = (37; 0,60)). Kromě krátkých sekvencí jsme také mapovali 26 párů primerů [1, Supplementary Table 8. PCR primers used in this study] pro sekvenování celého genomu SARS-CoV-2 (GenBank: MN908947.3) na sledované referenční genomy. Následná analýza byla provedena prostřednictvím programu Tablet a tabulkového procesoru Excel.

Nejprve se podíváme na náhodně vygenerovanou referenci rnd_uniform. Srovnatelná pozorování platí pro náhodně generované referenční genomy rnd_wuhan, rnd_wh_mk_1 a rnd_wh_mk_2 (tabulky a obrázky: obrázky 14-16)

Pokrytí (rnd_uniform_reads) kolísá náhodně a relativně homogenně ve všech nukleotidových pozicích. Struktura je srovnatelná s náhodně generovaným pokrytím (Exponential distributed coverage – exponenciálně rozložené pokrytí), přičemž rozptyl se zdá být o něco nižší. Na několika izolovaných nukleotidových pozicích vykazuje pokrytí v porovnání s průměrem vysoké pokrytí, které se však vztahuje pouze na několik souvisejících nukleotidových oblastí. Souvislost s polohami primerů je jen slabě pozorovatelná. Čistě náhodně se projevující pokrytí krátkými sekvenčními čteními koreluje s nespojitě-mapovatelnou konsenzuální sekvencí a vysokou chybovostí F1 38,60 %. Náhodná (vnitřní) nukleotidová struktura stochasticky simulované referenční sekvence „rnd_uniform“ tedy není v sekvenačních datech, která jsou zde zkoumána, obsažena.

Naproti tomu se nyní podíváme na referenční genom SARS-CoV-2 (GenBank: MN908947.3).

Na rozdíl od obrázku 13 vykazuje rozložení pokrytí spíše vlnovitý vzorec s pravidelným, jasně zvýšeným pokrytím nukleotidů. 26 primerových párů je rozmístěno rovnoměrně na všech nukleotidových pozicích referenční sekvence. V blízkosti nukleotidových pozic s nukleotidovým pokrytím, které je ve srovnání s průměrem vysoké, se často nacházejí také pozice primerů. To ukazuje, že ne všechny části referenčního genomu byly amplifikovány stejnoměrně. Za předpokladu, že se všech 29 903 nukleotidových pozic ve čteních asociovaných se SARS-CoV-2 vyskytuje se stejnou pravděpodobností, by mělo pokrytí každé nukleotidové pozice ležet mezi oběma čarami s 99,5% pravděpodobností. To ale neplatí pro přibližně 90 % nukleotidových pozic. A priori by se dalo očekávat, že pokud je ve vzorku přítomno dostatečné množství virové RNA a je přečten dostatečný počet fragmentů celého řetězce, bude dosaženo homogenního pokrytí nukleotidů ve virovém genomu.

Následující graf umožňuje sledovat rozložení délek čtení právě zmíněných referencí (rnd_uniform a MN908947.3).

Obrázek 23e) ukazuje rozložení délek čtení v případě reference „rnd_uniform“. Průměrná délka čtení leží se svými 41,96 nt jen mírně napravo od maxima rozdělení. Pro srovnání, rozložení pro referenci MN908947.3, ukazuje obrázek 23a) nápadný (náhodný) rozsah podobný obrázku 23e) a výrazný rozsah se čteními o délce přibližně 150 nt. Průměrná délka čtení je více než 110 nt. Všechny referenční sekvence se srovnatelným, a tedy spíše náhodným rozložením délek čtení jako u stochasticky simulované referenční sekvence „rnd_uniform“ (tabulky a obrázky: obrázek 23d), f); Obrázek 24d), e), f); Obrázek 25a)-c)) vykazují rovněž vysokou chybovost F1 a F2 (tabulky a obrázky: tabulka 4).

Toto zjištění podtrhuje i následující analýza. Abychom lépe pochopili vnitřní strukturu publikovaných přibližně 56 milionů sekvencí, zohlednili jsme pro sekvenci MN908947.3 po provedeném mapování pomocí BBMap při tvorbě podskupiny kromě M1 a M2 také dodatečnou podmínku maxlength=100.

Vyloučením všech mapovatelných sekvencí delších než 100 nukleotidů bylo v podstatě odstraněno přibližně 120 000 čtení asociovaných se SARS-CoV-2. Rozložení pokrytí zbývajících krátkých sekvencí se nyní jeví jako náhodné, obdobně jako na obrázku 13. To opět koreluje s vysokou chybovostí R1 (29,90 %) a R2 (29,96 %). To ukazuje, že s výjimkou přibližně 120 000 asociovaných krátkých čtení (tabulky a obrázky: tabulka 1) není v publikovaných sekvencích obsažena žádná významná struktura referenčního souboru MN908947.3.

Než se budeme podrobněji věnovat některým zkoumaným referenčním genomům, rádi bychom se nejprve podívali na pokrytí dvou dalších kontigů (k141_12253 a k141_20796). Zatímco kontig s identifikací k141_12253 se vyznačuje poměrně vysokým pokrytím, patří k141_20796 k jedněm ze tří nejdelších vypočtených kontigů.

Kontig k141_12253 vykazuje vysokou podobnost s bakterií Leptotrichia (GenBank: CP012410.1). Z 52 publikovaných sekvencí primerů jich bylo možné 38 mapovat na referenci k141_12253 s poměrně vysokou chybovostí 37,30 %. Rozložení pokrytí je extrémně nehomogenní a vykazuje extrémně vysoké pokrytí nukleotidy ve srovnání s průměrem, zejména v prvních 500 nukleotidech. Oblasti s vysokým pokrytím korelují se zjištěnými primerovými pozicemi. To by mohlo naznačovat, že ve velkém množství nebyla amplifikována výhradně čtení, která jsou asociována se SARS-CoV-2. Vzhledem k poměrně vysoké chybovosti 37,30 % by to znamenalo poměrně nespecifickou amplifikaci. To vyvolává otázku, zda čtení získaná amplifikací cDNA pomocí specifických sekvencí primerů byla obsažena již v původním vzorku, nebo zda byla vytvořena samotným postupem.

Kontig k141_20796, který vykazuje vysokou podobnost s bakterií Veillonella parvula (GenBank: LR778174.1), ukazuje ve srovnání s kontigem s identifikací k141_12253 nižší pokrytí asociovanými čteními. Struktura nukleotidového pokrytí je podobná struktuře SARS-CoV-2 (GenBank: MN908947.3). Zejména pokrytí je opět nehomogenní, což svědčí o nerovnoměrné amplifikaci. Kvůli větší délce nukleotidů bylo nyní možné na referenční kontig mapovat 47 z 52 publikovaných sekvencí s průměrnou chybovostí 35,8 % primerů. Pozice primerů opět dobře korelují s oblastmi s vysokým pokrytím nukleotidů. To by také mohlo naznačovat nespecifickou amplifikaci sekvencí, které nejsou asociovány se SARS-CoV-2 (GenBank: MN908947.3).

V následující části se budeme podrobněji zabývat referenčními sekvencemi „Human immunodeficiency virus 1“ (GenBank: LC312715.1) a „Measles virus genotyp D8 strain MVi/Muenchen“ (GenBank: KJ410048.1). Všechny ostatní obrázky najdete v doplňujících materiálech (tabulky a obrázky: obrázky 1–22 a obrázky 23.–25).

Již předchozí část ukázala vysokou strukturní shodu publikovaných sekvencí s referenční sekvencí LC312715.1. Vypočtená konsenzuální sekvence vykázala ve srovnání např. s referencemi asociovanými se SARS relativně nižší chybovost R1 = 8,60 % a R2 = 8,83 %. Obrázek 6 ukazuje zřetelné rozdíly oproti obrázku 13. Rozložení pokrytí také vykazuje spíše vlnovitý vzorec s relativně pravidelnými oblastmi obzvláště vysokého pokrytí, a proto se zřetelně liší od rozložení pokrytí náhodného referenčního vzorku „rnd_uniform“. Rozložení délek čtení (obrázek 23b), srovnej také c)) se také výrazně liší od spíše náhodných rozložení a ukazuje signifikantní počet mapovatelných čtení s délkou do cca 110 nt. Průměrná délka čtení s 51,84 nt je také vyšší než například u „rnd_uniform“.

Zajímavá je zde také poloha primerových sekvencí ve vztahu k oblastem s vysokým pokrytím nukleotidy ve srovnání s průměrným pokrytím. Celkem 46 z 52 sekvencí primerů bylo možné s chybovostí 38,00 %přiřadit ke sledované referenci. Obrázek 6 naznačuje, že při PCR konfirmaci byla amplifikována i krátká sekvenční čtení asociovaná s referencí LC312715.1, přestože sekvence primerů je možné přiřadit k referencím pouze s poměrně vysokou chybovostí.

Nakonec se podívejme na referenci KJ410048.1 (virus spalniček).

Rozložení pokrytí se výrazně odlišuje od rozložení na obrázku 6 a vykazuje jistou podobnost s rozložením asociovaných sekvenačních čtení pro „rnd_uniform“, přičemž výkyvy v oblastech s nižším pokrytím jsou méně výrazné. Rozložení délek čtení (tabulky a obrázky: obrázek 24d)), stejně jako průměrná délka čtení 42,38 jsou srovnatelné s údaji „rnd_unifom“ a korelují také s poměrně vysokou chybovostí F1=28,70 % a F2=28,79 %.

Diskuze a závěr

Pomocí jednoduchého bioinformatického protokolu jsme prozkoumali publikované sekvenační údaje (přístupové číslo projektu BioProject PRJNA603194 v databázi NCBI Sequence Read Archive (SRA)) o sekvenci genomu SARS-CoV-2 (genová banka: MN908947.3). Námi použité metody nejsou specifické pro SARS-CoV-2 a lze je bez zvláštních úprav použít i na jiná sekvenační data.

Nejprve jsme zopakovali generování kontigů pomocí programu Megahit (verze 1.2.9) s použitím dostupných sekvenačních dat a získali jsme výrazně odlišné výsledky ve srovnání s daty uvedenými v [1]. Zejména se nám nepodařilo reprodukovat nejdelší kontig o délce 30 474 nt, který podle [1] zahrnuje téměř celý virový genom a posloužil jako základ pro design primerů. Naopak, nejdelší námi vygenerovaný kontig (29 802 nt) vykazoval téměř úplnou shodu s referencí MN908947.3. Z toho plyne, že publikovaná sekvenační data nemohou být původními krátkými čteními, která byla použita pro generování kontigů. A to je v kontextu vědeckých publikací velmi problematické, protože to znamená, že zveřejněné výsledky již nelze ověřit. Přičemž podstatou živé vědy je právě možnost ověřit publikované vědecké hypotézy.

Na rozdíl od skutečností uvedených v [1], jsme našli kontigy s vysokou hloubkou pokrytí, které jsou pravděpodobně asociovány s (ribozomálními) ribonukleovými kyselinami lidského původu. Při konstrukci SARS-CoV-2 tedy nebyly eliminovány všechny nukleové kyseliny spojené s lidmi. Dále nebyl předložen žádný důkaz o přítomnosti virových nukleových kyselin ve vzorku pacienta, a proto existuje možnost, že k vytvoření deklarované virové sekvence MN908947.3 byly do značné míry nevědomě použity fragmenty lidských nukleových kyselin. Tuto možnost by bylo třeba vyloučit kontrolními pokusy.

Ve všech publikacích týkajících se referenčních genomů analyzovaných v této studii nebyly rovněž poskytnuty potřebné důkazy o přesném původu sekvenovaných fragmentů, které byly použity pro konstrukci a nebyly zveřejněny ani potřebné kontrolní experimenty.

Rádi bychom na tomto místě také zmínili, že kontrolní experimenty již mohly být provedeny mnohokrát, a sice nepozorovaně, a poukazují na možnost konstrukce genomů SARS-CoV-2 z neinfekčních lidských vzorků. Například v [5] a [17] je zmínka o celogenomovém sekvenování ze vzorků s výchozí hodnotou Ct vyšší než 35. Což by mohla být skutečnost, která vyvrací virový model SARS-CoV-2.

Analýza rozložení nukleotidového pokrytí a rozložení délek mapovatelných sekvenčních čtení pro příslušné referenční sekvence vede k hypotéze o možné neúmyslné amplifikaci sekvenačních čtení, která nejsou asociována se SARS-CoV-2. Spolu s tím je třeba vzít v úvahu možnost neúmyslného generování sekvencí, které nebyly obsaženy v původním vzorku, ale vznikly pouze v důsledku podmínek amplifikace, jakými jsou např. použité sekvence primerů a počet provedených cyklů. Tato možnost proto vyžaduje provedení vhodných kontrolních experimentů.

Kromě pokusu o replikaci sestavy publikované v [1] se zveřejněnými sekvenačními čteními jsme sledovali jednoduchý protokol pro analýzu vnitřní struktury velkých souborů dat krátkých sekvenačních čtení. Díky dostupným sekvenačním datům jsme byli schopni vypočítat konsenzuální sekvence pro referenční genomy LC312715.1 (HIV) a NC_001653.2 (hepatitis delta) s vyšší kvalitou než pro referenční sekvence asociované s koronaviry, které jsme sledovali. To platí zejména pro bat-SL-CoVZC45 (GenBank: MG772933.1), který vedl k původní hypotéze o SARS-CoV-2. Tím se nám podařilo potvrdit naši hypotézu, že deklarované sekvence virového genomu jsou chybnou interpretací v tom smyslu, že byly nebo jsou nepozorovaně konstruovány z nevirových fragmentů nukleových kyselin. Naše výsledky zejména zdůrazňují naléhavou potřebu provést vhodné kontrolní experimenty. Pro každou genomovou sekvenci, u níž je podezření, že patří patogennímu viru, by zcela samozřejmý protokol spočíval v tom, že by se z odpovídajících nepodezřelých vzorků sestavily pomocí identických protokolů genomové sekvence.

U referenčních genomů spalniček, eboly nebo Marburgu, u nichž byly fragmenty nukleových kyselin, použité později ke konstrukci, namnoženy ve Vero buňkách, jsme zaznamenali vysokou chybovost R1 a R2. Zatím zůstává nezodpovězená otázka, zda je to způsobeno samotnými zdroji nukleových kyselin, nebo použitými amplifikačními podmínkami (např. primerovými sekvencemi a počtem cyklů) nebo sekvenačními protokoly (např. použitými polymerázami a reverzními transkriptázami).

V souvislosti s našimi výsledky doporučujeme kromě zveřejnění konečných použitých sekvenačních dat vždy také ta, která vznikla pouze amplifikací náhodnými hexamery a nízkým počtem cyklů, aby poskytla co nejméně zkreslená data pro strukturní analýzu.

Materiál a metody

Hloubka pokrytí jedné referenční sekvence krátkými sekvenačními čteními

Nechť G označuje délku referenční sekvence, ØL průměrnou délku čtení, n počet čtení krátké sekvence a N náhodnou průměrnou hloubku pokrytí referenční sekvence čteními krátké sekvence. Pak platí

Výraz lze chápat jako pravděpodobnost pokrytí nukleotidu v referenční sekvenci krátkým sekvenačním čtením.

Generování náhodných referenčních sekvencí

Následující věta umožňuje simulovat náhodnou veličinu X s kumulativní distribuční funkcí F.

Věta (princip inverze) [28]. Nechť U je náhodná veličina rovnoměrně rozdělená na intervalu (0,1). Dále nechť X je náhodná veličina s kumulativní distribuční funkcí , stejně jako

Pak platí

Nechť Ui,i=1,…,29,903 jsou nezávisle stejné rovnoměrně rozdělené náhodné veličiny na intervalu (0,1). Nechť pnt,nt ∈{A,T,C,G} označuje pravděpodobnost pro nukleotid nt. Pak nukleotid Ni,i=1,…,29.903 náhodně vygenerované referenční sekvence vznikne prostřednictvím

Pro referenční sekvenci „rnd_unifom“ bylo použito rovnoměrné rozdělení na množině {A,T,C,G}. Pro simulaci náhodné referenční sekvence „rnd_wuhan“ byl jako rozložení nukleotidů zvolen relativní výskyt nukleotidů A, T, C a G v sekvenci genomu SARS-CoV-2 (GenBank: MN908947.3). Při konstrukci randomizovaných referenčních sekvencí „rnd_wh_mk_1“ a „rnd_wh_mk_2“ byla zvolena podmíněná pravděpodobnost, která byla podmíněná posledním, resp. posledními dvěma nukleotidy podle odpovídajících empirických četností v sekvenci pro SARS-CoV-2 (GenBank: MN908947.3).

Stochastická simulace náhodných pokrytí jedné referenční sekvence

Kumulativní distribuční funkce exponenciálního rozdělení s parametrem λ je [28],

Nechť X je náhodná veličina s distribuční funkcí F. Pak platí

Bioinformatické metody (strukturní analýza)

1. Mapování pomocí BBMap

2. Výběr mapovaných sekvencí v závislosti na M1 a M2 s pomocí BBMap (reformat.sh)

3. Výpočet konsenzuální sekvence

• 3.1. Příprava pomocí Samtools

• 3.2 Hrubá verze konsenzuální sekvence

• 3.3 Finální konsenzuální sekvence (min. Q20)

4. Mapování konsenzuální sekvence na referenční sekvenci s pomocí BWA

5. Hodnocení pomocí programu Tablet

Hodnocení bylo provedeno pomocí vizualizačního softwaru Tablet, určeného pro vizualizaci sekvenačních dat.

---

Zdroje

• [1] Fan Wu u. a. A new coronavirus associated with human respiratory disease in China. In: Nature 580.7803 (2020). DOI: 10.1038/s41586-020-2202-3.

• [2] Na Zhu u. a. A Novel Coronavirus from Patients with Pneumonia in China, 2019. In: New England Journal of Medicine 382.8 (2020), s. 727–733. DOI:10.1056/nejmoa2001017.

• [3] Divinlal Harilal u. a. SARS-CoV-2 Whole Genome Amplication and Sequencing for Effective Population-Based Surveillance and Control of Viral Transmission. In: Clinical Chemistry 66.11 (2020), s. 1450–1458. DOI: 10.1093/clinchem/hvaa187.

• [4] Jalees A. Nasir u. a. A Comparison of Whole Genome Sequencing of SARSCoV-2 Using Amplicon-Based Sequencing, Random Hexamers, and Bait Capture. In: Viruses 12.8 (2020), s. 895. DOI: 10.3390/v12080895.

• [5] Clinton R. Paden u. a. Rapid, sensitive, full-genome sequencing of severe acute respiratory syndrome coronavirus 2. In: Emerging Infectious Diseases 26.10 (2020), s. 2401–2405. DOI: 10.3201/eid2610.201800.

• [6] Sureshnee Pillay u. a. Whole Genome Sequencing of SARS-CoV-2: Adapting Illumina Protocols for Quick and Accurate Outbreak Investigation during a Pandemic. In: Genes 11.8 (2020), s. 949. DOI: 10.3390/genes11080949.

• [7] Dan Hu u. a. Genomic characterization and infectivity of a novel SARS-like coronavirus in Chinese bats. In: Emerging Microbes & Infections 7.1 (2018), s. 1–10. DOI: 10.1038/s41426-018-0155-5.

• [8] Davaalkham Jagdagsuren u. a. The second molecular epidemiological study of HIV infection in Mongolia between 2010 and 2016. In: Plos One 12.12 (2017). DOI: 10.1371/journal.pone.0189605.

• [9] J. A. Saldanha, H. C. Thomas und J. P. Monjardino. Cloning and sequencing of RNA of hepatitis delta virus isolated from human serum. In: Journal of General Virology 71.7 (1990), s. 1603–1606. DOI: 10.1099/0022-1317-71-7-1603.

• [10] Jernej Mlakar u. a. Zika Virus Associated with Microcephaly. In: New England Journal of Medicine 374.10 (2016), s. 951-958. DOI: 10.1056 /nejmoa1600651.

• [11] Christopher L. Parks u. a. Comparison of Predicted Amino Acid Sequences of Measles Virus Strains in the Edmonston Vaccine Lineage. In: Journal of Virology 75.2 (2001), s. 910–920. DOI: 10.1128/jvi.75.2.910-920.2001.

• [12] Konstantin M. J. Sparrer u. a. Complete Genome Sequence of a Wild-Type Measles Virus Isolated during the Spring 2013 Epidemic in Germany. In: Genome Announcements 2.2 (2014). DOI: 10.1128/genomea.00157-14.

• [13] Paul A. Rota u. a. Characterization of a Novel Coronavirus Associated with Severe Acute Respiratory Syndrome. In: Science 300.5624 (2003), s. 1394–1399. DOI: 10.1126/science.1085952.

• [14] Runtao He u. a. Analysis of multimerization of the SARS coronavirus nucleocapsid protein. In: Biochemical and Biophysical Research Communications 316.2 (2004), s. 476–483. DOI: 10.1016/j.bbrc.2004.02.074.

• [15] Tracey Goldstein u. a. The discovery of Bombali virus adds further support for bats as hosts of ebolaviruses. In: Nature Microbiology 3.10 (2018), s. 1084–1089. DOI: 10.1038/s41564-018-0227-2.

• [16] Jonathan S. Towner u. a. Marburgvirus Genomics and Association with a Large Hemorrhagic Fever Outbreak in Angola. In: Journal of Virology 80.13 (2006), s. 6497–6516. DOI: 10.1128/jvi.00069-06.

• [17] Annika Brinkmann u. a. Amplicov: Rapid whole-genome sequencing using multiplex PCR amplication and real-time Oxford Nanopore minion sequencing enables rapid variant identication of SARS-COV-2. In: Frontiers in Microbiology 12 (2021). DOI: 10.3389/fmicb.2021.651151.

• [18] SARS-COV-2. url: https://artic.network/ncov-2019.

• [19] Ncbi. ncbi/sra-tools: SRA Tools. URL: https://github.com/ncbi/sra-tools.

• [20a] Dinghua Li u. a. MEGAHIT: an ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph. In: Bioinformatics 31.10 (2015), s. 1674–1676. DOI: 10.1093/bioinformatics/btv033.

• [20b] Voutcn. voutcn/megahit: Ultra-fast and memory-ecient (meta-)genome assembler. URL: https://github.com/voutcn/megahit.

• [21a] Shifu Chen u. a. fastp: an ultra-fast all-in-one FASTQ preprocessor. In: Bioinformatics 34.17 (2018), s. i884–i890. DOI: 10.1093/bioinformatics/bty560.

• [21b] OpenGene. OpenGene/fastp: An ultra-fast all-in-one FASTQ preprocessor (QC/adapters/trimming/ltering/splitting/merging…) URL: https://github. com/OpenGene/fastp.

• [22a] Ben Langmead u. a. Scaling read aligners to hundreds of threads on generalpurpose processors. In: Bioinformatics 35.3 (2018), s. 421–432. DOI: 10. 1093/bioinformatics/bty648.

• [22b] Ben Langmead. BenLangmead/bowtie2: A fast and sensitive gapped read aligner. URL: https://github.com/BenLangmead/bowtie2.

• [23a] Brian Bushnell. BBMap: A Fast, Accurate, Splice-Aware Aligner. In: (March 2014). URL: https://www.osti.gov/biblio/1241166.

• [23b] BBMap. url: https://sourceforge.net/projects/bbmap/.

• [24a] Li H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. In: (May 2013). URL: https://arxiv.org/abs/1303.3997.

• [24b] lh3. lh3/bwa: Burrow-Wheeler Aligner for short-read alignment (see mini-map2 for long-read alignment). URL: https://github.com/lh3/bwa.

• [25a] H. Li u. a. The Sequence Alignment/Map format and SAMtools. In: Bioinformatics 25.16 (2009), S. 2078-2079. DOI: 10.1093/bioinformatics/btp352.

• [25b] Samtools. url: http://www.htslib.org/.

• [25c] Petr Danecek u. a. Twelve years of SAMtools and BCFtools. In: GigaScience 10.2 (2021). DOI: 10.1093/gigascience/giab008.

• [25d] P. Danecek u. a. The variant call format and VCFtools“. In: Bioinformatics 27.15 (2011), s. 2156–2158. DOI: 10.1093/bioinformatics/btr330.

• [26] Tablet. URL: https://ics.hutton.ac.uk/tablet/.

• [27a] Wei Shen u. a. SeqKit: A Cross-Platform and Ultrafast Toolkit for FASTA/Q File Manipulation. In: Plos One 11.10 (2016). DOI: 10.1371/journal.pone.0163962.

• [27b] lh3. lh3/seqtk: Toolkit for processing sequences in FASTA/Q formats. URL: https://github.com/lh3/seqtk.

• [28] Albrecht Irle. Wahrscheinlichkeitstheorie und Statistik: Grundlagen – Resultate – Anwendungen. Teubner, 2010.