Information

Was ist der Unterschied zwischen Contig- und Read-basiertem Sequenz-Alignment?

Was ist der Unterschied zwischen Contig- und Read-basiertem Sequenz-Alignment?


We are searching data for your request:

Forums and discussions:
Manuals and reference books:
Data from registers:
Wait the end of the search in all databases.
Upon completion, a link will appear to access the found materials.

Ich versuche, den Unterschied zwischen lesebasierter und Contig-basierter Ausrichtung zu verstehen. Bezieht sich die Ausrichtung auf Contig-Basis auf eine De-novo-Assemblierung und dann auf ein Referenzgenom. Ich bin verwirrt, wenn Sie die basierte Ausrichtung gelesen haben, welche Bedeutung es dann hat, Contigs an der Referenz auszurichten.


Ich habe den Begriff „Contig-basiertes Alignment“ noch nie gehört, und Ihre Frage ist der einzige Google-Treffer dieser genauen Abfrage (abgesehen von einer Patentanmeldung von 2012).

Trotzdem und ohne den genauen Kontext zu kennen, gehe ich davon aus, dass Sie im Wesentlichen Recht haben: Contig-basiertes Alignment bezieht sich wahrscheinlich auf die de novo Zusammenstellung von Reads in Contigs, die dann mit Hilfe einer Referenz zu einem Scaffold ausgerichtet werden.


Vergleiche zwischen SRST2, ARIBA und KmerResistance

SRST2 1 , ARIBA 2 und KmerResistance (Webservice, Code) 3 sind drei weit verbreitete eigenständige Software für den lesebasierten Nachweis von Zielgenen in bakteriellen Genomen. SRST2 wurde 2014 veröffentlicht und gilt als Pionier unter diesen drei Tools 2, 3 . In diesem Beitrag vergleiche ich die diesen Tools zugrunde liegenden Methoden kurz und bündig, um die Auswahl geeigneter Software für den Gennachweis zu beleuchten. Insbesondere gehe ich hier davon aus, dass der Nachweis von antimikrobiellen Resistenzdeterminanten der einzige Anwendungsfall ist. Sofern nicht angegeben, sind in diesem Beitrag folgende Softwareversionen: SRST2 v0.2.0, ARIBA 2.14.4 und KmerResistance v2.2.


Hintergrund

Die Metagenomik ist ein leistungsstarker und sich schnell entwickelnder Ansatz, der verwendet werden kann, um die unkultivierte mikrobielle Vielfalt zu entschlüsseln und den Baum des Lebens zu erweitern und neue biologische Einblicke in die Mikroben zu geben, die in wenig erforschten Umgebungen leben [1]. Metagenomik liefert sowohl im gastrointestinalen (GI) als auch im fäkalen Mikrobiom des Hundes Informationen über Gesundheit und Krankheit sowie wesentliche Hinweise zur Vorbeugung oder Behandlung bestimmter Pathologien.

Frühere Studien haben Ähnlichkeiten zwischen dem gastrointestinalen Mikrobiom von Hunden und Menschen berichtet. Im Allgemeinen beziehen sich verschiedene GI-Erkrankungen auf ein verändertes GI-Mikrobiom, das andererseits durch Ernährung und Nahrungsergänzungsmittel (wie Prä- und Probiotika) moduliert werden kann (siehe [2,3,4,5] für ausführliche Übersichten) . Abgesehen vom tierärztlichen Interesse selbst gelten Hunde als nähere Modelle für den Menschen als andere Tiermodelle für GI-Mikrobiomstudien [6, 7].

Mikrobiomstudien sind überwiegend entweder markerspezifisch (z. B. 16S-rRNA-Gen für Bakterien) oder Sequenzierung des gesamten Metagenoms [8]. Bis heute verwenden die verfügbaren GI-Mikrobiomstudien des Hundes Next-Generation-Sequenzierung – Short-Read-Sequenzierung – oder frühere Technologien und sind meist amplikonbasierte Strategien (16S rRNA-Gen). Nur drei Studien verwendeten Shotgun-Metagenomik mit Short-Read-Sequenzierung, um die gesamte mikrobielle Gemeinschaft und den Gengehalt im Hundekot zu charakterisieren [7, 9, 10].

Die Anwendung der Long-Read-Sequenzierung auf die Metagenomik ermöglicht das Abrufen von Metagenom-Assembly-Genomen (MAGs) mit hoher Vollständigkeit. Die neueste Strategie in der Long-Read-Metagenomik verwendet die Long-Reads, um den Entwurf der Metagenom-Assembly zu erhalten – um die größte Kontiguität von MAGs zu gewährleisten – und Short-Reads, um die Gesamtgenauigkeit zu verbessern und zu verbessern. Diese Strategie wurde unter anderem angewendet, um das menschliche GI-Mikrobiom [11] zu untersuchen – wie zum Beispiel Scheingemeinschaften [12], Kuhpansen [13], natürliche Molkestarterkulturen [14] oder Abwasser [15]. Einige Autoren schlagen vor, dass wir die Notwendigkeit von kurzen Lesevorgängen überwinden können, um lange gelesene Daten zu polieren, indem wir entweder Korrektursoftware verwenden, wie z. .

In unserer vorherigen Arbeit haben wir langgelesene Metagenomik verwendet, um die Taxonomie zu bewerten und die Artidentifikation im Hundekotmikrobiom zu erreichen. Obwohl wir einen Sequenzierungsansatz mit geringer Tiefe verwendeten, stellten wir ein kreisförmiges Contig zusammen, das einem entspricht unkultivierter CrAsphage [17].

In der vorliegenden Studie verwenden wir Nanoporen-Long-Read-Metagenomik und Frameshift-bewusste Korrektur, um die Notwendigkeit des Polierens mit Short-Reads zu überwinden. Als Ergebnis gewinnen und charakterisieren wir acht hochwertige MAGs und gewinnen neue biologische Einblicke in das Hundekot-Mikrobiom.


Ergebnisse

Genetische Vielfalt und Populationseigenschaften

Wir sammelten DNA-Resequenzierungsdaten für 1961er Baumwolle für eine genomische Variationsanalyse mit einer durchschnittlichen Tiefe von

14,8× für jeden [3,4,5,6, 16, 33, 34]. Nachdem doppelte Akzessionen verworfen wurden, wurden insgesamt 1913 Baumwollakzessionen für die SNP- und InDel-Analyse verwendet, darunter 256 G. hirsutum Landrassen (Ghlandrassen), 438 verbessert G. hirsutum Sorten aus den USA und anderen Ländern (GhImpUSO), 929 verbessert G. hirsutum Sorten aus China (GhImpCHN), 261 G. barbadense Beitritte und 29 weitere Gossypium Arten, die als Fremdgruppe verwendet wurden (Zusatzdatei 1: Tabelle S1). Wir haben diese Daten mit dem Referenzgenom von abgeglichen G. hirsutum gem. „TM-1“ [12] und identifizierte 63.084.975 SNPs und 12.354.432 kleine Insertionen oder Deletionen (InDels-Länge ≤ 20 bp), wobei der Kernvariationsdatensatz 19.246.497 SNPs und 4.815.125 InDels mit einer Minor Allel Frequency (MAF) ≥ 0,01 und mehr umfasst als fünf Akzessionen mit homozygoten Variationen (Tabelle 1 Zusatzdatei 1: Tabellen S2-S6 Zusatzdatei 3). Basierend auf SNP-Kerndaten untersuchten wir die Bevölkerungsstruktur von G. hirsutum und G. barbadense. Die Analyse von Nachbarbäumen zeigte, dass die Beitritte von 1913 in 12 Kladen eingeteilt wurden. G. hirsutum Beitritte bilden 8 Kladen, G. barbadense Akzessionen bilden 3 Kladen und andere Arten bilden 1 Klade (Abb. 1a Zusatzdatei 2: Abbildung S1). Die Bevölkerungsanalyse hat gezeigt, dass G. barbadense Beitritte wurden von den G. hirsutum Landrassen, GhImpUSO und GhImpCHN (Abb. 1b, c Zusatzdatei 2: Abbildung S2). G. hirsutum Nukleotidvielfalt (π) wird auf 1,07 × 10 – 3 in Landrassen, 3,74 × 10 – 4 in GhImpUSO, 3,34 × 10 – 4 in GhImpCHN und 1,01 × 10 – 3 in . geschätzt G. barbadense (Zusatzdatei 2: Abbildung S3), ähnlich den neueren Studien an Baumwolle [3,4,5,6, 34] (Abb. 1d).

Bevölkerungsstruktur und genetische Vielfalt in G. hirsutum und G. barbadense Beitritte. ein Der ungewichtete Neighbor-Joining-Stammbaum von 1913 Baumwollakzessionen wurde basierend auf 20.000 zufälligen SNPs von Kern-SNPs konstruiert. Die G. tomentosum (ANZEIGE3), G. mustelinum (ANZEIGE4), G. darwinii (ANZEIGE5), G. ekmanianum (ANZEIGE6), G. stephensii (ANZEIGE7) von tetraploiden Arten, G. arboreum (EIN2) und G. davidsonii (D3-d) diploider Arten dienen als Fremdgruppe. B Hauptkomponentenanalyse (PCA) Diagramm der ersten beiden Komponenten für alle Beitritte. C STRUKTURanalyse aller Baumwollakzessionen mit unterschiedlicher Clusteranzahl K = 6 und K = 12 (K = 12 ist der optimale Wert). Die x-axis listet die Fremdgruppenarten auf (grau), G. barbadense (Blau), G. hirsutum Landrassenzugänge (orange) und G. hirsutum verbesserte Beitritte (grün) bzw. die ja-axis quantifiziert die genetische Vielfalt in jedem Beitritt. Die weiteren Strukturergebnisse sind in der Zusatzdatei 2 dargestellt: Abbildung S2. D Nukleotiddiversität (π) und Fixationsindexdivergenz (Fst) in den fünf Gruppen. e Die Anzahl der Deletionen, Duplikationen, Inversionen und Translokationen in fünf Populationen (zweiseitiger Wilcoxon-Rangsummentest für benachbarte Gruppen, P < 0,001). Jeder Knoten repräsentiert einen Beitritt. In dieser Analyse wurde die Anzahl der SVs mit dem TM-1-Referenzgenom gezeigt

Wir verwendeten 742 Baumwollakzessionen mit einer hohen Sequenzierungstiefe (> 10×) gegen die G. hirsutum „TM-1“-Referenzgenom (Zusatzdatei 1: Tabelle S1 Zusatzdatei 3) und identifizierte 32.099 Deletionen, 7576 Duplikationen, 1112 Inversionen und 357 Translokationen (Zusatzdatei 1: Tabelle S7). Es gibt mehr SVs in Ghlandrace als in den GhImpUSO- und GhImpCHN-Gruppen (Abb. 1e). Darüber hinaus wurden 173.166 (MAF ≥ 0,01) Kopienzahlvariationen (CNVs) in den 742 Akzessionen identifiziert, darunter 82.431 in den Landrassen, 59.309 in der GhImpUSO-Gruppe und 38.057 in der GhImpCHN-Gruppe (Zusatzdatei 1: Tabelle S8). Populationsgenetische Eigenschaften von CNVs in 742 Akzessionen zeigten, dass G. hirsutum Landrassen wurden deutlich von den verbesserten Akzessionen getrennt, ähnlich dem SNP-basierten Ergebnis, wurden jedoch zusammen mit den Akzessionen GhImpUSO und GhImpCHN geclustert (Zusatzdatei 2: Abbildung S4). Diese Ergebnisse deuten darauf hin, dass CNVs mit hohem Vertrauen eine starke Divergenz zwischen G. hirsutum Landrasse und verbesserte Population und kann verwendet werden, um komplexe quantitative Trait-Loci (QTLs) zu entdecken. Dieser umfassende Variom-Datensatz bietet eine genomische Ressource für Baumwollpopulationsgenetik, Domestikationsanalyse und agronomische Allelidentifizierung (Zusätzliche Datei 2: Abbildung S5).

Beweise für genomische Divergenz während der Domestikation und Verbesserung

Domestikationsbezogene Merkmale ergeben sich aus ausgewählten genetischen Variationen bei Wildarten, die sich auf Samengröße, Blütezeit, Ertrag, Qualität und Pflanzenanpassung auswirken [35,36,37]. Um potenzielle Selektionssignale während der Baumwolldomestikation zu identifizieren, scannten wir genetische Variationen mit Allelfrequenzdifferenzierung in der Nukleotiddiversität, indem wir jede kultivierte Gruppe mit ihrer entsprechenden Wildgruppe verglichen. Wir identifizierten 76 Domestikations-Sweep-Regionen (DSRs) mit πLandrasseVerbessert (Ratio ≥ 15) und eine Likelihood-Methode (XP-CLR, Top 5%) (Zusatzdatei 2: Abbildung S6a), die 66,8 Mb im A-Subgenom und 51,4 Mb im D-Subgenom belegt und mit 837 und 1272 Genen assoziiert ist, einschließlich 274 homologe Genpaare (Abb. 2a). Im Vergleich zu früheren Studien mit kleinen Akzessionszahlen [3,4,5] identifizierte diese Domestikationsselektionsanalyse 31 neue DSRs mit einer Größe von 43,6 Mb (Zusatzdatei 1: Tabelle S9). Einige faserbezogene und bekannte domestizierte Gene wurden zwischen Wild-/Landrassen und verbesserten Sorten unterschiedlich exprimiert (Zusatzdatei 2: Abbildung S6b, c). Die ausgewählten Domestikationsgene waren an der Stressreaktion, der Zellwandregulierung, der Jasmonsäure, dem Ethylen und dem zirkadianen Rhythmus beteiligt (Zusatzdatei 2: Abbildung S7). Eine weitere Manipulation dieser Gene im Pflanzenhormon- und Stressreaktionsweg kann dazu beitragen, ihre mutmaßliche regulatorische Rolle bei der Verbesserung der Faserqualität und der Umweltanpassung während der Baumwolldomestikation zu veranschaulichen [3, 38, 39]. Wir haben auch 120 Mb (πGhImpUSOGhImpCHN ≥ 2) mit Verbesserungssignalen, darunter 1006 selektierte Gene im A-Subgenom und 2369 im D-Subgenom mit 353 homologen Genpaaren (Abb. 2a Zusatzdatei 2: Abbildung S6d) und 79,5% (95,4 Mb) der Verbesserungsselektionsregionen wurden bisher nicht identifiziert [5] (Zusatzdatei 1: Tabelle S10). Bemerkenswert ist die Beobachtung, dass 19 Mb Sequenz sowohl mit Domestikations- als auch Verbesserungsselektionssignalen gescreent wurden, wobei das D-Subgenom (441 Gene) mehr Gene als das A-Subgenom (50 Gene) aufweist (Zusätzliche Datei 1: Tabelle S11). Diese Daten legen nahe, dass das D-Subgenom sowohl bei Domestikations- als auch bei Verbesserungsprozessen stärkere SNP-basierte Selektionssignale aufweist.

Multiskalenvariation für subgenomische Divergenz und GWAS auf agronomische Merkmale während der Baumwolldomestikation. ein Circos-Plot, der die SNP- und SV-basierten Selektionssignale und QTLs während der Baumwolldomestikation und -verbesserung zeigt. Die Selektionsregion wurde in einem 1-Mb-Gleitfenster mit einer Schrittweite von 200 kb berechnet. I–VIII, Circos-Plot von den äußeren zu den Zwischenspuren mit Gendichte (I), snpQTLs (II), cnvQTLs (III), dem Verhältnis der Nukleotiddiversität (π) basierend auf SNPs zwischen 256 Landrassen und 1364 verbesserten Akzessionen für die Domestikation (IV .) ( VII). Die Spur (VIII) repräsentiert das domestizierte Homologe. Obere und untere Felder (VI) repräsentieren den Unterschied des Deletions- bzw. Duplikationsvariations-Allels. Die snpQTLs wurden anhand der Meta-GWAS-Analyse von 890 Baumwollakzessionen identifiziert. Der äußerste Kreis der violetten und gelben Schrift des Circos-Plots zeigt pleiotrope snpQTLs (psnpQTLs) bzw. pleiotrope cnvQTLs (pcnvQTLs). Bi Selektive Signale von Kopienzahlvariationen (CNVs) zwischen den A (B) und D (F) Subgenom während der Domestikation. Die horizontalen grau gestrichelten Linien zeigen die Domestikationssignalschwelle mit dem Verhältnis der Nukleotiddiversität zwischen Wild-/Landrasse und verbesserten Baumwollakzessionen (πLandrasseVerbessert >200). c–e und g–i Sechs CNV-basierte GWAS-Treffer, die sich mit Domestikationsauswahlsignalen überlappten, werden für den Samenindex (SI) angezeigt (C), Faserlänge (FL) (D), Kapselgewicht (BW) (e), Fasergleichmäßigkeit (FU) (g), Faserdehnung (FE) (h) und Blütedatum (FD) (ich). Der Schwellenwert der cnvQTL-Linie war -log10 P = 4.4. Der Violinplot zeigte eine phänotypische Variation mit dem führenden CNV-Genotyp. Die Zahlen im Violinplot zeigen die Anzahl der Akzessionen für jede Kopie. Der Signifikanzunterschied wurde mit dem zweiseitigen Wilcoxon-Rangsummentest berechnet (**P < 0,01, *P < 0,05)

Die Domestikation ist ein Treiber für den CNV-Allel-Frequenzunterschied zwischen Wild-/Landrassen- und domestizierten Gruppen [37]. Insgesamt wurden 286 nicht-redundante CNV-basierte Regionen mit Selektionssignalen während der Baumwolldomestikation identifiziert, bestehend aus 297 Mb im A-Subgenom (Abb. 2b) und 105 Mb im D-Subgenom (Abb. 2f). Ungefähr 55% (65 Mb von 118 Mb) der SNP-basierten Domestikationssignale überlappten CNV-basierte Domestikations-Sweeps (Zusätzliche Datei 1: Tabelle S12). Insgesamt wurden 217 CNV-Regionen mit Verbesserungsselektionssignalen identifiziert, die 156 Mb im A-Subgenom und 133 Mb im D-Subgenom umfassen. Etwa 44% (52 Mb von 120 Mb) der SNP-basierten Verbesserungssignale überlappten die CNV-basierten Verbesserungssignale (Zusätzliche Datei 1: Tabelle S13). Insgesamt identifizierten wir 329 Mb (die 6339 Gene abdecken) von Sequenzen im A-Subgenom und 127 Mb (4955 Gene) im D-Subgenom mit sowohl SNP- als auch CNV-basierten Domestikationssignalen. Insgesamt 173 Mb (5526 Gene) und 184 Mb (8405 Gene) an Sequenzen weisen Verbesserungssignale in den A- und D-Subgenomen auf. Die Identifizierung von Selektionssignalen während der Domestikation und Verbesserung kann die weitere Identifizierung von genetischen Loci wichtiger agronomischer Merkmale erleichtern.

Um QTLs für Selektionssignale im Zusammenhang mit agronomischen Merkmalen zu identifizieren, führten wir eine genomweite Assoziationsstudie (GWAS)-Metaanalyse von 890 . durch G. hirsutum Zugänge aus drei unabhängigen experimentellen Fällen mit mehreren Umgebungen (Zusatzdatei 3) [3, 5, 6]. Unter Verwendung der genotypischen Daten von 2.291.437 hochwertigen SNPs mit MAF ≥ 0,05 in 890 Akzessionen identifizierten wir 2952 signifikante SNPs (0,05/2.291.437 P < 2,18 × 10 – 8 ) assoziiert mit faserqualitätsbezogenen Merkmalen. Nach strenger Filterung wurden 91 wichtige faserbezogene QTLs gefunden, darunter 11 für die Faserlänge (FL), 17 für die Faserdehnung (FE), 15 für die Faserfestigkeit (FS), 19 für die Faserlängengleichmäßigkeit (FU), 10 für die Faser micronaire (FM), 7 für Faserreife (MAT) und 12 für Spinnkonsistenzindex (SCI) (Zusatzdatei 1: Tabelle S14 und Zusatzdatei 2: Abbildung S8). Wir identifizierten auch 31 ertragsbezogene und 3 blütendatumsbezogene (FD) QTLs. Insgesamt wurden 125 Haupt-QTLs mit 4751 Kandidatengenen für 15 agronomische Merkmale identifiziert, von denen 78 mit früheren Studien übereinstimmten [3, 5, 6, 15, 40, 41] und die anderen 47 in Metaanalysen neu entdeckt wurden ( Zusatzdatei 1: Tabelle S14). In den 125 QTLs haben 14 Selektionssignale während der Domestikation und Verbesserung (Zusatzdatei 1: Tabelle S15). Darüber hinaus zeigten einundzwanzig QTL-Loci pleiotrope Effekte auf Faserqualität, Ertrag und Blühdatum (Abb. 2a Zusatzdatei 1: Tabelle S16). Zum Beispiel sind der Flusenanteil (LP), das Fasergewicht pro Boll (FWPB) und der Flusenindex (LI) Komponenten des Ertragsmerkmals, wobei die wichtigsten QTLs auf Chromosom D02 kolokalisiert sind (zusätzliche Datei 2: Abbildung S9a). Die LP, FD und Whole Growth Period (WGP) für Blütezeitmerkmale haben QTLs auf Chromosom D03 (zusätzliche Datei 2: Abbildung S9b).

Wir konzentrierten uns auf neue QTLs im Zusammenhang mit der Faserdehnung, die in meta-GWAS identifiziert wurden. Ein neuer QTL (mqFE253) wurde auf dem D05-Chromosom (bei 11,3–12,5 Mb der genomischen Region) lokalisiert. Die 64 Kandidatengene wurden durch Integration von Haplotypanalyse, Genexpression und funktioneller Annotation vorhergesagt (zusätzliche Datei 2: Abbildung S10). Ein Kandidatengen (Ghir_D05G013680, GhIDD7), das für einen Transkriptionsfaktor der unbestimmten Domäne 7 kodiert, wurde in vier Faserentwicklungsstadien unterschiedlich exprimiert (zusätzliche Datei 2: Abbildung S10f). Akzessionen, die zwei Haupthaplotypen der 5′-UTR-Region repräsentieren, zeigten einen signifikanten Unterschied in der Faserdehnung und Faserlänge (Zusatzdatei 2: Abbildung S11a-b). Nach KO GhIDD7, war die reife Faser signifikant kürzer als bei Wildtyp-Pflanzen (25,8 ± 0,3 vs. 27,1 ± 0,1) (Zusätzliche Datei 2: Abbildung S11c, d, e). Diese Ergebnisse zeigten, dass GhIDD7 war ein zuvor nicht charakterisiertes Gen, das zu einem qualitätsbezogenen Merkmal der Faser beitrug.

GWAS-Analyse von 26.831 CNVs mit hoher Konfidenz (MAF ≥ 0,05) in 419 G. hirsutum Akzessionen ergaben 370 signifikante CNVs für 50 QTLs (cnvQTLs) (Zusatzdatei 1: Tabelle S17), von denen 5 pleiotrope Effekte sowohl auf die Faserqualität als auch auf die Flusenausbeute zeigten (Abb. 2a). Dreizehn cnvQTLs überlappten mit SNP-basierten QTLs (snpQTLs), und die anderen 37 cnvQTLs werden nur von CNVs identifiziert. Von diesen cnvQTLs überlappten 15 mit Domestikations-Sweeps und 10 mit Verbesserungs-Selektionssignalen (zusätzliche Datei 1: Tabelle S18). Die phänotypischen Daten zeigen einen signifikanten Unterschied bei Baumwollakzessionen mit unterschiedlichen Kopienzahlen der Blei-CNV (Abb. 2c–e, g–i Zusatzdatei 2: Abbildung S12). Zum Beispiel wurde eine Assoziation des Samenindex (SI) mit dem Domestikationssignal auf dem A06-Chromosom identifiziert (Fig. 2c). Eine Assoziation der Faserlänge (FL) mit dem Domestikationssignal war auf dem A10-Chromosom lokalisiert, und FL mit 2 Duplikationskopien war signifikant länger als die mit 0 Kopien (Referenz) Allel (P < 0,01) (Fig. 2d). Die führende CNV-involvierte LD-Region hat 78 kodierende Kandidatengene, von denen einige an der Baumwollfaserentwicklung beteiligt sind, wie die UDP-Glucose-Pyrophosphorylase 3 (Ghir_A10G024310, UGP3) und AP2/B3-ähnlicher Transkriptionsfaktor (Ghir_A10G023950). Ein weiteres Beispiel zeigt, dass eine Assoziation der Faserreife (MAT) mit einem Verbesserungsselektionssignal auf dem A12-Chromosom lokalisiert wurde (Zusätzliche Datei 2: Abbildung S13a, b, c). Diese Assoziation hat ein Kandidatengen, das für Xyloglucan-Endotransglucosylase/Hydrolase 5 kodiert (Ghir_A12G008500, XTH5). Im D-Subgenom wurden drei cnvQTLs mit starken Selektionssignalen mit FD, FWPB und FS auf den Chromosomen D03, D06 und D07 assoziiert (Zusatzdatei 2: Abbildung S13d, e, f, g). Diese Ergebnisse liefern eine Reihe von cnvQTL-Kandidaten, die verwendet werden können, um wünschenswerte Merkmale in der zukünftigen Züchtung zu kultivieren.

Pan-Genome von G. hirsutum und G. barbadense Spezies

Wir verwendeten einen referenzgesteuerten Assemblierungsansatz [21], um Pan-Genome von . zu konstruieren G. hirsutum und G. barbadense. Die Sequenzierungsdaten von 1581 G. hirsutum (251 Landrassen, 424 GhImpUSO und 906 GhImpCHN) und 226 G. barbadense verbesserte Akzessionen wurden auf die Referenzgenome „TM-1“ bzw. „3–79“ ausgerichtet [12]. Etwa 5800 Millionen nicht zugeordnete Lesevorgänge von G. hirsutum und 1127 Millionen nicht zugeordnete Lesevorgänge von G. barbadense wurden de novo-Assembly (Zusatzdatei 2: Abbildung S14, S15) unterzogen, wobei 5.047.083.790 bp bzw. 1.517.253.311 bp Contig-Sequenz mit einer Mindestlänge von 500 bp hergestellt wurden (Zusatzdatei 1: Tabelle S19). Nach Entfernen von Redundanzen wurden 3704 Mb und 1422 Mb Nicht-Referenzsequenzen mit einem Contig N50 von 1530 bp (G. hirsutum) und 1108 bp (G. barbadense) alle Filterschritte für die endgültigen Nicht-Referenzgenome bestanden (Zusatzdatei 1: Tabelle S20). Die letzten 1041 Mb und 309 Mb Nicht-Referenzsequenzen in G. hirsutum und G. barbadense mit einer Contig-Länge von mehr als 1000 bp wurden zur Vorhersage proteinkodierender Gene verwendet (Zusatzdatei 2: Abbildung S16). Wir haben 32.569 . erhalten G. hirsutum Gene (65.679 Transkripte) und 8851 G. barbadense Gene (12.076 Transkripte) (Zusätzliche Datei 1: Tabellen S21-S22). Der endgültige G. hirsutum pan-Genom (Ghpan-Genom) ist 3388 Mb mit 102.768 Genen (2347 Mb mit 70.199 Genen im Referenzgenom „TM-1“) und G. barbadense (Gbpan-Genom) ist 2575 Mb mit 80.148 Genen (2266 Mb mit 71.297 Genen im Referenzgenom „3–79“) (Zusatzdatei 2: Abbildung S17).

Die Abdeckung des Ghpan-Genoms wurde mit PacBio-Reads von 10 repräsentativen Akzessionen untersucht, darunter G. hirsutum yucatanense, G. hirsutum richmondi, G. hirsutum morrilli aus den Wild-/Landrassen, Acala, Paymaster 54, Stoneville 2B von der GhImpUSO-Gruppe und Simian 3, CRI 7, Xinluzao 42 und Xuzhou 142 von der GhImpCHN-Gruppe (Zusatzdatei 1: S23-S25 Zusatzdatei 2: Abbildung S18 ). Nach der De-novo-Assemblierung (Zusatzdatei 3) wurden mehr als 93% der assemblierten Contigs dem TM-1-Referenzgenom zugeordnet. Ungefähr 18,9 Mb nicht kartierter Contigs (insgesamt 641 Mb Contigs von 10 Akzessionen, die nicht auf dem TM-1-Referenzgenom kartiert wurden) wurden mit den Nicht-Referenzsequenzen von 1581 . abgeglichen G. hirsutum Akzessionen (die durchschnittliche Länge der Nicht-Referenzsequenz beträgt

655 KB 1041 MB/1581 MB). Die PacBio-basierten Assemblies liefern Beweise für Nicht-Referenz-Genomsequenzen in G. hirsutum, was darauf hindeutet, dass unsere Pipeline der Pan-Genom-Konstruktion PAVs in einer großen Keimplasma-Population abrufen kann. Einige hochfrequente PAVs wurden auch in 23 repräsentativen Akzessionen mittels PCR verifiziert (Zusatzdatei 2: Abbildung S19).

Für die G. hirsutum Population, kartierten wir Resequenzierungs-Reads gegen 102.768 pan-Gene, was zu 17.100 Genen (16,64%, Singleton) in 561 Akzessionen (Sequenzierungstiefe < 5) und 85.667 Genen in 1020 Akzessionen (Tiefe > 5) führte. Der 1020 G. hirsutum Zu den Akzessionen gehören 63.489 Kerngene, die von allen geteilt werden G. hirsutum Akzessionen, 5941 (5,78 %) Softcore-Gene in 990–1019 Akzessionen (97–100 %), 3803 (3,7 %) Shell-Gene in 11–989 Akzessionen (1–97 %) und 12.434 (12,1 %) Wolken in weniger als 10 Akzessionen (0–1 %) (Abb. 3a, b). Für die G. barbadense Pan-Genom traten die 1536 Singleton-Gene nur in 49 Akzessionen mit geringer Tiefe auf. Wir verwendeten 78.612 pan-Gene, die in 177 Akzessionen auftraten, für die weitere PAV-Analyse. Die 177 G. barbadense Akzessionen umfassen 68.789 (85,8 %) Core-Gene, 1796 (2,24 %) Softcore-Gene in 172–176 Akzessionen (97–100 %), 5867 (7,32 %) Shell-Gene in 4–171 Akzessionen (2–97 %) und 2160 (2,75 %) Wolken in weniger als 3 Akzessionen (0–2 %) (Abb. 3c, d). Die Modellierung der Pan-Genomgröße mit iterativ zufälligen Stichproben legt nahe, dass das Ghpan-Genom durchschnittlich 81.688 pan-Gene und durchschnittlich 65.595 Core-Gene in 398 Akzessionen aufweist (Abb. 3e). Das Gbpan-Genom weist durchschnittlich 78.607 pan-Gene und 69.563 Core-Gene in 59 Akzessionen zur Modellierung der Sättigung auf (Abb. 3f). Daher nahm die Größe des Core-Genoms ab und das Pan-Genom nahm mit der Zunahme der Populationsgröße zu. Die GO-Analyse zeigte, dass Kerngene am zellulären Stoffwechselprozess und der Entwicklung beteiligt sind, während die variablen Gene an der „Abwehrreaktion“, der „Reaktion auf Stress“ und der „Signaltransduktion in der Umgebungsfitness“ beteiligt sind (Zusatzdatei 2: Abbildung S20).

Pan-Genome von G. hirsutum und G. barbadense Spezies. ein Gennummer und Präsenzhäufigkeit in G. hirsutum pan-Gene. Das Tortendiagramm entspricht den Kern- (in allen Beitritten vorhanden), Softcore-, Shell- und Cloud-Genen. Die Singleton-Gene in Akzessionen mit geringer Tiefe (< 5) wurden für eine weitere PAV-Analyse ausgeschlossen. Die variablen Gene sind in Zusatzdatei 2: Abbildung S17 in Referenz- und Nicht-Referenzgene unterteilt. B 1020 G. hirsutum Die Heatmap der Beitritte zeigte das Vorhandensein und Fehlen von variablen PAVs. C Gennummer und Präsenzhäufigkeit in G. barbadense pan-Gene. D 177 G. barbadense Die Heatmap der Beitritte zeigte das Vorhandensein und Fehlen von variablen PAVs. e, f Sättigungskurve, die die Zunahme der Pan-Genomgröße und die Abnahme der Kern-Genomgröße in 1020 . modelliert G. hirsutum (e) und 177 G. barbadense (F). Der Fehlerbalken wurde basierend auf 1000 zufälligen Kombinationen mit fünf Replikaten von Baumwollgenomen berechnet. Die oberen und unteren Ränder in Lila und Rot repräsentieren die maximale und minimale Genzahl. Die durchgezogenen Linien stellen die Anzahl der pan-Gene und der Core-Gene dar

Als nächstes untersuchten wir die genomischen Eigenschaften von Kern- und variablen Genen zwischen A- und D-Subgenom. Core-Gene haben in beiden Fällen höhere Expressionslevel als variable Gene G. hirsutum und G. barbadense (Zusätzliche Datei 2: Abbildung S21). Interessanterweise haben subgenomische variable Gene von A höhere Expressionsniveaus als subgenomische D-Gene ( 4a ). Variable Gene haben eine höhere benachbarte (2 kb) TE-Insertionswahrscheinlichkeit als Core-Gene, insbesondere für die Zigeuner Klasse (Zusatzdatei 2: Abbildung S22). Die variablen Gene im D-Subgenom haben ein höheres Verhältnis als die im A-Subgenom (Abb. 4b). Evolutionäre Selektionsanalysen zeigten, dass in beiden Fällen mehr variable Gene positiv selektiert wurden als Kerngene G. hirsutum und G. barbadense, insbesondere im D-Subgenom (Fig. 4c). Darüber hinaus haben variable Gene eine größere Nukleotid-Diversität als Core-Gene, und variablere Gene im D-Subgenom haben eine höhere Diversität (P < 0,001) (Fig. 4d Zusätzliche Datei 2: Fig. S23). Diese Daten zeigten, dass subgenomische variable D-Gene eine schnellere Evolutionsrate aufwiesen als subgenomische A-Gene.

Vergleich von Kern- und variablen Genen in A- und D-Subgenomen. ein Expressionsniveaus von Kern- und variablen Genen in G. hirsutum und G. barbadense. Die Softcore-Gene werden durch „Soft“ repräsentiert. B Verhältnis der Insertionshäufigkeit des transponierbaren Elements (TE) in den stromaufwärts gelegenen 2 kb der Kern- und der variablen Gene in den A- und D-Subgenomen. C Verhältnis nicht synonym/synonym (Kein/KS) Mutationen von Kern- und variablen Genen. D SNP-Diversität von Kern- und variablen Genen. Der Vergleich von Genexpression, TE- und SNP-Diversität zwischen Kern- und variablen Genen wurde mit einem zweiseitigen Kolmogorov-Smirnov-Test (*P < 0,05, **P < 0,01, ***P < 0,001)

PAV-Auswahl während der Domestikation und Verbesserung

Um eine Landschaft selektiver PAVs zwischen Landrasse und verbesserter Baumwolle zu etablieren, verglichen wir die PAV-Häufigkeit zwischen den Gruppen Landrasse, GhImpUSO und GhImpCHN. Die Landrassengruppe hat mehr variable Gene als verbesserte Sorten, was auf einen allgemeinen Trend zum Genverlust während der Baumwolldomestikation hindeutet (Abb. 5a). PCA und phylogenetische Analyse von PAVs legen nahe, dass die Gruppe der Landrassen von der Gruppe der verbesserten Sorten getrennt wurde (Abb. 5b, c). Die aus den Ureinwohnern Amerikas stammenden Landrassen wiesen in ihrer genetischen Zusammensetzung eine Populationsmischung mit amerikanischer Kulturbaumwolle auf, die mit der Clustering-Analyse von High-Confidenz-SNPs übereinstimmt (Zusatzdatei 2: Abbildung S24). Um die Falsch-Positiv-Rate zu kontrollieren, wurden acht Landrassen und vierunddreißig GhImpUSO-Zugänge in einer gemischten Populationsstruktur mit ungewisser Herkunft von der weiteren Analyse ausgeschlossen.

PAV-Selektionssignale während der Baumwolldomestikation und -verbesserung. ein Gennummer unter den G. hirsutum Landrasse und verbesserte Beitritte. Der Wilcoxon-Rangsummentest (P < 0,001) wurde für die signifikante Statistik verwendet. B PCA-Analyse von 1020 Akzessionen basierend auf Shell-PAVs. C Stammbaum und Populationsstruktur mit maximaler Wahrscheinlichkeit mit unterschiedlicher Anzahl von Clustern (K = 2, 3 und 4) in 1020 G. hirsutum Beitritte mit 3803 Shell-PAVs. Die Bevölkerungsstruktur ist nach dem phylogenetischen Baum sortiert. d, e Vergleich der signifikanten Genpräsenzhäufigkeit zwischen der Landrasse gegenüber der GhImpUSO-Gruppe (Domestikation) und der GhImpUSO gegenüber der GhImpCHN-Gruppe (Verbesserung) (FDR < 0,001, zweiseitiger exakter Fisher-Test). F Anzahl günstiger und ungünstiger Gene während der Domestikation und Verbesserung. g, ha PAV-Präsenzhäufigkeit günstiger und ungünstiger Gene während der Domestikation und Verbesserung. ich, j GO-Anreicherungsanalyse des günstigen Gens (ich) und ungünstiges Gen (J) Gewinn und Verlust während der Domestikation und Verbesserung

Um PAV-bezogene Gene mit Selektionssignalen während der Domestikation und Verbesserung zu identifizieren, führten wir zwei Vergleiche zwischen 182 Landrassen und 206 GhImpUSO-Akzessionen unter Verwendung der Anwesenheitshäufigkeit variabler Gene für die „Domestikation“ durch (Abb. 5d Zusätzliche Datei 2: Abbildung S25) und zwischen 206 GhImpUSO und 592 GhImpCHN-Zugängen zur „Verbesserung“ (Abb. 5e). Die Gene mit einer signifikanten Änderung der Anwesenheitsfrequenz (FDR < 0,001 und Häufigkeitsänderung > 2 für „ungünstiges Gen“ bzw. < 0,5 für „günstiges Gen“) wurden als ausgewählte Gene angesehen. Gene mit einer höheren Vorkommenshäufigkeit in Landrassen als in GhImpUSO und einer höheren Vorkommenshäufigkeit in GhImpUSO als in GhImpCHN waren potenziell „ungünstige Gene“, während Gene mit umgekehrten Mustern der Vorkommenshäufigkeit „günstige Gene“ waren. Wir identifizierten 2785 und 7867 günstige Gene mit Allelgewinn und 6753 und 3866 ungünstige Gene mit Allelverlust während der Domestikation bzw. Verbesserung (Zusatzdatei 1: Tabellen S26, S27). Die GO-Anreicherungsanalyse zeigte, dass günstige Gene im Oxidations-Reduktionsprozess angereichert wurden, während ungünstige Gene in der Fettsäurebiosynthese und Genregulation angereichert wurden. Die günstigen und ungünstigen Gene wurden in vier Vergleiche entsprechend der Anwesenheitshäufigkeit in drei Gruppen während der Domestikation und Verbesserung eingeteilt (Abb. 5f). Die kontinuierliche Selektion von 337 günstigen Genen mit Domestikations- und Verbesserungssignalen kann Elitekandidaten für die Züchtung sein, während 308 ungünstige Gene mit niedrigeren Präsenzfrequenzen in der GhImpCHN-Gruppe Verlustallele darstellen (Abb. 5g Zusatzdatei 1: Tabelle S28). Bei der Baumwollzüchtung wurden mehr ungünstige als günstige Gene eliminiert (Abb. 5h). Gene für günstigen Gewinn waren am Transmembrantransport und am Oxidations-Reduktionsprozess beteiligt, während Gene mit günstigem Verlust an der Elektronentransportkette und dem sekundären Stoffwechselprozess beteiligt sind (Fig. 5i, j). Ungünstige Verstärkungsgene hatten keinen signifikant angereicherten Prozess während der Verbesserung (Fig. 5j). Diese Analysen zeigten, dass viele ungünstige Gene während der Domestikation verloren gingen und beträchtliche günstige Gene während des Verbesserungsprozesses erhalten blieben.

Gene für verwandte Merkmale unter Verwendung eines Pan-Genom-Datensatzes

Basierend auf den obigen Daten schlagen wir ein zusammenfassendes Diagramm für die natürliche Selektion, Domestikation und Verbesserung von Baumwolle vor (Abb. 6a). Über die integrierten SNP-, CNV- und PAV-Karten identifizierten wir fast 456 Mb (19,4% des assemblierten Referenzgenoms) und 357 Mb (15,2%) von Sequenzen mit Domestikations- und Verbesserungssignalen (Zusatzdatei 1: Tabelle S29). Es gibt 21.169 Gene in Domestikationsregionen, von denen einige nachweislich an der Regulierung des Blühdatums, der Morphologie und der Faserentwicklung beteiligt sind. Für das Blütedatum weist ein signifikanter GWAS-Peak auf Chromosom D03 zwei Kandidatengene auf, die für a . kodieren COP1-interaktives Protein [6] (CIPI, Ghir_D03G008950) und ein CONSTANS-ähnliches Protein [42] (COL2, Ghir_D03G011010), die für die Anpassung der Landrassenbaumwolle an Kultursorten in verschiedenen geografischen Gebieten mit unterschiedlichen Photoperioden erforderlich sind. Eine weitere Untersuchung der kausalen SNP-Allele zeigt, dass die Vorfahren-Allele hauptsächlich in Landrassen verbreitet sind, mit geringeren Allelfrequenzen in verbesserten Sorten (Abb. 6b). In ähnlicher Weise fanden wir, dass Landrassen und verbesserte Gruppen Alleldifferenzierung in aufwiesen SPÄTE MERISTEM-IDENTITÄT1 [43] (LMI1, Ghir_D01G021810), das die Blattform reguliert, und im grundlegenden Helix-Loop-Helix-Protein-Gen GRF (Ghir_A12G025340), das ein Kandidatengen für Baumwolldrüsen-QTL ist [44] (Abb. 6b). Einige Gene, die für die Faserentwicklung verantwortlich sind und Domestikation und Verbesserungsselektion erfahren haben, wurden auch durch die geografische Differenzierungsanalyse nachgewiesen. KCS2 (Ghir_D10G015750) und CesA6 (Ghir_D03G004880), die für die Faserdehnung verantwortlich sind [45,46,47,48], einer Domestikations- und Verbesserungsauswahl unterzogen (Abb. 6b). The domestication gene PRF3 (Ghir_D13G021640) has a strongly mutated allele in improved cultivars [49].

An available pan-genome dataset for cotton breeding. ein A four-step model of variation during cotton domestication and breeding. B The spectrum of gene allele frequencies at the causal SNP polymorphisms of COL2, CIP1, PRF3, LMI1, GRF, KCS2, und CesA6 in landrace and two geographic groups. C The spectrum of domesticated PAV allele frequ encies of seven genes in landrace and two geographic groups. D An example of functional PAV located on the A08 chromosome. The dashed line in Manhattan plot indicates the threshold for GWAS signals (P < 2.62 × 10 − 8 −log P > 7.6). This locus includes four QTLs (lint percentage (LP), fiber weight per boll (FWPB), fiber micronaire (FM), fiber strength (FS)). e Four QTLs were displayed in a panel of multiple accessions. The two dashed lines represent GWAS thresholds for CNV (−log P > 6.45) and SNP (−log P > 4.42), respectively. F the phenotypic difference between presence and absence groups. The numbers below the violin plots show the accession numbers. The significance difference was calculated with a two-sided Wilcoxon rank-sum test (***P < 0.001, **P < 0,01). g Presence frequencies of Ghir_A08G006710 in 182 landrace, 206 GhImpUSO, and 592 GhImpCHN accessions

Pan-genome analysis uncovered favorable and unfavorable gene alleles during domestication and improvement, providing novel candidate genes for functional investigation (Fig. 5). For genes favorable to cotton improvement selection, SCD (short chain dehydrogenase, GhirPan.00056999), NS (sugar transporter, GhirPan.00054328), und RbfA (ribosome-binding factor A, GhirPan.00033905) have the lowest frequency in wild population and highest in domesticated cultivars (Fig. 6c Additional file 2: Figure S26). Some favorable genes exhibiting a decrease of frequency in the improvement process could be eliminated (308 genes), having almost the same allele frequency between wild and cultivated accessions, such as DXS (deoxyxylulose-5-phosphate synthase, Ghir_Scaffold1882G000030) und COX3 (cytochrome oxidase subunit 3, Ghir_Scaffold1273G00008). Genes unfavorable during domestication showed increased (182 genes) or decreased (5405 genes) frequency in the GhImpCHN group, such as RLP9 (receptor like protein 9, Ghir_D13G022380) und ZBD (Zinc-binding dehydrogenase, GhirPan.00044196) (Fig. 6c).

To determine the contribution of PAV to agronomic traits, we identified PAV-associated SNPs for 1196 PAVs (MAF ≥ 0.02) in 415 accessions (4 accessions were discarded from 419) using 1,904,926 SNPs and obtained 56,486 significant SNPs (P < 2.62 × 10 − 8 ) associated with 864 (72.2%) PAVs. Of these PAVs, 124 were overlapped with 89 trait-QTLs (Additional file 1: Table S30 Additional file 2: Figure S27). One representative PAV (Ghir_A08G006710, 543 bp, an uncharacterized gene in G. hirsutum) is located on chromosome A08 (Fig. 6d, Additional file 2: Figure S28). This hotspot region contained two yield-related (LP, FWPB) QTLs and two fiber quality-related (FM, FS) QTLs (Fig. 6e). These accessions with the presence haplotype of this gene showed significantly increased appearance of LP and FWPB traits than those with the absence haplotype, but no difference for FS and FM traits (Fig. 6f). Further presence frequency analysis showed that Ghir_A08G006710 was present in nearly all landrace and GhImpUSO accessions, but was absent in only a few GhImpCHN accessions (Fig. 6g). Interestingly, in the population RNA-Seq data of 15 DPA fiber [15], absence of this gene in 18 accessions was accompanied by significant low expression of an adjacent gene Ghir_A08G006730 (locating at upstream

61 kb, encoding an AUX/IAA transcriptional regulator family protein) compared with that representing presence of this gene in 233 accessions, supported by the change of IAA content in fibers of representative accessions (Additional file 2: Figure S29, S30). These results implied that this gene represented a recent loss event with a potential regulatory role in other gene expression during cotton improvement. These PAV localization and QTL analyses may improve the efficiency of identifying favorable genes associated with desirable agronomic traits.


Einführung

Phylogenetic reconstructions have traditionally used only a fraction of the sequence data of an organism’s genome, but due to the widespread application of Next Generation Sequencing (NGS) to phylogenetics the quantity of data continues to increase. Phylogenomic studies have therefore heavily relied on a handful of reduced representation approaches including transcriptome sequencing (RNASeq), DNA-based reduced representation techniques, and genome skimming. RNASeq was among the early, still fairly expensive, techniques to obtain large numbers of loci that are informative for deep phylogenetic divergences. Recently, the more cost-effective sequencing of targeted genomic DNA, enriched via hybrid capture, became popular and is at the core of widely used approaches including Ultra Conserved Element (UCE) (McCormack et al., 2012) and Anchored Hybrid Enrichment (Lemmon, Emme & Lemmon, 2012) methods. As sequencing costs have dropped during the past decade, genome skimming (low coverage whole genome sequencing) has become a viable alternative to target enrichment, at least for taxa with relatively small (1 Gbp) genomes. This technique is less challenging with respect to sample quality, involves less complicated lab protocols and does not require expensive probe synthesis. This last point is critical for sampling phylogenetically diverse taxa because the recovery of target sequences is not bound by limitations of the probe design.

While genome skimming does confer these potential benefits, the resulting data can be difficult to parse or integrate into a phylogenetic dataset and can pose substantial problems for analysis. For example, assembled sequences may differ from deep-sequenced model taxon genomes in being much less contiguous as well as unannotated. Genome skimming data also differ from RNASeq data, most notably by the presence of untranslated highly variable regions such as introns. As opposed to typical target capture data, where targeted loci have much higher coverage than non-target ones (Knyshov, Gordon & Weirauch, 2019), genome skimming produces more uniform coverage across the genome (Zhang et al., 2019), with differences associated primarily with sequence properties such as GC content (Barbitoff et al., 2020). Also unlike hybrid capture methodologies, where probes are typically designed for a particular set of taxa based on a related reference taxon (Faircloth, 2017 Young et al., 2016), genome skimming can be applied to taxa with or without available reference genomes or transcriptomes. Nevertheless, hybrid capture-based bioinformatic solutions are most commonly applied to the phylogenetic analysis of genome skimming data (Chen et al., 2018 Zhang et al., 2019).

Phylogenetically-oriented hybrid capture and genomic pipelines are subdivided into two main groups of approaches. Software in the first group identifies reads of interest with the help of reference sequences and subsequently assembles this limited pool of reads (aTRAM (Allen et al., 2015, 2018), HybPiper (Johnson et al., 2016), Assexon (Yuan et al., 2019), Kollector (Kucuk et al., 2017), and HybPhyloMaker (Fér & Schmickl, 2018)). The search for reads that match target regions typically makes use of read aligners (HybPiper, Kollector, HybPhyloMaker) or local similarity search algorithms on both the nucleotide and protein levels (aTRAM, HybPiper, Assexon). After reads are gathered, they are fed to an assembler, and assembled contigs are further processed. A benefit of this group of approaches is that there is no need to assemble the entire read pool, making them potentially faster and less memory demanding than approaches that use the whole read pool. Some drawbacks are the need to perform new read searches and assemblies for each new set of baits and the inability to work with assembled data.

The second group of approaches uses an assembly compiled from the total read pool. The assembly is queried for target sequences, which are then extracted and processed. Post-assembly dataset-specific target searches can be performed relatively quickly. However, especially for highly divergent taxa, the assembly process itself may be both a memory- and time-demanding procedure. Generating a set of contigs from transcriptomic assemblies can be relatively straightforward, because they mostly consist of spliced protein coding sequences. This approach is utilized in HaMStR (Ebersberger, Strauss & Von Haeseler, 2009), Orthograph (Petersen et al., 2017), Orthofinder (Emms & Kelly, 2019), and FortyTwo (Simion et al., 2017), among other applications. However, unannotated genomic assemblies may have contigs comprised of multiple genes or untranslatable introns of varying size. Gene prediction and protein extraction may be complicated when a target gene is fragmented into many small contigs. Recently, Zhang et al. (2019) suggested using Phyluce (Faircloth, 2016) for UCE extraction and Benchmarking Using Single Copy Orthologs (BUSCO) (Simão et al., 2015 Waterhouse et al., 2017) for OrthoDB Single Copy Ortholog (SCO) extraction from genomes at shallow phylogenetic levels, that is, from relatively closely related taxa. Between these two solutions, only BUSCO is specifically designed for genomic assemblies and has the capability to search for and predict genes de novo, but it is only feasible for a few predetermined sets of proteins. Phyluce was originally designed for short, conserved fragments and it is unclear how well it performs on longer multiexon genes. The recently published Assexon software (Yuan et al., 2019) is capable of searching for and retrieving sequences from genomic assemblies, but this module has not yet been extensively tested.

To address issues with commonly-used techniques for including genome-skimming data in phylogenies, we have developed a software, named ALiBaSeq (ALignment Based Sequence extraction), that is designed for sequence extraction based on a local alignment search and is applicable to all types of assembled data and a wide range of assembly qualities. The software is flexible with respect to both input and output, which will facilitate its incorporation into existing bioinformatics pipelines. Any read processing technique and assembler are supported to generate the input for the software, while the resulting sequences are output in FASTA format and can be grouped in several ways (per target locus, per sample, etc.) depending on what is required in downstream analyses. The software also allows for the integration of different types of datasets (e.g., transcriptomic and sequence capture data) allowing phylogenies with more complete taxon sampling as these various phylogenomic datasets become more and more available (Kieran et al., 2019). One of the software’s particular strengths is its ability to efficiently obtain orthologous regions from unannotated genome skimming data. Existing tools frequently rely on a particular type of sequence aligners (BLAST (Altschul et al., 1990) for aTRAM and FortyTwo, both BLAST and BWA (Li & Durbin, 2009) for HybPiper, Usearch for Assexon, LASTZ (Harris, 2007) for Phyluce). Our software supports several commonly utilized similarity search programs and their outputs. While we provide utility scripts for some of the tools, the aforementioned search programs can be run on their own, thus giving the user full control over search program settings if needed. Finally, compared to other programs, we offer greater customization of parameters, including different alignment score cutoff criteria, specification of number of alternative matches, and sequence output structure. The software is available for download at https://github.com/AlexKnyshov/alibaseq.

We here describe the implementation of this software, assess its performance, and benchmark it against other commonly utilized algorithms. Tests are conducted on (1) both conserved and variable loci as determined by average pairwise sequence distance, on (2) contiguous whole genome assembly, short read assemblies of variable depth of coverage, and a hybrid capture sample. We focus testing on the insect samples (see below), but also perform a subset of tests on a plant system to verify the software’s versatility, the details of which are available in the Text S1. Overall, we find that our software matches or outperforms other techniques applied to genome skimming data in recovering the most orthologous genes with the lowest amount of error in low-coverage, fragmented and unannotated genome assemblies. Furthermore, we determine that it works as well or better than other tools on high coverage genome assemblies and target capture assemblies especially at relatively deep phylogenetic levels (100–200 Mya). Thus, ALiBaSeq is a valuable tool for compilation of phylogenomic datasets across diverse taxa and diverse data types.


Diskussion

Our data demonstrate the complex interaction between heterozygosity, genome assembler, and length thresholding effects with some problems becoming evident only after extensive comparison to a high-quality reference sequence. For example, from the 200 bp size cutoff assemblies, LAST showed an average of 10% sequence missing across the SOAPdenovo2 assemblies when compared to the reference, yet they were an average of 50% larger than the reference, in total assembly size. This suggests regional expansions account for a 60% excess of genomic sequence for these assemblies over the reference (S1 Table). To state this another way, an average of 40% of SOAPdenovo2 assemblies consist of expanded sequence (S1 Table). This may be an underestimate given that some regions have undergone sequence collapse (discussed below) which is also compensated by regional expansion. For the multigene pgp family we showed lower heterozygosity for the SOAPdenovo2 assemblies and one Platanus assembly (Fig 6). We interpret the lower heterozygosity in SOAPdenovo2 assemblies as evidence that these regions are not properly resolved and likely expanded regionally--consistent with duplicate genes observed throughout the phylogenetic tree in isolog clusters (Fig 5).

Confirming this, we performed PANTHER analysis of specific GO categories, yielding highly significant enrichment or depletion of 237 specific categories even after correction for false discovery rate to 0.01 (S5 Table and Fig 7A). These discrepancies can be at least partly explained by a complex interplay between regional heterozygosity and assembly parameters. While the reference genome does not display unusual heterozygosity or coverage of these regions (Fig 10) we documented in four categories that the assemblies of these regions diverge from the reference genome in terms of coverage, heterozygosity, and length assembled (Fig 7B, 7C and 7D). We would predict that if an assembler maximally “spreads out” the variation within a dataset into distinct contigs, length assembled would go up, while coverage and heterozygosity would go down as the reads are able to find their perfect match. In many cases this is precisely what we see: the assemblies shown for Oxidoreductase and Dehydrogenase behave in this way (Fig 7B, 7C and 7D) and are examples of ‘regional expansion’ (Fig 9). Somewhat surprisingly, this regional expansion appears to be far greater than one would expect for separation of alleles, which should lead to a doubling of the sequence length in most cases we saw well over 3-fold expansion of length and in one extreme case 7-fold (Fig 7C). Even Platanus, algorithmically optimized for heterozygous genome assembly, was prone to this artifact under specific parameter settings (Fig 7B and 7C). While Platanus step-size 1 performs particularly poorly with our dataset, step-size 3 and 7 both showed artifacts in our PANTHER analysis (Fig 7, see Oxidoreductase, Dehydrogenase, and Response to Heat) while yielding reasonable N50 values (step-size 3, N50 = 74 kb step-size 7, N50 = 70 kb). Therefore, our data highlight a potentially worrisome problem for genome assembly algorithms when confronted with moderate to highly heterozygous datasets.

The Amino Acid Transport category appears to violate the expectation that heterozygosity will behave similarly to coverage it is increased, not decreased, in two of three SOAPdenovo2 assemblies where coverage was decreased (Fig 7D). Hypothesizing that this might reflect collapsed repetitive elements that are intronically located within these genes, we ran RepeatMasker over the corresponding extracted genomic regions from the reference, SOAPdenovo2 23, 47, and 63, along with Platanus 20 (control). We found that while the reference assembly encodes a highly repetitive component (34.6%), the repetitive content of SOAPdenovo2 23, 47, and 63 were dramatically reduced (4.6%, 9.2%, and 9.2%, respectively). Platanus 20 (control) was 30.3% repetitive. Thus, while the Amino Acid Transport coding regions were expanded in length (Fig 7C) leading to PANTHER enrichment (Fig 7A), these genomic regions encode repeats which are collapsed leading to higher heterozygosity (Fig 7D). Thus, rather than reflecting a simple expansion or contraction (Fig 9), Amino Acid Transport-related genomic regions reflect a combination of expansion and collapse. The reasons for this anomaly remain to be investigated in future work, especially given that the repetitive elements included in these regions are unclassified by RepeatMasker. It is worth noting that the expansion of sequence encoding Amino Acid Transport-related genes, and the collapse of repetitive elements should lead to compensatory changes in coverage and heterozygosity (i.e., increased lengths should decrease the apparent heterozygosity, while collapsed repeats should increase the coverage) but overall deviations from reference are detectable (Fig 7). Indeed, the extreme length extension (7-fold, Fig 7C) of the k-mer 23 assembly may have created the apparent low heterozygosity, offsetting the effect of its highly collapsed repeat (Fig 7D). These data suggest that taken together, coverage and heterozygosity offer better information on genome assembly quality than coverage alone.

The extreme enrichment of heterozygosity for the category ‘response to heat’ for the SOAPdenovo2 23 assembly is particularly striking. While it would suggest the collapse of the genes in this category relative to the reference genome, the expected decrease in sequence length was not observed (Fig 7C). However, to construct Fig 7C we required a 98 percent identical BLASTn match or better between sequences, using blast_analysis.py (Fig 8). By relaxing this requirement to 80 percent identity we found a 3.57-fold contraction (43,083 bp from SOAPdenovo2 23 corresponding with 153,958 bp in the reference genome) which agrees with the 3.49-fold enrichment in heterozygosity (Fig 7D). (Read-mapping was performed with BWA-MEM and does not invoke a percent identity threshold). Platanus step-size 7 represents a curious case: it also is depleted for the ‘response to heat’ category but the increase in heterozygosity was only minor and coverage did not increase, suggesting these regions simply did not assemble well and were likely lost from the assembly when we filtered out contigs smaller than 1 kb, leaving the corresponding reads without a suitable target in the mapping step.


Data availability

Raw sequencing data used in this study can be found in the NCBI database under the following Bioproject accession numbers: PRJNA603155 (genome sequencing dataset of Harukei-3 melon), PRJNA624817 (genome sequencing dataset of seven melon accessions), PRJNA603146 (ONT cDNA RNA-seq), PRJNA603129 (ONT direct RNA-seq), PRJNA603204 (tissue-wide RNA-seq of Harukei-3 melon), or PRJNA603202 (leaf RNA-seq in the greenhouse). Genome assembly and annotation of Harukei-3 melon (ver. 1.41 genome reference) is available on Melonet-DB (https://melonet-db.dna.affrc.go.jp/ap/dnl).


Genom-Annotation

To harness the full potential of a genome sequence, it needs to be annotated with biologically relevant information that can range from gene models and functional information, such as gene ontology (GO) terms (Gene Ontology Consortium 2004 Primmer et al. 2013 ) or ‘Kyoto encyclopedia of genes and genomes’ (KEGG) pathways (Kanehisa and Goto 2000 ), to microRNA and epigenetic modifications (The ENCODE Project Consortium 2012 ). In the context of genetic nonmodel organisms, annotation is often confined to protein-coding sequence (CDS) or transcripts more generally. Despite the considerable challenge to annotate genes in newly sequenced species where preexisting gene models are mostly lacking, automated gene annotation has in principle become possible for individual research groups (Yandell and Ence 2012 ). Still, a complete genome annotation constitutes a considerable effort and requires bioinformatic proficiency. We describe only the general workflow and refer the interested reader to a comprehensive review by Yandell and Ence ( 2012 ) for more details (Box 2). Before starting, it should be noted that successful annotation strongly depends on the quality of the genome assembly. Only contiguous near-complete (

90%) genomes interrupted only by small gaps will yield satisfying results. As a rule of thumb, large genomes have longer genes and thus need more contiguous assemblies for successful annotation (cf. Figure 1 in Yandell and Ence 2012 ).

The annotation process can be conceptually divided into two phases: a ‘computational phase’ where several lines of evidence from other genomes or from species-specific transcriptome data are used in parallel to create initial gene and transcript predictions. In a second ‘annotation phase’, all (sometimes contradicting) information is then synthesized into a gene annotation, following a set of rules determined by the annotation pipeline.

Prior to gene prediction, it is of vital importance to mask repetitive sequences including low-complexity regions and transposable elements. As repeats are often poorly conserved across species, it is advisable to create a species-specific repeat library using tools like RepeatModeler or RepeatExplorer (Novák et al. 2013 ). Once repeats are masked (e.g. with RepeatMasker http://www.repeatmasker.org), ab initio algorithms trained on gene models from related species can be used for baseline prediction of coding sequence (CDS) (e.g. AUGUSTUS Stanke et al. 2006 ). Protein alignments (using e.g. tblastx) and syntenic protein lift-overs from a variety of other species provide a valuable resource to complement the predicted gene models. Arguably, the best evidence comes from detailed EST or RNA-seq data, which in addition to CDS, provides gene models with information on splice sites, transcription start sites and untranslated regions (UTRs). If possible, mRNA should be sequenced strand-specifically, as this helps resolve gene models, facilitates transcriptome assembly and eventually aids in the evaluation of the genome assembly.

In a next step, all the evidence from ab initio prediction and protein-, EST- or RNA-alignments need to be synthesized into a final set of gene annotations. As the evidence is mostly incomplete and sometimes contradicting, this is a difficult task that often benefits from manual curation. Still, several automated annotation tools like MAKER (Cantarel et al. 2008 ) or PASA (Haas et al. 2003 ) exist that incorporate, and weigh the evidence from, several sources. Although these tools generally provide good results, qualitative validation is important (e.g. by assessing the length of open-reading frames). Visual inspection of the annotation is another vital component to detect systematic issues such as intron leakage (introns being annotated as exons due to the presence of pre-mRNA) or gene fusion. Tools like WebApollo (Lee et al. 2013 ) from the GMOD project are particularly useful, as they allow the user to edit the annotation directly through the visual interface.

Publishing the genome

Draft genome sequences are now being produced at an ever-increasing rate. Traditional databases such as ENSEMBL from the European Molecular Biology Labs (EMBL) and the Wellcome Trust Sanger Institute, or genomic databases from the National Center for Biotechnology Information (NCBI) providing access to genomes and meta-information can no longer annotate and curate all incoming genomes. NCBI therefore already provides the possibility to upload draft genome sequences and user-generated annotation. To allow other users to improve the assembly and its annotation, all available raw data should be uploaded, together with the assembled genome and all relevant meta-data, for example as a BioProject on NCBI.


Computational analysis of next generation sequencing data and its applications in clinical oncology

Rucha M. Wadapurkar , Renu Vyas , in Informatics in Medicine Unlocked , 2018

1.5.2 Aligning sequences

After assessing the quality of NGS reads, the reads are aligned to the reference genome . For that UCSC (University of Santa Cruz) and GRC (Genome Reference Consortium) are mainly used as sources of human reference genome [ 59–61 ]. There are some issues in selecting alignment software, the first is solving the problem of ambiguity in mapping short reads to the reference genome, which can be solved by considering paired-end reads as a better option [ 62 ]. Secondly, mutations generated from reads with many mismatches have to be discarded from further analysis steps.


This work was supported by the Netherlands Organization for Scientific Research [Vidi grant 864.14.004] to [B.E.D.] and the Conselho Nacional de Desenvolvimento Científico e Tecnológico [Science Without Borders program] to [D.D.C.] and [F.H.C.].

F. A. Bastiaan von Meijenfeldt and Ksenia Arkhipova contributed equally to this work.

Affiliations

Theoretical Biology and Bioinformatics, Science for Life, Utrecht University, Utrecht, The Netherlands

F. A. Bastiaan von Meijenfeldt, Ksenia Arkhipova, Diego D. Cambuy & Bas E. Dutilh

Centre for Molecular and Biomolecular Informatics, Radboud University Medical Centre, Nijmegen, The Netherlands

Felipe H. Coutinho & Bas E. Dutilh

Instituto de Biologia, Universidade Federal do Rio de Janeiro, Rio de Janeiro, RJ, Brazil

Present Address: Evolutionary Genomics Group, Departamento de Produccíon Vegetal y Microbiología, Universidad Miguel Hernández, Campus San Juan, San Juan, 03550, Alicante, Spain


Schau das Video: Elektriker UV video: Relæteknik - princippet bag sekvensstyring (Juli 2022).


Bemerkungen:

  1. Ridley

    Abhängig von der Art der Arbeit

  2. Ludwig

    Ich denke, dass Sie sich irren. Lassen Sie uns darüber diskutieren. Schreib mir per PN, wir kommunizieren.

  3. Porfirio

    Sie haben den Punkt getroffen. Da ist was dran und ich halte das für eine gute Idee. Ich stimme mit Ihnen ein.

  4. Abdulla

    die Bossy-Sicht, amüsant ...

  5. Mingan

    Sie liegen falsch. Ich bin sicher. Schreib mir per PN, diskutiere es.

  6. Charybdis

    Ehrlich gesagt sind die Kommentare hier viel unterhaltsamer als die Nachrichten selbst. (Natürlich keine Beleidigung für den Autor :))

  7. Styles

    Ich entschuldige mich, aber ich denke, Sie liegen falsch. Schreiben Sie mir in PM.



Eine Nachricht schreiben