Costruzione di SILVA-derivati albero
Mentre l’originale SILVA albero è ben curata tassonomicamente, è principalmente destinato a essere utilizzato come una guida albero, e ri-calcolo del ramo lunghezza è generalmente consigliato per il downstream analisi filogenetiche ., Qui, per costruire un albero filogenetico con lunghezze di ramo più significative usando OTUs nel database SILVA non-redundant (NR99) 16S (release 128; ), abbiamo proceduto come segue. Le sequenze di SSU rappresentative allineate in SILVA sono state ridotte rimuovendo prima le posizioni dei nucleotidi con >95% di lacune e quindi rimuovendo il 5% più alto delle posizioni dei nucleotidi entropici. Le identità tassonomiche fornite da SILVA per OTUs a livello di dominio, phylum e classe sono state utilizzate per creare vincoli di divisione per FastTree, vincolando ogni taxon ad essere su un singolo lato di uno split e monofiletico., I taxa con meno di 10 OTU sono stati omessi dai vincoli. Sono stati così definiti un totale di 354 vincoli. Usando i vincoli generati tassonomicamente e prendendo l’albero SILVA originale come albero di partenza, abbiamo costruito un albero filogenetico dagli allineamenti ridotti con FastTree v2.1.10 (opzioni “-spr 4-gamma-fastest-no2nd-constraintWeight 100”). L’albero filogenetico è stato riavviato in modo che batteri e Archaea siano divisi alla radice. Il nostro albero derivato da SILVA è fornito come file aggiuntivo 2., Per tutte le analisi a valle, i cloroplasti, i mitocondri e gli Eucarioti sono stati omessi dall’albero. Nell’articolo principale, descriviamo le nostre analisi usando questo albero derivato da SILVA (Fig. 1); risultati analoghi per l’albero SILVA originale sono mostrati nel file aggiuntivo 1: Figura S1.
Distribuzione filogenetica di 16S GCNs
Per esaminare come 16S GCNs sono distribuiti filogeneticamente e per valutare la loro prevedibilità generale utilizzando vari metodi filogenetici, abbiamo proceduto come segue., Un totale di 8.767 genomi batterici e archaeal annotati con stato di completamento “Genoma completo” sono stati scaricati dal database NCBI RefSeq il 4 gennaio 2018. I genomi scaricati sono stati controllati per la potenziale contaminazione usando checkM 1.0.6 (opzione “reduced_tree”), che si basa sulla rilevazione di geni marcatori conservati (riassunti di assemblaggio e checkM nel file aggiuntivo 3). I genomi trovati per esibire un livello di contaminazione superiore a 1% o un’eterogeneità del ceppo superiore a 1% sono stati scartati, lasciandoci con 6.868 genomi completi per l’analisi a valle (file aggiuntivo 4).,
Per ogni genoma, i GCN 16S sono stati determinati utilizzando due approcci: In primo luogo, abbiamo contato il numero di sequenze rRNA 16S annotate nelle annotazioni NCBI (file rna_from_genomic.Fna). In secondo luogo, abbiamo usato modelli di covarianza con il programma cmsearch (come parte della versione INFERNALE 1.1.2, opzioni “no noali cut cut_nc” ” per cercare sequenze rRNA 16S all’interno dei genomi assemblati (file genomici.Fna). Modelli di covarianza separati per i geni rRNA 16S archaeal e batterici sono stati ottenuti dal database Rfam (accessions RF00177 e RF01959)., Una tabella che elenca i GCN calcolati utilizzando entrambi i metodi viene fornita come file aggiuntivo 5. Solo i genomi per i quali i due metodi hanno prodotto gli stessi GCN 16S sono stati considerati per l’analisi successiva, producendo GCN 16S per i genomi 6,780 (“genomi di alta qualità”, file aggiuntivo 6). L’accuratezza di questi GCN è stata ulteriormente verificata attraverso il confronto con il Ribosomal RNA Operon Copy Number Database (rrnDB, accessibile il 7 giugno 2017; ) ogni volta che un’adesione all’assemblaggio del genoma era presente nell’attributo rrnDB (rrnDB “Data source record id”)., Attraverso 5.616 genomi di alta qualità testati, abbiamo trovato un accordo quasi perfetto con rrnDB (R2>0.999; File aggiuntivo 1: Figura S2). I riassunti di qualità checkM per il set di genoma di alta qualità sono forniti come file aggiuntivo 7.
I suggerimenti sull’albero derivato da SILVA sono stati mappati su genomi di alta qualità, quando possibile, come segue: In primo luogo, le sequenze rappresentative 16S di SILVA OTUs sono state allineate alla sequenza rRNA 16S più lunga di ciascun genoma utilizzando vsearch 2.3.,4 al massimo (100%) somiglianza (vsearch opzioni “–strand sia –usearch_global –maxaccepts 0 –top_hits_only –iddef 0 –id 1.0”). Se un OTU allineato a più genomi, tutti i genomi sono stati inizialmente tenuti. Successivamente, per ogni coppia OTU-genoma allineata, abbiamo confrontato l’ID del taxon NCBI (”taxid”) dell’OTU con quello del genoma. I taxid OTU sono stati ottenuti da una tabella di ricerca fornita da SILVA (https://www.arb-silva.de/fileadmin/silva_databases/release_128/Exports/taxonomy/taxmap_embl_ssu_ref_128.txt.gz). I tassidi del genoma sono stati ottenuti da tabelle di ricerca fornite da NCBI (ftp://ftp.ncbi.nlm.nih.gov/genomes/genbank/* / assembly_summary.txt, dove “*” è “batteri” o “archaea”)., Qualsiasi coppia OTU-genoma allineata con tassidi non identici è stata omessa. Delle restanti coppie OTU-genoma con tassidi identici, abbiamo mantenuto solo il primo genoma allineato per ogni OTU. Un totale di 9.395 OTU potrebbe quindi essere mappato su uno dei genomi. Per ogni OTU mappato, abbiamo assunto un GCN uguale al GCN contato per il genoma corrispondente. Per tutti gli altri OTU, abbiamo assunto un GCN sconosciuto.
Tutte le analisi filogenetiche sono state eseguite utilizzando il pacchetto R castor , disponibile presso il CRAN (Comprehensive R Archive Network). NSTD per tutti i suggerimenti rispetto ai suggerimenti mappati su un genoma sequenziato (Fig., 1b) sono stati calcolati utilizzando la funzione di ricino find_nearest_tips. La funzione filogenetica di autocorrelazione (ACF)dei noti 16S GCN attraverso l’albero derivato da SILVA (Fig. 1a) è stato calcolato utilizzando la funzione castor get_trait_acf basata su 108 coppie di punte (opzioni “Npairs=1e8, Nbins=100”), scelte casualmente tra punte con GCN noto. La funzione get_trait_acf sceglie casualmente le coppie OTU sull’albero, le mette in uno dei tanti intervalli di distanza filogenetica e calcola l’autocorrelazione di Pearson tra GCN delle coppie OTU all’interno di ciascun bin., Si noti che questa analisi non presuppone che GCNs scala linearmente con la distanza filogenetica. Invece, l’ACF misura semplicemente la correlazione statistica tra GCN su punte distinte, a condizione che le punte siano entro una certa distanza filogenetica l’una dall’altra.,
GCNs sono stati ricostruiti su SILVA-derivati albero utilizzando Sankoff massima parsimonia (funzione hsp_max_parsimony, con opzione transition_costs sia impostato su “esponenziale” “proporzionale” o “all_equal”), filogenetica indipendente contrasti (funzione hsp_independent_contrasts), ponderato a-squared-cambiare parsimonia (funzione hsp_squared_change_parsimony), sottostruttura in media (funzione hsp_subtree_averaging), e di massima verosimiglianza di Mk modelli con rerooting (funzione hsp_mk_model_rerooting con opzioni root_prior= “empirica”, optimization_algorithm=‘nlminb’, Ntrials=5, rate_model=‘ER’).,
Per calcolare la frazione validata incrociata della varianza prevista da (aka. coefficiente di determinazione cross-validated di) ogni metodo (\(R^{2}_{\text {cv}}\);) in funzione del NSTD (Fig. 1c), abbiamo proceduto come segue. Abbiamo scelto casualmente il 2% dei suggerimenti con GCN 16S noti da escludere dall’input alle ricostruzioni e da utilizzare come “set di test” indipendente in seguito. A seconda del taglio NSTD considerato (ad esempio 10% di sostituzioni per sito), abbiamo anche escluso tutti i suggerimenti la cui distanza filogenetica dal set di test era inferiore al taglio NSTD., I restanti suggerimenti con GCN noti (”training set”) sono stati utilizzati come input per le ricostruzioni e i GCN previsti per il set di test sono stati quindi confrontati con i GCN noti del set di test. Questo processo è stato ripetuto tre volte e la R2 risultante è stata calcolata in media su tutte le ripetizioni, ottenendo un \(R^{2}_{\text {cv}}\) per ogni cutoff NSTD considerato. Lo script R per analizzare e ricostruire GCN 16S attraverso l’albero derivato da SILVA è disponibile come file aggiuntivo 8. Per confronto, tutte le analisi di cui sopra sono state eseguite anche utilizzando l’albero guida SILVA originale (file aggiuntivo 1: Figura S1).,
Valutazione di strumenti di previsione GCN di terze parti su genomi sequenziati
Per testare l’accuratezza predittiva di CopyRighter , PICRUSt e PAPRICA per genomi con GCN noti , abbiamo confrontato le loro previsioni con i GCN contati nei genomi sequenziati (di alta qualità). Per valutare l’accuratezza predittiva di CopyRighter sui genomi, abbiamo proceduto come segue: Abbiamo prima scaricato la tabella di ricerca precalcolata che elenca le previsioni di CopyRighter per il database rRNA di Greengenes 16S (rilascio ottobre 2012, “GG2012”; ), dal Github del progetto il 6 giugno 2017 (v0.,46):https://github.com/fangly/AmpliCopyRighter (CopyRighter-0.46/data/201210/ssu_img40_gg201210.txt). Poi abbiamo allineato la più lunga sequenza 16S rRNA di ogni genoma OTUs (cluster al 99% di similarità) in Greengenes database utilizzando vsearch (vsearch opzioni “–strand sia –usearch_global –maxhits 1 –maxaccepts 10 –top_hits_only”), scegliendo sempre i migliori match in Greengenes e tenendo solo i genomi che mappato a un Greengenes voce da parte di almeno il 99% di somiglianza (5688 genomi mappati)., Per ogni genoma mappato, abbiamo preso il GCN previsto da CopyRighter per la voce Greengenes corrispondente come previsione del copyright per il genoma. Questa previsione è stata quindi confrontata con il GCN contato dalla sequenza del genoma. Un istogramma delle previsioni di CopyRighter attraverso genomi mappati è mostrato nel file aggiuntivo 1: Figura S4B. L’accuratezza predittiva di CopyRighter è stata misurata in termini di frazione di varianza spiegata (R2), in funzione di NSTD di un genoma (Fig. 1 bis). NSTD di genomi sono stati calcolati come descritto in una sezione separata di seguito.,
Un approccio simile è stato utilizzato per PICRUSt : la tabella di ricerca precalcolata che elenca le previsioni di PICRUSt per il database Greengenes (release May 2013; “GG2013”) è stata scaricata dal sito Web del progetto il 6 giugno 2017 (v1.1.1):https://picrust.github.io/picrust/picrust_precalculated_files.html(16S_13_5_precalculated.tab.gz). Un totale di genomi di alta qualità 5,708 potrebbe essere mappato a un OTU (somiglianza 99%) in GG2013. Un istogramma delle previsioni di PICRUSt su tutti i genomi mappati è mostrato nel file aggiuntivo 1: Figura S4C. L’accuratezza predittiva di PICRUSt è stata misurata in termini di R2 in funzione dell’NSTD di un genoma (Fig., 1b), analogamente a CopyRighter.
Per valutare l’accuratezza predittiva di PAPRICA sui genomi, abbiamo proceduto come segue: Abbiamo prima scaricato e installato PAPRICA dal Github del progetto il 6 giugno 2017 (v0.4.0 b): https://github.com/bowmanjeffs/paprica. Questa versione include alberi di riferimento precalcolati (uno per archaea e uno per batteri) e tabelle che elencano 16 GCN per i genomi di calibrazione dello strumento rappresentati negli alberi di riferimento. Abbiamo usato la più lunga sequenza 16S rRNA da ogni genoma come input per la pipeline PAPRICA (comando “paprica-run.sh”), separatamente per archaea e batteri., La pipeline produce, tra gli altri, una tabella che elenca l’abbondanza non corretta di ogni sequenza di input univoca (questa può essere maggiore di 1 se più genomi condividono la stessa sequenza rRNA 16S) e l’abbondanza corretta corrispondente (dopo aver diviso per il GCN 16S previsto). Abbiamo usato questa tabella per ottenere i GCN 16S previsti da PAPRICA per le sequenze 16S uniche (che rappresentano 3473 sequenze 16S), dividendo il non corretto per l’abbondanza corretta. Abbiamo quindi confrontato questi GCN previsti con i GCN contati nelle sequenze del genoma, in modo simile a quanto sopra., Un istogramma delle previsioni di PAPRICA su tutti i genomi rappresentati è mostrato nel file aggiuntivo 1: Figura S4D. L’accuratezza predittiva di PAPRICA è stata misurata in termini di R2 in funzione dell’NSTD di un genoma (Fig. 1a), in modo simile a CopyRighter.
Confronto degli strumenti di previsione GCN di terze parti su Greengenes
Per confrontare le previsioni di CopyRighter con quelle di PICRUSt su tutti gli OTU in Greengenes (Fig. 3a), per prima cosa abbiamo mappato tutti gli OTU in GG2013 a OTU in GG2012 usando vsearch (con le opzioni ” strand strand entrambi us usearch_global”)., Abbiamo mantenuto solo le partite al 100% di somiglianza (153.375 su 203.452 OTU in GG2013). Per ogni OTU mappato in GG2013, abbiamo confrontato il GCN corrispondente previsto da PICRUSt con il GCN previsto da CopyRighter per l’OTU corrispondente in GG2012. Per calcolare le distribuzioni di frequenza di GCN previste da CopyRighter e PICRUSt su tutte le OTU in Greengenes (istogrammi nel file aggiuntivo 1: Figura S3A, B), abbiamo usato le GCN elencate nelle loro tabelle di ricerca precalcolate.
Per confrontare PAPRICA a PICRUSt attraverso Greengenes (Fig., 3b), abbiamo proceduto come segue: Le sequenze rappresentative di OTU in GG2013 sono state suddivise in sequenze archaeal e batteriche. Ogni file fasta risultante è stato utilizzato come input per la pipeline PAPRICA per prevedere il corrispondente GCN 16S, come descritto sopra per i genomi. Ciò ha prodotto un GCN previsto per tutte le voci di Greengenes. Queste previsioni sono state confrontate con i valori GCN precalcolati forniti da PICRUSt. Queste previsioni sono state anche utilizzate per calcolare la distribuzione di frequenza dei GCN predetti da PAPRICA attraverso i Greengenes (File aggiuntivo 1: Figura S3C). Per confrontare il copyright con PAPRICA (Fig., 3c), abbiamo proceduto come descritto sopra per il confronto di CopyRighter a PICRUSt.
Confronto degli strumenti di previsione GCN di terze parti tra comunità microbiche
Per confrontare CopyRighter, PICRUSt e PAPRICA tra OTU in varie comunità microbiche, abbiamo proceduto come segue. I dati di sequenza di amplicon 16S rRNA disponibili pubblicamente da vari campioni ambientali sono stati scaricati dall’Archivio nucleotidico europeo (http://www.ebi.ac.uk/ena). Sono stati presi in considerazione solo i dati di sequenza Illumina ottenuti da ampliconi ottenuti utilizzando primer sensibili ai batteri e/o agli archaea., I campioni sono stati scelti per coprire una vasta gamma di ambienti, tra cui i sedimenti oceanici, marini e lacustri, il suolo, i laghi salini e ipersalini, le bocche idrotermali, le sorgenti calde, i bioreattori e i microbiomi associati agli animali. Tutti i dati di sequenziamento sono stati elaborati in modo simile, ove possibile, come segue. Le letture accoppiate sovrapposte sono state unite usando flash v1.2. 11 (opzioni –min-overlap=20 –max-overlap=300 –max-mismatch-density 0.25 –phred-offset=33 –allow-outies) e le letture accoppiate non sovrapposte sono state omesse. Le letture single-end sono state mantenute invariate., Tutte le letture single-end e le letture accoppiate unite sono state quindi filtrate dalla qualità utilizzando vsearch v2.4.3 (opzioni –fastq_ascii 33 –fastq_minlen 120 –fastq_qmin 0 –fastq_maxee 1 –fastq_truncee 1 –fastq_maxee_rate 0.005 –fastq_stripleft 7). I campioni con più di 20.000 letture filtrate dalla qualità sono stati rarefatti fino a 20.000 letture per ridurre il tempo di calcolo, selezionando casualmente le letture senza sostituzione., Le sequenze filtrate dalla qualità sono state raggruppate in unità tassonomiche operative (OTU; con una somiglianza del 97%) allineando il riferimento globale chiuso contro il database di riferimento SILVA SSU non ridondante (NR99) (release 128;), usando vsearch. Entrambi i trefoli sono stati considerati per l’allineamento (opzione vsearch both strand entrambi). Le sequenze che non corrispondono a nessuna voce del database con una somiglianza del 97% o superiore sono state scartate. Si noti che gli OTU erano quindi rappresentati dalle voci SILVA, vale a dire quelle utilizzate per seminare i cluster. Sono stati omessi cloroplasti, mitocondri e qualsiasi Eucariota., OTU rappresentato da meno di cinque letture in tutti i campioni sono stati omessi. Infine, tutti i campioni con meno di 2.000 letture contabilizzate da OTU sono stati omessi. Ciò ha prodotto una tabella OTU con 635 campioni e 65.673 OTU rappresentati da 4.827.748 letture (in media 734 OTU per campione). I numeri di adesione del campione, le coordinate, le date di campionamento, le pubblicazioni originali, le piattaforme di sequenziamento, le lunghezze di lettura filtrate dalla qualità e i conteggi di lettura e le regioni di primer coperte (se disponibili) sono forniti nel file aggiuntivo 9.,
Per prevedere GCNs per OTUs in ogni campione utilizzando CopyRighter, abbiamo utilizzato lo stesso approccio per i genomi: Rappresentante 16S sequenze di OTUs sono stati allineati a GG2012 utilizzando vsearch (opzioni “–strand sia –usearch_global –iddef 0 –id 0.99 –maxhits 1 –maxaccepts 10 –top_hits_only”), omettendo qualsiasi OTUs non abbinato a una Greengenes voce da parte di almeno il 99% di somiglianza. Per ogni OTU mantenuto, il GCN elencato da CopyRighter per la voce Greengenes abbinato è stato preso come previsione del copyright. Per PICRUSt, abbiamo proceduto in modo analogo, utilizzando GG2013 invece di GG2012., Per PAPRICA, abbiamo proceduto in modo analogo, utilizzando le previsioni GCN di PAPRICA calcolate in precedenza per GG2013 (vedi sezione precedente).
Per confrontare due strumenti dati (CopyRighter vs. PICRUSt, PICRUSt vs. PAPRICA, o CopyRighter vs. PAPRICA) per un campione specifico, sono stati considerati solo OTU con almeno una lettura nel campione e con una previsione GCN da entrambi gli strumenti. Abbiamo misurato l’accordo tra due strumenti in termini di frazione di varianza nelle previsioni del 1 ° strumento che è stato spiegato dalle previsioni del 2 ° strumento (R2)., Abbiamo calcolato l’NSTI del campione (nearest sequenced taxon index) secondo , cioè, come la media aritmetica NSTD su tutti gli OTU considerati nel confronto e ponderati per le frequenze OTU relative. I dettagli su come sono stati calcolati gli NSTD sono forniti nella sezione seguente. Per ogni coppia di strumenti confrontati, abbiamo così ottenuto 635 NSTIs e 635 R2 su 635 campioni, mostrati in Fig. 4. I coefficienti di correlazione di Pearson (r2) tra NSTIs e R2 sono stati calcolati per ciascuna coppia di strumenti, separatamente per campioni associati agli animali e non agli animali., Le significanze statistiche (valori P) dei coefficienti di correlazione sono state stimate utilizzando un test di permutazione con 1000 permutazioni. File aggiuntivo 1: Figure S6 e S7 mostrano GCN previsti da ogni strumento per varie comunità microbiche. Mostriamo anche deviazioni relative tra gli strumenti (|A−B|/((A+B)/2), dove A e B sono GCN previsti da due strumenti per lo stesso OTU) e NSTD per OTU in vari campioni (File aggiuntivo 1: Figura S8).,
Valutazione e confronto degli strumenti di previsione GCN a seconda di NSTD
Per esaminare l’accuratezza predittiva di CopyRighter, PICRUSt e PAPRICA in funzione dell’NSTD di un OTU o del genoma, abbiamo proceduto come segue. Per ogni OTU in SILVA, e separatamente per ogni strumento, abbiamo calcolato il NSTD come la distanza filogenetica al genoma sequenziato più vicino utilizzato dallo strumento per fare previsioni (“genomi di calibrazione”). Per PAPRICA, è stato ottenuto un elenco di 5.628 genomi di calibrazione dai file precalcolati di PAPRICA (PAPRICA/ref_genome_database/*/genome_data.finale.,csv, dove ” * ” è o batteri o archaea). I genomi di calibrazione sono stati abbinati a SILVA OTUs tramite l’allineamento globale del gene 16S a una soglia di somiglianza del 99%, utilizzando vsearch. Si presume che gli OTU abbinati abbiano un NSTD uguale a zero, e per tutti gli altri OTU SILVA, l’NSTD è stato calcolato in base all’albero derivato da SILVA e utilizzando il castor del pacchetto R. Una corrispondenza approssimativa di genomi a OTU (cioè.,, al 99% di somiglianza) è stato scelto per garantire che il maggior numero possibile di genomi di calibrazione siano inclusi; si noti che SILVA OTUS sono essi stessi raggruppati a tale somiglianza e che l’errore potenzialmente introdotto per NSTD e NSTIs è trascurabile (< 1% sostituzioni nucleotidiche per sito). Per PICRUSt, è stata scaricata una tabella dal sito web del progetto che elenca gli ID IMG (Integrated Microbial Genomes) per 2.887 genomi di calibrazione (https://github.com/picrust/picrust/tree/master/tutorials/picrust_starting_files.zip, file GG_to_IMGv350.txt). Gli ID IMG sono stati tradotti in ID sequenza GG2013 utilizzando gg_13_5_img.,tabella di ricerca txt scaricata dal sito web di Greengenes (http://greengenes.secondgenome.com/downloads). Gli ID GG2013 abbinati sono stati quindi mappati a SILVA OTUs tramite l’allineamento della sequenza 16S globale con vsearch, a una soglia di somiglianza del 99%. Gli NSTD di SILVA OTUs sono stati quindi calcolati allo stesso modo di PAPRICA. Per CopyRighter, una tabella di ricerca è stata scaricata dalla pagina Github del progetto che mappa i genomi di calibrazione alle sequenze GG2012 (https://github.com/fangly/AmpliCopyrighter, file AmpliCopyrighter-0.46/preprocessing/ data/img_to_gg.txt)., Le sequenze GG2012 elencate in quella tabella sono state mappate su SILVA OTU e gli NSTD sono stati calcolati per tutti SILVA OTU, in modo simile a PICRUSt. Per determinare gli NSTD per i genomi esaminati in questo studio (separatamente per CopyRighter, PICRUSt e PAPRICA), i genomi sono stati mappati a SILVA OTUs tramite l’allineamento globale della loro sequenza 16S più lunga disponibile al 99% di somiglianza. Per ogni genoma, l’NSTD del SILVA OTU più strettamente abbinato è stato preso come NSTD del genoma. Per determinare NSTD per tutti Greengenes OTUs, abbiamo mappato Greengenes OTUs a SILVA OTUs tramite allineamento globale al 99% di somiglianza., Per determinare gli NSTD per OTU recuperati dalle comunità microbiche campionate, abbiamo utilizzato direttamente gli NSTD di SILVA OTUS usati come semi durante la raccolta di OTU di riferimento chiuso. Quando si confrontano due strumenti di previsione GCN su un OTU (ad esempio, Figs. 3 e 4 e file aggiuntivo 1: Figura S8), nei casi in cui i due NSTD differivano, abbiamo usato la loro media aritmetica. Per calcolare l’R2 tra due strumenti di previsione GCN o tra uno strumento di previsione GCN e i “GCN veri”, in funzione dell’NSTD (Figg., 2 e 3d-f), abbiamo binned gli OTU o genomi utilizzati nel confronto in intervalli NSTD di dimensioni uguali e calcolato il R2 separatamente per ogni intervallo. Sono stati considerati solo intervalli NSTD con almeno 10 OTU o genomi.
Lascia un commento