Informatie

Frequentie van nucleotidesequenties van specifieke lengte in DNA


  1. Hoe vaak zou je de nucleotidesequentie GGATATCCCC (5' naar 3' richting) bij toeval vinden in een DNA-molecuul?
  2. Hoe vaak verwacht u gemiddeld een specifieke sequentie van 20 nucleotiden te vinden in een genoom met een totale grootte van 4 x 109 basenparen?

Aangezien die kans om elk basenpaar in zijn gegeven positie te vinden 1 op 4 is, zou de totale kans om die bepaalde reeks te vinden dus zijn $0.25^{10}$

Maar hoe pak ik het tweede probleem aan?

TIA


DNA-sequentiestatistieken (2)¶

In het hoofdstuk over het installeren van R heb je geleerd over variabelen in R, zoals scalairen, vectoren en lijsten. U hebt ook geleerd hoe u functies kunt gebruiken om bewerkingen op variabelen uit te voeren, bijvoorbeeld door de functie log10() te gebruiken om de log naar de basis 10 van een scalaire variabele te berekenen x, of gebruik de functie mean() om het gemiddelde van de waarden in een vectorvariabele te berekenen mijnvector:

Je hebt ook geleerd dat je een element van een vector kunt extraheren door de vectornaam te typen met de index van dat element tussen vierkante haken. Om bijvoorbeeld de waarde van het 3e element in de vector te krijgen mijnvector, wij typen:

Een handige functie in R is de functie seq(), die kan worden gebruikt om een ​​reeks getallen te maken die van een bepaald getal naar een ander bepaald getal lopen. Als we bijvoorbeeld de reeks getallen van 1-100 willen maken in stappen van 1 (dwz 1, 2, 3, 4, .97, 98, 99, 100), kunnen we typen:

We kunnen de stapgrootte wijzigen door de waarde van het argument “by” dat aan de functie seq() wordt gegeven, te wijzigen. Als we bijvoorbeeld een reeks getallen van 1-100 willen maken in stappen van 2 (dwz 1, 3, 5, 7, .97, 99), kunnen we typen:

In R is het, net als in programmeertalen zoals Python, mogelijk om a . te schrijven for loop om dezelfde opdracht meerdere keren uit te voeren. Als we bijvoorbeeld het kwadraat van elk getal tussen 1 en 10 willen afdrukken, kunnen we de volgende for-lus schrijven:

In de for loop hierboven, de variabele l is een teller voor het aantal cycli door de lus. In de eerste cyclus door de lus, de waarde van l is 1, en zo l * l =1 wordt afgedrukt. In de tweede cyclus door de lus, de waarde van l is 2, en dus l * l =4 wordt afgedrukt. In de derde cyclus door de lus, de waarde van l is 3, en dus l * l =9 wordt afgedrukt. De lus gaat door totdat de waarde van l is 10. In de tiende cyclus door de lus is de waarde van l is 10, en dus l * l =100 wordt afgedrukt.

Merk op dat de commando's die moeten worden uitgevoerd bij elke cyclus van de for loop moet tussen accolades staan ​​(“<” en “>”).

Je kunt ook een geven for loop een vector van getallen die de waarden bevat waarvan u de teller wilt hebben l volgende cycli in te nemen. U kunt bijvoorbeeld een vector avector met de getallen 2, 9, 100 en 133, en schrijf a for loop om het kwadraat van elk getal in vector af te drukken avector:

Hoe kunnen we een gebruiken for loop om het kwadraat van elk tweede getal tussen 1 en 10 af te drukken? Het antwoord is om de functie seq() te gebruiken om de for loop om elk tweede getal tussen 1 en 10 te nemen:

In de eerste cyclus van deze lus is de waarde van l is 1, en zo l * l =1 wordt afgedrukt. In de tweede cyclus door de lus, de waarde van l is 3, en dus l * l =9 wordt afgedrukt. De lus gaat door totdat de waarde van l is 9. In de vijfde cyclus door de lus is de waarde van l is 9, en zo l * l =81 wordt afgedrukt.

R maakt de productie van een verscheidenheid aan plots mogelijk, waaronder scatterplots, histogrammen, piecharts en boxplots. Als je bijvoorbeeld twee vectoren van getallen hebt mijnvector1 en mijnvector2, kunt u een spreidingsdiagram plotten van de waarden in mijnvector1 tegen de waarden in mijnvector2 met behulp van de functie plot(). Als u de assen op de plot wilt labelen, kunt u dit doen door de functiewaarden plot() te geven voor de optionele argumenten xlab en ylab:

Als je naar de help-pagina voor de functie plot() kijkt, zul je zien dat er veel optionele argumenten (invoer) zijn die daarvoor nodig zijn. Een optioneel argument is bijvoorbeeld de type argument, dat het type plot bepaalt. Standaard tekent plot() een punt op elk gegevenspunt, maar als we instellen type “b” is, dan zal het ook een lijn trekken tussen elk volgend gegevenspunt:

We hebben ingebouwde R-functies gebruikt, zoals mean(), length(), print(), plot(), enz. We kunnen ook onze eigen functies in R maken om berekeningen uit te voeren die u heel vaak wilt uitvoeren op verschillende invoergegevenssets. We kunnen bijvoorbeeld een functie maken om de waarde van 20 plus het kwadraat van een invoergetal te berekenen:

Deze functie berekent het kwadraat van een getal (x), en voeg vervolgens 20 toe aan die waarde. De instructie return() retourneert de berekende waarde. Nadat u deze functie hebt ingetypt, is de functie beschikbaar voor gebruik. We kunnen de functie bijvoorbeeld gebruiken voor verschillende invoernummers (bijv. 10, 25):

U kunt de code zien waaruit een functie bestaat door de naam te typen (zonder haakjes). We kunnen dit bijvoorbeeld proberen door “myfunction'8221 te typen:

Wanneer u R typt, kunt u, als u dat wilt, opmerkingen schrijven door de opmerkingstekst achter het teken “#” te schrijven. Dit kan handig zijn als u enkele R-commando's wilt schrijven die andere mensen moeten lezen en begrijpen. R zal de opmerkingen negeren wanneer het de opdrachten uitvoert. U kunt bijvoorbeeld een opmerking schrijven om uit te leggen wat de functie log10() doet:


Invoering

Studies van genomen op basis van taalkundige benaderingen dateren van enkele decennia terug (Brendel et al. 1986 Pevzner et al. 1989 Searls 1992 Botstein en Cherry 1997 Gimona 2006 Faltýnek et al. 2019 Ji 2020). Een samenspel met methoden van statistische fysica en theorie van complexe systemen bracht nieuwe inzichten in de biologie (Dehmer et al. 2009 Qian 2013). Studies variëren van poging tot op n-gram gebaseerde classificatie van genomen (Tomović et al. 2006 Huang en Yu 2016) tot algoritmen voor optimale segmentatie van RNA's in voorspellingen van secundaire structuren (Licon et al. 2010) en analyse van substitutiesnelheden van coderende genen tijdens evolutie (Lin et al. 2019), om er maar een paar te noemen. Verschillende soorten sequenties in genomen zijn gerelateerd aan meerdere genetische codes (Trifonov et al. 2012) en kunnen zowel vanuit kwantitatief taalkundig oogpunt worden bestudeerd (Ferrer-i-Cancho et al. 2013 Ferrer-i-Cancho et al. 2014) en vanuit een breder perspectief, binnen meer abstracte benaderingen (Neuman en Nave 2008 Barbieri 2012). Onlangs zijn neurale netwerken en deep learning-algoritmen naar voren gekomen als nieuwe hulpmiddelen om nucleotidesequenties te analyseren (Fang et al. 2019 Singh et al. 2019 Melkus et al. 2020 Ren et al. 2020) die bredere vooruitzichten bieden voor studies van genomen. Virussen, die balanceren op de vage grens tussen niet-levend en levend, en dus op de rand van het leven blijven (Villarreal 2004 Kolb 2007 Carsetti 2020), behoren tot de meest interessante onderwerpen van onderzoek.

Het doel van deze brief is om aandacht te vragen voor eenvoudige behandelingen van nucleotidesequenties in virale RNA's door middel van nieuwe parameters, die direct uit genoomgegevens kunnen worden geëxtraheerd. We verwachten dat dergelijke parameters mogelijk kunnen worden gebruikt als een hulpmiddel bij de classificatie van virussen (zie in het bijzonder Wang 2013). Het idee van deze studie houdt verband met de recente uitbraak van COVID-19 en de analyse begon met het vergelijken van menselijke coronavirussen (Su et al. 2016 Wu et al. 2020) en enkele andere virussen. Om relatieve homogeniteit van het materiaal te bereiken, beperken we ons monster alleen tot enkelstrengs RNA-virussen. Zowel positieve- als negatieve-sense-RNA's worden beschouwd. Voor toekomstig gebruik nemen we ook twee retrovirussen op, HIV-1 en HIV-2.

Het blad is als volgt ingedeeld. Samenvatting van gegevens en beschrijving van methoden worden gegeven in de sectie "Gegevens en methoden". De resultaten worden gepresenteerd in de sectie "Resultaten". Ten slotte wordt een korte discussie gegeven in de sectie "Discussie".


3. Resultaten en discussie

3.1. Palindroom optreden over de levensboom

We telden het voorkomen van de 16 palindroomwoorden van lengte 4 (Tabel 1), samen met een gelijk aantal niet-palindroomwoorden van lengte 4 (Tabel 1), in DNA-sequenties van geselecteerde genomen. Twintig verschillende soorten werden geanalyseerd voor elk van de 10 verschillende fylogenetische groepen, d.w.z. de gewervelde dieren, ongewervelde dieren, schimmels, planten, protozoa, mitochondriën, bacteriën, archaea, dubbelstrengs DNA-virussen en retrovirussen. Perfecte herhalingen werden uit de invoersequenties verwijderd om de introductie van een triviale vooringenomenheid uit regio's met een extreem lage complexiteit, zoals telomere of centromere herhalingen, te voorkomen. Voor elke ingevoerde DNA-sequentie en elk 4-meer-woord hebben we vervolgens de verhouding berekend: R van het daadwerkelijke voorkomen van het woord gedeeld door het verwachte aantal keren dat het woord voorkomt, gezien het GC-gehalte en dat van de ingevoerde DNA-sequentie. De meeste palindromen waren ondervertegenwoordigd (R < 1) over alle geanalyseerde genomen. Over het algemeen vertoonden de palindromen een gemiddelde R van 0,86, in tegenstelling tot een gemiddelde R van 1,08 voor de niet-palindroomcontroles (tabel 1). De ondervertegenwoordiging van palindromen was het meest uitgesproken in het genoom van gewervelde dieren, planten, dubbelstrengs DNA-virussen en retrovirussen (Fig. 1). In tegenstelling tot eerdere rapporten waren 20 palindromen zelfs in mitochondriale genomen ondervertegenwoordigd, wat aantoont dat de zeldzaamheid van palindromen in prokaryote genomen niet alleen kan worden verklaard door de selectieve druk die wordt uitgeoefend door restrictie-enzymen. Bijkomende selectieve krachten tegen palindromen kunnen hun impact op de DNA-structuur of hun rol als transcriptiefactor-bindingsplaatsen omvatten. 17 Wat de onderliggende krachten ook zijn, korte palindromen zijn ondervertegenwoordigd in allerlei genomen (Fig. 1). Welke palindromen precies en hoe sterk hangt af van de bron van het DNA. Interessant is dat de intergenoomfrequenties van korte palindromen meer dan tweemaal de variantie van de niet-palindroomcontrolesequenties vertonen (22 versus 9% tabel 1), terwijl de intra-genoomfrequenties, b.v. tussen verschillende chromosomen van hetzelfde organisme, zijn uniform (Fig. 2-4). Hierdoor zijn korte palindromen optimaal geschikt voor het typeren van DNA.

Frequentie van palindromen in een diverse selectie van genomen. Palindroomfrequentie wordt uitgedrukt als de verhouding (R) van optreden tot verwachting. Palindromen zijn ondervertegenwoordigd (R < 1, stippellijn) in alle soorten genomen, het sterkst in gewervelde dieren, planten en virussen, en ze vertonen ongeveer tweemaal de variantie tussen soorten in frequentie (foutbalken) dan niet-palindromen. Per groep werden twintig verschillende genomen geanalyseerd (zie paragraaf 2).

Frequentie van palindromen in een diverse selectie van genomen. Palindroomfrequentie wordt uitgedrukt als de verhouding (R) van optreden tot verwachting. Palindromen zijn ondervertegenwoordigd (R < 1, stippellijn) in alle soorten genomen, het sterkst in gewervelde dieren, planten en virussen, en ze vertonen ongeveer tweemaal de variantie tussen soorten in frequentie (foutbalken) dan niet-palindromen. Per groep werden twintig verschillende genomen geanalyseerd (zie paragraaf 2).

Voorbeelden van palindroom frequentiepatronen. Frequentie van de 16 palindromen van lengte 4 in geselecteerde genomen, uitgedrukt als log2 van verhouding (R) van werkelijk tot verwacht optreden. Hiërarchische clustering werd uitgevoerd op basis van de afstand tussen steden. 23 (Boven) Gemiddelde en variantie door palindroom. (Onder) Ter vergelijking worden de signalen van drie willekeurige reeksen getoond.

Voorbeelden van palindroom frequentiepatronen. Frequentie van de 16 palindromen van lengte 4 in geselecteerde genomen, uitgedrukt als log2 van verhouding (R) van werkelijk tot verwacht optreden. Hiërarchische clustering werd uitgevoerd op basis van de afstand tussen steden. 23 (Boven) Gemiddelde en variantie door palindroom. (Onder) Ter vergelijking worden de signalen van drie willekeurige reeksen getoond.

Variantie van palindroomfrequenties in willekeurige DNA-sequenties van verschillende lengtes (N = 20 voor elk). De gemiddelde variantie voor elk palindroom van lengte 4 over de 20 verschillende sequenties wordt vergeleken met die over de eerste 20 menselijke chromosomen (grijze stippellijn) en over de 20 verschillende chromosomen van gewervelde dieren die zijn geanalyseerd in Fig. 1 (zie aanvullende tabel S1).

Variantie van palindroomfrequenties in willekeurige DNA-sequenties van verschillende lengtes (N = 20 voor elk). De gemiddelde variantie voor elk palindroom van lengte 4 over de 20 verschillende sequenties wordt vergeleken met die over de eerste 20 menselijke chromosomen (grijze stippellijn) en over de 20 verschillende chromosomen van gewervelde dieren die zijn geanalyseerd in Fig. 1 (zie aanvullende tabel S1).

Casestudy's over Caenorhabditis spp. (A), zoogdierchromosomen (B), en sensu stricto gisten (C). De meeste chromosomen worden correct opgelost door clustering op basis van palindroomfrequentie. Perfecte tandemherhalingen werden voorafgaand aan de analyse verwijderd om triviale verschillen van repetitieve regio's te voorkomen. Let op het opvallende verschil tussen DNA van gewervelde dieren en ongewervelde dieren.

Casestudy's over Caenorhabditis spp. (A), zoogdierchromosomen (B), en sensu stricto gisten (C). De meeste chromosomen worden correct opgelost door clustering op basis van palindroomfrequentie. Perfecte tandemherhalingen werden voorafgaand aan de analyse verwijderd om triviale verschillen van repetitieve regio's te voorkomen. Let op het opvallende verschil tussen DNA van gewervelde dieren en ongewervelde dieren.

3.2. Clustering van DNA op basis van palindroomfrequentie

Hier stellen we een gegeven DNA-sequentie voor door een vector van 16 getallen: voor elk van de 16 palindromen van lengte 4, de log2 van de verhouding R van werkelijke tot verwachte frequentie (gezien het GC-gehalte van het geanalyseerde DNA en dat van het palindroom). Wanneer dergelijke vectoren, gegenereerd uit een diverse selectie van DNA-sequenties, werden uitgelijnd en hiërarchisch geclusterd op basis van de stadsblokafstand, werden verschillende DNA-sequenties van dezelfde soort gemakkelijk gegroepeerd (zie Fig. 2 voor een representatieve set van verschillende genomen). De clustering werkte voor alle soorten geteste genoomsequenties - eukaryoot, prokaryoot, plastide of virus - maar de topologie van de resulterende boom was niet fylogenetisch zinvol (figuur 2). Het ontbreken van een grootschalig fylogenetisch signaal bleek evenzeer uit de analyse van de volledige set van 200 genomen (aanvullende figuur S1). De resolutie van palindroomfrequentieclustering zou verder toenemen door gebruik te maken van de 64 verschillende palindromen met lengte 6. Dit zou echter ook vereisen dat de invoersequenties langer zijn. Op basis van de willekeurige sequenties die zijn opgenomen in Fig. 2, bleek de huidige benadering te werken voor sequenties langer dan ongeveer 10 kb. Om een ​​betere schatting te krijgen van de minimaal vereiste grootte van input-DNA, analyseerden we willekeurig gegenereerde sequenties van toenemende lengte (Fig. 3). Boven 9 kb is de gemiddelde variantie van R per palindroom daalde tot onder de waarde die werd verkregen voor verschillende chromosomen van gewervelde dieren (0,025, grijze stippellijn in Fig. 3). Ter vergelijking: de gemiddelde variantie van R per palindroom over menselijke chromosomen was 0,0008 (grijze stippellijn in Fig. 3), wat opnieuw aantoont dat de variantie van de palindroomfrequentie veel lager is binnen het genoom dan tussen het genoom.

Ongewervelde dieren die de kleinste intergenoomvariantie van palindroomfrequentie vertonen (Fig. 1), kozen we: Caenorhabditis soort om zijn discriminerende kracht uit te dagen. De volledige nucleaire genomen van C. briggsae en C. elegans werden vergeleken zoals hierboven beschreven en alle chromosomen werden correct opgelost ondanks de zwakke patronen (Fig. 4A). Clustering op basis van palindroomfrequentie scheidde ook verschillende zoogdierchromosomen die, in tegenstelling tot ongewerveld DNA, het karakteristieke patroon vertoonden dat werd veroorzaakt door een sterke ondervertegenwoordiging van palindromen die een CG-dinucleotide bevatten (ACGT, TCGA, CCGG, GCGC en CGCG Fig. 4B). Dit is in overeenstemming met het model dat bij gewervelde dieren DNA-methylering beperkt is tot cytosines gevolgd door guanine (CpG), terwijl bij ongewervelde dieren cytosines in een bredere context worden gemethyleerd. 25 Spontane mutatie van het palindromische CG naar het niet-palindromische TG door deaminering van gemethyleerd cytosine elimineert dus korte palindromen uit het DNA van gewervelde dieren. De resolutielimiet van palindroomfrequentieclustering werd bereikt met een dataset van zeer vergelijkbare sensu stricto gisten. 26 De verschillende chromosomen van de nauw verwante soorten Saccharomyces cerevisiae, S. bayanus, S. mikatae, en S. kudriavzevi scheidde niet perfect die van de meer in de verte verwante S. castellii deed (Fig. 4C).

Clustering op basis van palindroomfrequentie werkte ook voor prokaryoten, en genereerde soortspecifieke patronen voor zowel archaea als bacteriën. Prokaryote genomen vertoonden zeer diverse patronen (aanvullende figuur S1). Natuurlijke plasmiden van Escherichia coli duidelijk geclusterd met het gastheer-DNA (Fig. 5A). Hetzelfde gold voor bepaalde dsDNA-bacteriofagen zoals Lambda of P2. Andere dsDNA-fagen zoals T3, evenals alle geanalyseerde ssDNA-fagen, vertoonden echter niet dezelfde palindroomfrequentiepatronen als E coli (Fig. 5A). Een interessant beeld kwam naar voren bij het vergelijken van menselijke virussen: terwijl alle ssRNA minus-strengs virussen en de retro-transcriberende HIV geclusterd met menselijk DNA, deden dsDNA-virussen en ssRNA plus-strengs virussen dat niet (Fig. 5B).

Palindroom frequentiepatronen van genomisch DNA van de gastheer (A, E coli B, Homo sapiens zwart gelabeld) en bijbehorende virussen (kleurgecodeerd volgens nucleïnezuurtype van het genoom) of plasmiden (grijs).

Palindroom frequentiepatronen van genomisch DNA van de gastheer (A, E coli B, Homo sapiens zwart gelabeld) en bijbehorende virussen (kleurgecodeerd volgens nucleïnezuurtype van het genoom) of plasmiden (grijs).

3.3. Mogelijke toepassing op metagenomica

Het zich snel ontwikkelende veld van jachtgeweersequencing in de omgeving maakt metagenomische analyses mogelijk van gemeenschappen van micro-organismen, waarvan de meeste niet in het laboratorium kunnen worden gekweekt en daarom tot voor kort onopgemerkt zijn gebleven. 27 Een belangrijke uitdaging bij het interpreteren van gegevens over de sequentiebepaling van jachtgeweren in de omgeving is het samenvoegen van niet-overlappende DNA-structuren in groepen die idealiter overeenkomen met de verschillende soorten aanwezige micro-organismen. 28 Standaardmethoden zoals zoeken naar overeenkomsten met bekende genomen of fylogenetische analyse van merkergenen zijn van beperkt nut bij DNA-fragmenten die zijn bemonsterd van voorheen onbeschreven soorten. 28 Di-, tri- en tetra-nucleotidefrequenties zijn voorgesteld om DNA-handtekeningen te verschaffen. 29–31 Palindroomfrequenties die een soortspecifiek signaal dragen (Fig. 2 en 4), de verhoudingen van voorkomen tot verwachting zoals hier toegepast, kunnen ook nuttig zijn voor het opslaan van omgevingsgegevens over shotgun-sequencing, op voorwaarde dat de te analyseren contigs langer zijn dan 9 kb (Afb. 3). Vanaf de 2007 Tovenaar II Global Ocean Sampling Expedition, die op dat moment voornamelijk nieuwe sequenties produceerde, 32 de honderd grootste contigs, met een grootte tussen 11 en 59 kb, werden geanalyseerd zoals hierboven beschreven. Dit onthulde een divers beeld van palindroomfrequentiepatronen met verschillende grote clusters (aanvullende figuur S2). De geanalyseerde sequenties leverden echter nog steeds geen treffers van hoge kwaliteit op wanneer ze doorzochten met blastn 33 tegen de niet-redundante nucleotideverzameling van de NCBI, met slechts één uitzondering van 99% identiteit met Prochlorococcus faag P-SSM4 (GenBank toegangsnummer AY940168). Het was dus niet mogelijk om het voordeel van palindroomfrequentieclustering met deze dataset te beoordelen. Om toch het potentieel van de methode te testen, hebben we willekeurig 10 niet-overlappende fragmenten met een lengte van 10 kb geselecteerd uit elk van de 20 verschillende bacteriële genomen die zijn geanalyseerd in Fig. 1 (aanvullende tabel S1). Toen deze 200 sequenties werden geclusterd volgens palindroomfrequentiepatronen, werd meer dan 90% correct geassembleerd volgens soort van herkomst.


Resultaten

HT-SELEX experimentele gegevens bieden nauwkeurige m-woordscores voor diverse TF-families

We analyseerden HT-SELEX-gegevens, waaronder 548 experimenten met 410 menselijke en muizeneiwitten uit 40 verschillende TF-families, om m-woord bindende scores. Door de grotere sequentiëringsdiepte konden we langer nauwkeurige scores afleiden m-woorden. Dit aspect is vooral belangrijk omdat de DNA-vorm wordt beïnvloed door de flankerende gebieden van TFBS's. Daarom hebben we de originele dataset (Jolma et al, 2013 ) met extra sequencing om de leesdiepte van de experimenten met bijna een factor 10 te vergroten (van een gemiddelde van

168.000 leest per sequencing-bestand naar

1.656.000 keer gelezen). Experimentele gegevens werden gefilterd door strenge kwaliteitscontrolecriteria (QC) om gevallen met voldoende bibliotheekcomplexiteit en leesaantallen te identificeren om het bouwen van multiparametrische modellen mogelijk te maken. Een totaal van 218 TF's uit 29 families slaagden voor het eerste filter op basis van hoge variabiliteit en grote steekproefomvang van de gegevens, en in totaal 215 TF's uit 27 verschillende families slaagden voor de QC-stap op basis van regressieprestaties (Fig 1).

Figuur 1. Pijpleiding gebruikt om HT-SELEX te genereren m-woordscores en filterdatasets

Voor elke TF hebben we een kernbindend motief geselecteerd om identificatie van de meest waarschijnlijke bindingsplaats binnen m-woorden en filter oligonucleotiden uit die waarschijnlijk ongebonden zijn. De gebruikte motieven zijn ontleend aan een eerdere studie (Jolma et al, 2013 ). Deze motieven bevatten over het algemeen lange flanken naast de kernconsensussequentie, wat zou voorkomen dat we robuust worden m-woordscores vanwege lage leesdekking voor lange reeksen. Om deze moeilijkheid te overwinnen, gebruikten we motieven uit de catalogus samengesteld door Weirauch en Hughes (Weirauch & Hughes, 2011 ) om alleen de kernposities te identificeren en te gebruiken. We hebben de bindingsscore voor elk berekend m-woord dat het kernmotief in het midden bevatte (waardoor enkele mismatches mogelijk waren) en eventuele flankerende sequenties 5' en 3' van het motief. We hebben geprobeerd de mogelijkheid van coöperatieve TF . te vermijdenDNA-binding, waarbij meerdere kopieën van de TF verschillende DNA-bindingsplaatsen (BS's) op dezelfde sequentie bezetten, en om ruis te minimaliseren die wordt veroorzaakt door onnauwkeurige uitlijning van m-woorden gebaseerd op het kernmotief. We hebben dus HT-SELEX-lezingen uitgesloten die meerdere instanties van de kernmotieven bevatten.

Vervolgens hebben we afgeleid m-woordbindingsscores op basis van waargenomen experimentele verrijking. Elk HT-SELEX-experiment omvatte verschillende ronden van selectie van de bindingsplaats (BS) door de TF, waarbij de bindingsspecificiteit van geselecteerde DNA-sequenties in elke ronde toenam. We berekenden de m-woordscore als de verhouding van de frequentie van de m-woord in ronde l over de geschatte frequentie in de eerste ronde, met behulp van een vijfde-orde Markov-model (Slattery et al, 2011 ). De uiteindelijke output van dit proces was de m- woordscores van de kernsequentie en zijn flanken voor elk HT-SELEX-experiment (bijlage Fig. S1A).

Om de nauwkeurigheid van onze m-woordscoreschema en de waarde van diepere sequencing, vergeleken we scores afgeleid door HT-SELEX met die gemeten door genomische-context PBM's (gcPBM's). De gcPBM's gebruiken arrays die speciaal zijn ontworpen met de kernsequentie in het midden, geflankeerd door een genomische context (Gordân et al, 2013 ). Deze sondes zijn bedoeld om het effect van flankerende sequenties te meten en bieden daarom een ​​nauwkeurige gouden standaard voor lange tijd m-woord (m ≥ 12) bindende scores. Het enige eiwit waarvoor zowel gcPBM als HT-SELEX experimentele gegevens bestaan, was het Max homodimeer (Zhou et al, 2015 ). Bijlage Fig S1B toont de goede correlatie (R = 0,64) van 12-woordscores geproduceerd door de twee technologieën, wat de nauwkeurigheid van ons proces bij het produceren aantoont m-woordscores van HT-SELEX-gegevens. Om te testen hoeveel we winnen met betrekking tot gcPBM-bindingsscores door de nieuwe gegevens te gebruiken, hebben we drie verschillende m-woordscores: frequentie, verhouding ten opzichte van de beginronde en verhouding ten opzichte van de geschatte beginronde. Een diepere sequencing verbeterde de correlatie van deze drie scores met gcPBM-scores van 12 woorden, en de verhouding tot geschatte score bereikte de hoogste correlatie (bijlage figuur S1C). Met name bij het verwerken van de eerder gepubliceerde gegevens in (Jolma et al, 2013 ) met dezelfde pijplijn, passeerden slechts 22 eiwitten de kwaliteitscontrole, vergeleken met 218 met de hogere dekking, wat het voordeel van diepere sequencing aantoont.

Hoofdcomponentenanalyse (PCA) onthult TF-familiespecifieke DNA-bindingsspecificiteiten en heterogeniteiten binnen TF-families

We hebben PCA uitgevoerd om TF-familiespecifieke DNA-bindende specificiteiten te visualiseren. De DNA-bindingsvoorkeur van elke TF werd weergegeven door het DNA m-woord met de hoogste bindingsaffiniteit voor deze TF. We hebben dit gecodeerd m-woord in numerieke kenmerkvectoren die (i) alleen mononucleotide (d.w.z. 1-meer) kenmerken, en (ii) zowel 1-meer als DNA-vormkenmerken bevatten. Kenmerken van de DNA-vorm omvatten kleine groefbreedte (MGW), rol, propellertwist (ProT) en helix twist (HelT) en worden voorspeld met onze DNAshape-benadering (Zhou et al, 2013 ). Figuur 2A en B tonen de eerste twee hoofdcomponenten die zijn verkregen met behulp van elke kenmerkvector.

Figuur 2. PCA onthult verschillende DNA-bindingsspecificiteiten tussen TF-families

  1. PCA met behulp van 1-mer-functies. Elke stip vertegenwoordigt een TF. Stippen van dezelfde kleur behoren tot dezelfde TF-familie. Voor elke TF-familie werd een ellips getekend. De ellips is een contour van een gepaste twee-variate normale verdeling die een kans van 0,68 omsluit (R-pakket standaard).
  2. PCA met behulp van 1-mer en vormkenmerken, geannoteerd op dezelfde manier als beschreven in (A).
  3. Boxplots van inter- en intra-familie TF-afstanden afgeleid van (A). Het verschil tussen de medianen van afstanden tussen en binnen de familie is 2,02 (rood).
  4. Boxplots van inter- en intra-familie TF-afstanden afgeleid van (B). Het verschil tussen de medianen van afstanden tussen en binnen de familie is 3,68 (rood).

Verschillende TF-families hadden de neiging om verschillende clusters te vormen in de PCA-spreidingsdiagrammen. Om de clusteringkwaliteit in de twee plots te vergelijken, hebben we de tweedimensionale Euclidische afstanden tussen alle paren TF's uit Fig 2A en B verkregen. Afstanden werden geclassificeerd als intra- of interfamilie en gevisualiseerd als boxplots (Fig 2C en D). Afstanden tussen gezinnen waren over het algemeen groter dan afstanden binnen gezinnen. Toen we zowel 1-mer- als DNA-vormkenmerken gebruikten, was het verschil tussen de medianen van de inter- en intra-familiegroepen iets groter dan het verschil dat werd verkregen bij het gebruik van alleen 1-mer-kenmerken (Fig. 2C en D). Dit resultaat was consistent met figuur 2A en B, wat aangeeft dat meer variantie kon worden verklaard door het introduceren van DNA-vormkenmerken, deels vanwege de betere scheiding van de homeodomeinfamilie (figuur 2B). Om te testen of dergelijke effecten eenvoudigweg te wijten waren aan de hogere dimensionaliteit die werd geïntroduceerd door de extra DNA-vormkenmerken, hebben we willekeurig gegenereerde vormkenmerken toegevoegd op basis van Gauss-verdeling met gemiddelde en standaarddeviatie van de oorspronkelijke vormkenmerken. Zowel de verklaarde variantie als de afstand tussen intra- en intergezinsgroepen waren lager in deze test (bijlage figuur S2).

DNA-vormkenmerken verbeteren de modellering van DNA-bindende specificiteiten in TF-families

We hebben het belang van de herkenning van DNA-vorm door elke TF getest door middel van kwantitatieve modellering van DNA-bindingsspecificiteiten en vergelijking van modelprestaties in termen van de R 2 tussen voorspeld en experimenteel m-woordscores. Vergelijkbaar met de methodologie in Yang et al ( 2014 ) en Zhou et al (2015), hebben we regressiemodellen gebouwd die alleen DNA-mononucleotide-kenmerken gebruikten (d.w.z. 1mer-modellen) of die DNA-mononucleotide- en vormkenmerken combineerden (d.w.z. 1mer+shape-modellen). Een resultaat waarin het 1mer+shape-model beter presteert dan het 1mer-model, geeft aan dat het uitlezen van de DNA-vorm een ​​rol kan spelen bij TF-binding.

Op basis van een analyse van 215 TF's uit 27 verschillende families, ontdekten we dat 1mer+shape-modellen over het algemeen beter presteerden dan 1mer-modellen (Fig 3A), wat wijst op de prevalentie van DNA-vormuitlezing in verschillende TF-families (voor een volledige lijst van datasets die worden gebruikt in Fig. , zie tabel EV1). Omdat het uitlezen van DNA-sequenties een dominante rol speelt bij TF-binding, varieerde het belang van DNA-vormherkenning als extra bijdrage zowel tussen als binnen TF-families. De modelprestaties voor homeodomein-TF's waren bijvoorbeeld over het algemeen aanzienlijk verbeterd dan voor C2H2-TF's. Binnen de homeodomein TF-familie was er een grote variatie tussen individuele leden. Van homeodomein- en bHLH-TF's is eerder waargenomen dat ze gevoelig zijn voor DNA-vormkenmerken (Slattery et al, 2011 Gordon et al, 2013 Yang et al, 2014 Zhou' et al, 2015 ). Hier hebben we deze observatie bevestigd en uitgebreid tot de bZIP-, CENPB-, CP2-, CUT-, ETS-, HSF-, IRF-, MYB-, NFAT-, nucleaire receptor-, PAX-, POU-, PROX-, TBX- en TEA TF-families. Ten minste de helft van de leden in elk van deze families, gedekt door onze gegevens, vertoonde een prestatieverbetering van meer dan 10% wanneer DNA-vormkenmerken aan het model werden toegevoegd. Sommige families waren echter ondervertegenwoordigd in de gegevens met slechts één TF aanwezig (Tabel EV1 voor volledige namen en gedetailleerde informatie over de TF-families, zie Tabel EV2).

Afbeelding 3. Prestatievergelijkingen tussen modellen met verschillende functies

  1. Vergelijking tussen 1mer en 1mer+shape modellen.
  2. Vergelijking tussen vormmodellen die gebaseerd zijn op de originele DNAshape-methode (Zhou et al, 2013 ) en willekeurig geschudde pentameer-querytabellen.
  3. Vergelijking tussen 1mer+2mer+3mer en 1mer+shape modellen.
  4. Vergelijking tussen 1mer+2mer+3mer en 1mer+shape+3merE2 modellen. Het label 3merE2 vertegenwoordigt 3mer-kenmerken van de twee eindposities aan het 5'- en 3'-uiteinde van elke DNA-sequentie.
  5. Vergelijking tussen 1mer+2merNoE2+3merNoE2 en 1mer+shape-modellen. De labels 2merNoE2 en 3merNoE3 geven aan dat respectievelijk 2mer- en 3mer-kenmerken uit de eindposities zijn verwijderd.
  6. Vergelijking tussen 1mer+shape en 1mer+shape+3merE2 modellen.

Om de robuustheid van de experimentele gegevens en onze computationele pijplijn te testen, herhaalden we de bovenstaande analyse op herhaalde experimentele gegevens voor drie TF's uit de bHLH- en homeodomeinfamilies. Onze resultaten toonden consequent bijdragen van de uitlezing van de DNA-vorm voor deze twee families (bijlage Fig S3A). Om te testen of de prestatiewinst eenvoudigweg het resultaat is van het toegenomen aantal modelparameters als gevolg van de toegevoegde DNA-vormkenmerken, hebben we de querytabel voor DNA-vormkenmerken geschud. Vormmodellen die zijn gebaseerd op de geschudde querytabel hebben over het algemeen slechtere prestaties dan die op basis van de oorspronkelijke querytabel (Fig 3B). We hebben ook getest of de resultaten robuust waren voor de motiefzaden die werden gebruikt tijdens de voorverwerking van gegevens. We herhaalden de bovenstaande analyses met behulp van de Weirauch- en Hughes-zaden (Weirauch & Hughes, 2011) als de uiteindelijke zaden in plaats van ze te gebruiken voor het identificeren van de kernposities van de op HT-SELEX gebaseerde motieven gepubliceerd door Jolma et al (2013). We berekenden de correlatiecoëfficiënten van Pearson tussen de prestaties van modellen die waren gebaseerd op de zaden van Weirauch en Hughes (Weirauch & Hughes, 2011) en de Jolma et al (2013) zaden. De hoge correlatie tussen de twee sets motiefzaden gaf aan dat de resultaten robuust waren voor de keuze van motiefzaden (bijlage Fig. S3B). We hebben ook de robuustheid van de resultaten getest onder kleine veranderingen in de mismatch-drempel (zie 4) en de lengte van de flankerende regio's. Beide tests toonden een hoge correlatie tussen verschillende parameterinstellingen, wat voldoende robuustheid aantoont (bijlage Fig. S3C en D).

De homeodomein-TF's in deze studie binden vermoedelijk DNA als monomeren, terwijl onze eerdere studies het belang van de DNA-vorm voor Exd-Hox-heterodimeren (Slattery) aantoonden. et al, 2011 ). Röntgen- en nucleaire magnetische resonantie (NMR)-structuren van DNA-bindende domeinen van het homeodomein in een complex met DNA laten herhaaldelijk zien dat de N-terminale staart van het DNA-bindende domein van het homeodomein interageert met het DNA via kleine groeven en ruggengraatcontacten, wat een signature of DNA shape readout (Joshi et al, 2007 ).

DNA shape features in flanking regions are important for different TF families

We previously observed that 1mer+2mer+3mer models usually outperform 1mer+shape models (Zhou et al, 2015 ). Here, we gained additional clues for possible explanations of this observation. As noted previously (Zhou et al, 2015 ), both 2-mer and 3-mer features are indirect representations of DNA shape characteristics. The 2-mer features describe stacking interactions between adjacent base pairs, whereas 3-mer features describe short structural elements, such as A-tracts that tend to form narrow minor groove regions. Thus, it is not surprising that 1mer+2mer+3mer models can capture TFDNA binding specificities with high accuracy.

Using our high-quality HT-SELEX data, we observed that, for most TFs, 1mer+2mer+3mer models outperformed 1mer+shape models (Fig 3C). As our prediction of local DNA shape features was based on a sliding window of 5 base pairs (Zhou et al, 2013 ), we were unable to predict shape features for the two extreme positions at the 5′ and 3′ ends of each DNA sequence. This limitation could give an edge to 1mer+2mer+3mer models. However, we could encode 2-mer and 3-mer features for those terminal positions, which in turn would work as a proxy for DNA shape. To test this hypothesis, we added 3-mer features from only the two end (E2) positions (i.e., 3merE2 features) to the 1mer+shape model. Performance of the resulting 1mer+shape+3merE2 model was indeed comparable to that of the 1mer+2mer+3mer model (Fig 3D). As an additional test, we removed 2-mer and 3-mer features at the end positions from the 1mer+2mer+3mer model, which resulted in the 1mer+2merNoE2+3merNoE2 model that showed similar performance to the 1mer+shape model (Fig 3E).

We also hypothesized that if longer flanking sequences were available for predicting shape features, then 1mer+shape models would perform similar to 1mer+2mer+3mer models without adding 3merE2 features. To verify this possibility, we used an independent dataset generated by the gcPBM platform (Zhou et al, 2015 ). As expected, 1mer+shape models performed comparable to 1mer+2mer+3mer models for the data without additional 3merE2 features (Appendix Fig S3E). These results imply that DNA shape features in the flanking regions contribute to TFDNA binding specificities, which was previously known for bHLH TFs (Gordân et al, 2013 Yang et al, 2014 Zhou et al, 2015 ). Here, we showed for the first time that this phenomenon is of general nature, as adding 3merE2 features as proxy for missing DNA shape features consistently improved the model performance for various TF families (Fig 3F).

Beyond better interpretability of shape-augmented models, an important distinction between the models is the different number of features required to achieve similar performance. The 1mer+shape model requires 12 features (including second-order DNA shape features) per nucleotide position compared with the 84 features required by the 1mer+2mer+3mer model per nucleotide position (Zhou et al, 2015 ). Although we previously included lower-order 1-mers and 2-mers in our 1mer+2mer+3mer models for reasons of interpretability, nevertheless, the 3-mer features actually contain all of the information of the 1-mers and 2-mers. Thus, a 3mer model is equivalent to a 1mer+2mer+3mer model (4 and Appendix Fig S3F). This choice, however, would still leave the 3mer model with 64 required features per nucleotide position compared with a maximum of only 12 features in the 1mer+shape model.

Feature selection can provide insights into TF–DNA readout mechanisms

We performed feature selection to identify BS positions where DNA shape features contribute to TF-binding specificities. The method is similar to the one we previously introduced for the analysis of SELEX-seq data for Hox proteins (Abe et al, 2015 ). For each TF, we evaluated the R 2 performance of the baseline 1mer model, denoted . Next, we evaluated models that combined 1-mer features with DNA shape features individually at single nucleotide positions l, denoted 1mer+shapel modellen. We denoted the performance as . We calculated the difference in model performance for each nucleotide position l (Fig 4A). De ratio indicates the percentage change in performance due to the availability of DNA shape features at nucleotide position l, with a positive ratio suggesting performance gain. The ratio at position l compared with other positions reflects the relative importance of DNA shape features at different nucleotide positions. We visualized the ratio as a function of position l for each TF in the form of a heat map (Fig 5A and Appendix Fig S4).

Figure 4. Schematic representation of feature-selection process

  1. Feature-selection scheme for adding DNA shape features at one individual position to a sequence-only model.
  2. Feature-selection scheme for removing DNA shape features from one single position from a shape-only model.

Figure 5. Importance of DNA shape features as a function of nucleotide positions revealed by feature selection with machine learning

  1. Heat map based on adding DNA shape features to a sequence-only model.
  2. Heat map based on removing DNA shape features from a shape-only model.
  3. Combined heat map that takes cell-by-cell minimum of heat maps in (A and B).

To avoid interference from DNA sequence information, we devised a second feature-selection approach in which we removed DNA shape features at individual positions from a shape-only model. De ratio was then used for generating the heat map (Figs 4B and 5B, and Appendix Fig S4), where . These two different approaches can sometimes yield conflicting heat maps as discussed below. To address such cases and facilitate the use of these heat maps, we also generated a combined heat map based on the cell-by-cell minimum of the two heat maps (Fig 5C and Appendix Fig S4). Quantitative information about the importance of the position-dependent DNA shape in TFDNA recognition at single-base pair resolution provides the means to determine the structural proteinDNA readout mechanisms based on sequence data. To achieve this goal, we further expanded our feature-selection method to test each individual DNA shape feature category, which enabled us to gauge the importance of each DNA shape feature, that is, MGW, Roll, ProT, or HelT, at every position (Appendix Fig S5). To date, obtaining such information required experimentally solved structures.

Figure 5 shows the position-dependent DNA shape importance for homeodomain TFs that recognize a TAAT motif. For most of these TFs, DNA shape was more important at the 3′ side of the core motif, as indicated by the darkness of colors (Fig 5). Homeodomain TFs that recognize a different motif, for example, TCRTAAA, were shown to have a different positional DNA shape preference (Appendix Fig S4F). Positional preferences were also protein-family specific. For example, for bHLH TFs DNA shape features in both flanking regions were important, whereas for nuclear receptors that bind to an ACANNNTGT motif the central motif region was generally important (Appendix Fig S4A and H). In comparison, bZIP TFs that bind to a TTRCGC motif and homeodomain TFs were generally sensitive to DNA shape features at only one flanking side of the core motif (Appendix Fig S4B and F).

The exact positions where DNA shape features are important were not unambiguously pinpointed for the bHLH TFs and the nuclear receptors that bind to an ACANNNTGT motif (Appendix Fig S4A and H). Both Appendix Fig S4A and H relate to a scenario where the red heat map shows prominent shape effects in multiple consecutive positions, whereas the blue heat map shows almost no effects. We believe that this is due to false positives in the red heat map, that is, positions that are not important for shape readout but identified as such, and false negatives in the blue heat map, that is, positions that are important for shape readout that were not identified. We conclude in this case that DNA shape is important in some positions in the consecutively red regions, but we failed to locate it, even with the help of the blue heat map.

We illustrated the relevance of feature importance heat maps derived from feature-selection approaches by considering experimental structures of the homeodomain proteins PITX2 (PDB ID 2LKX) and GBX1 (PDB ID 2ME6) in complex with DNA (Fig 6A and B). These structures provide possible explanations for entries representing PITX3 and GBX1 on the heat maps (Fig 5). As no experimental structure for PITX3 is available, we used an NMR structure for PITX2 (Chaney et al, 2005 ), which shares the same DNA-binding domain as PITX3. In the heat maps, PITX3 has darker colors at the 3′ side of the TAAT motif, indicating a more important role of DNA shape at these positions. In the PITX2 structure, the N-terminal tail of the protein interacts with DNA in the minor groove of the TAAT motif. The structure contains a narrow minor groove region near the second A within the TAAT motif (Fig 6A). In this case, the protein might exploit the DNA structural characteristics at positions highlighted in the heat maps to achieve its binding specificity.

Figure 6. Three-dimensional structure and DNA sequence and shape logos for the homeodomain TFs PITX2/PITX3 and GBX1

  1. NMR structure of PITX2 in complex with DNA (PDB ID 2LKX) and the CURVES (Lavery & Sklenar, 1989 ) derived plot for the MGW of the bound DNA.
  2. NMR structure of GBX1 in complex with DNA (PDB ID 2ME6) and the CURVES (Lavery & Sklenar, 1989 ) derived plot for the MGW of the bound DNA.
  3. DNA sequence and shape logos for PITX3.
  4. DNA sequence and shape logos for GBX1.

We observed similar concurrence between heat map and structural analyses for the TF GBX1, where the structure has a narrow minor groove region at the 3′ flank (Fig 6B). Although the positions indicated by the heat maps do not match the positions in the structure in an exact way, the heat maps successfully highlighted those nearby positions. Moreover, the heat maps were consistent with our conclusion that DNA shape features in flanking regions are important for TFDNA binding specificities (Fig 3D–F). In addition to the homeodomain family, we used a structure of the human progesterone receptor (PDB ID 2C7A) from the nuclear receptor family to illustrate how the heat maps can provide hints to the structural mechanisms of proteinDNA binding. In the structure (Roemer et al, 2006 ), MGW, Roll, and ProT show distinct characteristics in the central region of the DNA-binding site, which potentially explains the central “red” regions in the heat maps (Appendix Fig S6).

DNA shape logos represent structural readout mechanisms

To visualize the detailed DNA shape preferences of individual TFs, we propose a new visualization, DNA shape logos, analogous to sequence logos for PWMs. In these logos, we used the letters H, M, P, and R to represent DNA shape features HelT, MGW, ProT, and Roll, respectively. The height of each letter indicates the importance derived from the feature-selection analysis for the corresponding DNA shape feature at a specific position (Fig 6). As an example, we used ΔR 2 , that is, the performance gain due to adding an individual DNA shape feature to a 1mer model, to generate shape logos for PITX3 and GBX1 (Fig 6C and D). For PITX3, a prominent M at positions 7, 8, 9, and 10 overlaps with the narrow minor groove region in the structure. Similarly, for GBX1, a prominent M at positions 7 and 8 overlaps with the narrow minor groove in the structure. DNA shape information was missing for the two nucleotide positions at each end of the TFBS thus, no letters are shown at these positions in the shape logo. DNA shape logos can facilitate the integration of structural information in motif finding tools. Sequence and shape logos for all the TFs studied in this work are provided as Datasets EV1 and EV2, respectively.


4 Contact:

Repetitive elements in DNA sequences consist two or more copies of approximate patterns of nucleotides and are abundant in both prokaryotic and eukaryotic genomes. Over two-thirds of the human genome and 5 - 10 % bacterial genomes are repetitive regions (de Koning et al., 2011) . Repetitive elements play important roles in genome structure and functions such as nucleoprotein complex formation, chromosome structure, and gene expression. Various diseases including cancer and neurodegentive disease can also arise from changes of repetitive elements. The distribution of repetitive DNA sequences can be used as fingerprints of bacterial genomes (Versalovic et al., 1991) and human individuals.

Repetitive elements are complex structures. They may exist as imperfect tandem repeats, insertion and deletions in repeats, interspersed repeats, and palindromic sequences, etc. These partial and hidden repeat signals in DNA sequences are difficult to analyze through straightforward observation and sequence comparison.

Currently, repetitive elements and hidden periodicities of DNA and protein sequences are primarily detected by digital signal processing and statistical approaches (Treangen and Salzberg, 2011) . In most signal processing methods, DNA sequences are converted to numerical sequences, and the hidden periodicities arising from repetitive elements can be identified by Fourier power spectrum at specific periodicities (Yin and Wang, 2016)

. Commonly used signal processing methods by Fourier transform include SRF maps

(Sharma et al., 2004) , spectral analysis (Buchner and Janjarasjitt, 2003) , Ramanujan-Fourier transform (Yin et al., 2015) , and the periodic power spectrum method (Yin and Wang, 2016) . The statistical methods are based on distribution analysis of nucleotides in DNA sequences. The common statistical methods for repeat findings are tandem repeats finder (Benson, 1999) and statistical spectrum (Epps et al., 2011)

(Arora and Sethares, 2007) , and information decomposition (Korotkov et al., 2003) . Besides signal processing and statistical approaches, sequence alignments such as RepeatMask are also used to identify repetitive patterns in genomes, and but require a known reference repeat sequence.

Despite significant advances in repeat finding, it is still difficult to precisely capture the essential features of repetitive elements such as consensus patterns, perfect levels and copy numbers of repeats. For example, while Fourier transform is the most common used approach for finding repeats, it may not exactly correlate the strength of Fourier power spectrum with the perfect level of repeat patterns. Furthermore, since Fourier power spectrum is weak for short DNA sequences and long harmonious periodicities are embedded in short periodicities, Fourier transform can not capture repeats in short DNA sequences and long harmonious periodicities. Moreover, the relationship between repetitive elements and periodicities of genomes is not fully understood. Thus there is a high potential for improving the accuracy for identifying repetitive elements and better understanding the relationship of periodicities and repeats in DNA sequences (Suvorova et al., 2014 Epps et al., 2011 Illingworth et al., 2008) .

In this paper, we present an ab initio method to quantitatively identify repetitive sequences and periodicities in DNA sequences. The method is based on the nucleotide distribution uniformity at periodic positions in DNA sequences or genomes. The distribution uniformity of nucleotides reflects the unbalance of nucleotide frequencies on periodic positions and thus can indicate the strength for periodic signals in DNA sequences. The method can also reveal the consensus repeat pattern for the major periodicity of DNA sequences, and quantitatively determine the perfect level and copy numbers of repetitive sequences. The proposed method also formulates the relationship between repetitive elements and the corresponding periodicities in DNA sequences.


Experimental procedures

Bacterial strains and growth conditions

Helicobacter pylori strains (Table S1) were grown on solid horse blood agar (HB) plates containing 4% Columbia agar base (Oxoid), 5% defibrinated horse blood (HemoStat Laboratories), 0.2% β-cyclodextrin (Sigma), 10 µg ml −1 vancomycin (Sigma), 5 µg ml −1 cefsulodin (Sigma), 2.5 U ml −1 polymyxin B (Sigma), 5 µg ml −1 trimethoprim (Sigma), and 8 µg ml −1 amphotericin B (Sigma) at 37°C either under a microaerobic atmosphere generated using a CampyGen sachet (Oxoid) in a gas pack jar or in an incubator equilibrated with 14% CO2 and 86% air. For liquid culture, H. pylori was grown in Brucella broth (Difco) containing 10% fetal bovine serum (BB10, Invitrogen) with shaking in a gas pack jar containing a CampyGen sachet. For resistance marker selections, bacterial media were additionally supplemented with 15 µg ml −1 chloramphenicol (Cm, Sigma), 25 µg ml −1 kanamycin (Kan, Fisher Scientific) 2.5 µg ml −1 erythromycin (Ery, Fisher Scientific) or 36 µg ml −1 metronidazole (Mtz, Sigma).

DNA manipulations

DNA manipulations, such as restriction digestion, PCR and agarose gel electrophoresis, were performed according to standard procedures ( Ausubel et al., 1997 ). H. pylori genomic DNA (gDNA) was prepared by Wizard genomic DNA preparation kits (Promega). Primers used for PCR and sequencing are described in Table S2. Plasmid DNA (Table S3) was isolated and prepared from E coli using Qiagen Maxiprep kit (Qiagen). The FHCRC Genomics Shared Resource performed the sequencing of plasmid DNA and PCR products and the resulting sequences were analysed using Sequencher (Gene Codes Corporation).

Generation of H. pylori knockout isogenic mutants

Knockout alleles were constructed in H. pylori NSH57 using a vector-free allelic replacement strategy to generate alleles in which a non-polar kanamycin resistance (aphA3) cassette ( Menard et al., 1993 ), an erm cassette conferring resistance to erythromycin ( Lampson and Parisi, 1986 Dailidiene et al., 2006 ), or a chloramphenicol acetyl transferase (kat) resistance cassette fused to a sucrose sensitivity marker (sacB) ( Copass et al., 1997 Humbert and Salama, 2008 ) replaced 80–90% of the coding sequence of the gene while preserving the start and stop codons. The primers used for this procedure are designated as 1 through 4 and are given in Table S2. After natural transformation with the appropriate PCR product and selection on Kan-, Ery- or Cm-containing media, four clones were evaluated by PCR to confirm replacement of the WT allele with the null allele. De ΔrecJ::kanΔaddA 852-2540 double mutant was generated by transforming strain ΔrecJ::kanΔaddA::catsacB with a PCR product digested with SspI (New England Biolabs) and ligated with T4 DNA ligase (Invitrogen) to delete a 1.7 kbp intergenic region in addA. Transformants were selected on sucrose-containing HB plates, screened on Cm-containing media and checked by PCR to confirm the addA deletion. Urease activity and flagella-based motility were confirmed for all the clones generated. Single clones were used for transformation experiments.

Generation of H. pylori complemented mutants

Constructs for chromosomal complementation at the rdxA locus were made by cloning each gene individually into pLC292 ( Terry et al., 2005 ), which were then introduced into H. pylori NSH57 by natural transformation and selection on Mtz-containing media ( Dailidiene et al., 2006 ). Each gene was amplified using primers -XbaI and -SalI (Table S2) from H. pylori NSH57 gDNA using high-fidelity Taq polymerase (Platinum Taq, Invitrogen). The resulting PCR product was digested with XbaI and SalI (New England Biolabs), ligated into pLC292, and electroporated in E coli strain DH10B or XA90 ( Ezaz-Nikpay et al., 1994 ) for pOH10 (Table S3). All inserted genes were fully sequenced and contained the expected nucleotide sequences.

Natural transformation

To generate knockout and complemented mutant strains of H. pylori, bacteria were freshly grown for 24–32 h on HB plates, transferred as patches onto fresh plates and grown for an additional 6–8 h. DNA (plasmid or PCR product) was diluted as appropriate in distilled water and 10 µl was added to each patch and incubated overnight. The mixture was harvested from the plate surface, resuspended in 350 µl phosphate-buffered saline (PBS) and plated onto selective HB plates.

To assess the frequency of natural transformation, recipient H. pylori bacteria freshly grown on HB plates were resuspended in 350 µl BB10 media and used to inoculate a 5 ml liquid culture grown for 6–8 h. The optical density at 600 nm (OD600) of this culture was measured and the culture was diluted back to OD600 0.015 to reach logarithmic phase of growth (OD600∼1) after overnight incubation. One hundred microlitres of recipient bacteria was dispensed in a flat-bottom 96-well plate and transformed in duplicates or triplicates with 10 µl of 1 ng µl −1 donor gDNA. Donor gDNA was constructed by inserting the kat resistance cassette at bp 483 in gene cagH van H. pylori strain NSH57 and J99 (hpG27-499 and jhp0489 respectively). To measure transformation of the ΔdprA mutant, donor gDNA was isolated from the G27 cag2::aphA3-sacB clone ( Pinto-Santini and Salama, 2009 ). After 3 h incubation, 50 µl and 5 µl of the mixture were plated on Cm or Kan HB plates and 20 µl of a 10 −5 dilution was plated on plain HB plates to determine the total number of viable bacteria. Transformation frequency was calculated as the number of Cm or Kan resistant colonies per colony-forming unit.

In the co-culture experiment, NSH57 and J99 ΔcomB10::ermΔcagH::cat were used as donor strains and to maximize DNA released in the culture media, we grew donor bacteria to stationary phase before mixing them with the recipient strain. ΔcomB10 strains show no detectable transformation ( Dorer et al., 2010 ) ensuring unidirectional transformation in the co-culture assay. Recipient strains NSH57 hp0203-hp0204::aphA3 ( Langford et al., 2006 ) and Δhpy188IIIR::aphA3ΔhpyCH4VR::erm were grown to logarithmic phase as described above and mixed at equal volume with the donor strains in a flat-bottom 96-well plate. After 3 h co-incubation, 100 µl of the mixture was plated on Cm + Kan HB plates to select for recombinant clones and 20 µl of a 10 −5 dilution was plated on Kan HB plates to determine the total number of recipient bacteria.

Mapping of integration end-points

Chromosomal DNA of the transformants was prepared and 5–7 kbp of the regions upstream and downstream of the kat marker were amplified by PCR using primer pairs -6FcagH/cagH::cat-3 and cagH::cat-4/5RcagH (Table S2) respectively. The resulting PCR products were purified with the DNA clean and concentrator-5 kit (Zymo Research) and digested with the appropriate restriction enzymes for a minimum of 4 h (New England Biolabs) or sequenced by the FHCRC Genomics Shared Resource.

Sensitivity to UV and antimicrobial agents

UV sensitivity assays were carried out as described previously ( Amundsen et al., 2008 ). For antimicrobial sensitivity testing, H. pylori were grown overnight in liquid culture to OD600 = 0.3, and 200 µl was plated on solid medium lacking all other antimicrobials, and incubated for 30 min in a CO2 incubator. E-test strips (AB Biodisk) were then placed on the plates, which were further incubated for two days and read according to the manufacturer's instructions.

Statistical analysis

EEN t-test was used to compare the mean of integration lengths or transformation frequency between WT bacteria and mutant clones and those comparisons resulting in a P-value of < 0.05 were considered significant. All statistical analyses were performed using the SAS version 9.1 software (SAS Institute, Cary, NC, USA).

In silico genomic analysis

Helicobacter pylori sequences were retrieved from the H. pylori genome browser http://hpylori.ucsc.edu/. Voor H. pylori strain NSH57, the sequence of the parent strain G27 was used ( Baltrus et al., 2009 ). The distribution of restriction sites and single nucleotide polymorphism was analysed with Sequencher (Gene Codes Corporation).


6.4: Restriction Mapping

  • Contributed by Michael Blaber
  • Professor (Biomedical Sciences) at Florida State University

De restriction/modification system in bacteria is a small-scale immune system for protection from infection by foreign DNA.

In the late 1960's it was discovered that E coli contains enzymes that will methylate specific nucleotide bases in DNA

· Different strains of E coli contained different types of these methylases

  • Typical sites of methylation include the N6 position of adenine, de N4 position of cytosine, of de C5 position of cytosine.

Figure 6.4.1:Methylation sites

  • In addition, only a fractional percentage of bases were methylated (i.e. not every adenine was methylated, for example) and these occurred at very specific sites in the DNA.
  • A characteristic feature of the sites of methylation, was that they involved palindroom DNA sequences.
  • Here is an example from a particular E. coli strain R1:

Figure 6.4.2:Palindromic DNA

(EcoR1 methylase specificity. Rubin and Modrich, 1977)

  • In addition to possessing a particular methylase, individual bacterial strains also contained accompanying specific endonuclease activities.
  • The endonucleases cleaved at or near the methylation recognition site.

Figure 6.4.3:Cleavage near methylation site

  • These specific nucleases, however, would niet cleave at these specific palindromic sequences if the DNA was methylated.

Thus, this combination of a specific methylase and associated endonuclease functioned as a type of immune system for individual bacterial strains, protecting them from infection by foreign DNA (e.g. viruses).

  • In the bacterial strain EcoR1, the sequence GAATTC will be methylated at the internal adenine base (by the EcoR1 methylase).
  • The EcoR1 endonuclease within the same bacteria will niet cleave the methylated DNA.
  • Foreign viral DNA, which is not methylated at the sequence "GAATTC" will therefore be recognized as "foreign" DNA and zullenbe cleaved by the EcoR1 endonuclease.
  • Cleavage of the viral DNA renders it non-functional.

Such endonucleases are referred to as "restriction endonucleases" because they restrict the DNA within the cell to being "self".

The combination of restriction endonuclease and methylase is termed the "restriction-modification" system.

Since different bacterial strains and species have potentially different R/M systems, their characterization has made available honderden of endonucleases with different sequence specific cleavage sites.

  • They are one of the primary tools in modern molecular biology for the manipulation and identification of DNA sequences.
  • Restriction endonucleases are commonly named after the bacterium from which it was isolated.

Arthrobacter luteus

"Four cutter". Leaves blunt ends to the DNA.

Bacteroides fragilis

"Four cutter". Leaves 5' overhang.

Neisseria cinerea

"Five cutter". Middle base can be either cytosine or guanine. Leaves 5' overhang. Different recognition sites may have non-complementary sequences.

"Six cutter". Leaves 5' overhang. Behaves like a "four cutter" ('star' activity) in high salt buffer. $44 for 10,000 units.

Haemophilusaegyptius

"Six cutter". Pu is any purine, Py is any pyrimidine. Leaves 3' overhang.

"Seven cutter". Pu is any purine, Py is any pyrimidine, N is any base. Leaves 5' overhang. Different recognition sites may have non-complementary sequences.

"Six cutter with interrupted palindrome". Leaves 5' overhang. Different recognition sites may have non-complementary sequences.

Bacillusstearothermophilus

"Six cutter". Different recognition sites zullen be complementary.

Acetobacter aceti

"Six cutter" with 3' overhang. Same recognition sequence as Bsa HI, but different cleavage position.

Non-palindrome, distal cleavage. Leaves 3' overhang. $50 for 50 units.

Nocardiaotitidiscaviarum

"Eight cutter". Leaves 5' overhang.

Bacillusstearothermophilus

  • The utility of restriction endonucleases lies in their specificity and the frequency with which their recognition sites occur within any given DNA sample.
  • If there is a 25% probability for a specific base at any given site, then the frequency with which different restriction endonuclease sites will occur can be easily calculated (0.25 n ):

Frequency of Occurrence

1 Alu site in every 256 bases (0.25 Kb)

1 Nci I site in every 1024 bases (1.0 Kb)

1 EcoR1 site in every 4,096 bases (4.1 Kb)

1 EcoO109I site in every 16,384 bases (16.4 Kb)

1 Not I site in every 65,536 bases (65.5 Kb)

Thus, on average, any given DNA will contain an Alu I site every 0.25 kilobases, whereas a Not I site occurs once about every 65.5 kilobases.

  • Not I is therefore a very useful enzyme for isolating large regions of DNA, typically in research involving genomic DNA manipulations.
  • Alu I would be expected to digest a DNA sample into lots of little pieces.

The assortment of DNA fragments would represent a specific "fingerprint" of the particular DNA being digested. Different DNA would not yield the same collection of fragment sizes. Thus, DNA from different sources can be either matched or distinguished based on the assembly of fragments after restriction endonuclease treatment. These are termed "Restriction Fragment Length Polymorphisms", or RFLP's. This simple analysis is used in various aspects of molecular biology as well as a law enforcement and genealogy. For example, genetic variations that distinguish individuals also may result in fewer or additional restriction endonuclease recognition sites.


Invoering

Comparative sequence analysis has had a major impact on molecular biology and genetics. Comparison of the sequences of protein-coding genes between multiple species has enabled prediction of gene function [1], identification of protein domains [2], prediction of functional amino acid residues [3,4], and detection of signals of natural selection at the level of whole genes [5] and individual codons [6,7]. Inferring non-neutral sequence elements in the human genome is of considerable interest even without a specific a priori hypothesis concerning their possible functional role(s). On a general level, for example, sequence conservation may considerably inform human genetic studies seeking to identify allelic variants associated with disease phenotypes, particularly in noncoding regions [8]. The effect of human SNPs at the level of molecular function and phenotype depends on the importance of the individual nucleotide position, whereas the information of the sequence region as a whole is not necessarily relevant. For example, about half of human SNPs within protein coding genes are represented by synonymous variants, which are likely to be of limited importance, even though they are embedded within highly conserved exonic sequences. In addition, a subset of individual nucleotides conserved in four mammalian genomes were shown to be under selective pressure [9]. A position-specific measure of selective constraint is therefore highly suitable for analysis of positions that are polymorphic within the human population.

Several algorithms have been developed for detection and scoring of sequence conservation in the context of a multispecies sequence alignment. However, to date these approaches have been applied almost exclusively to detect discrete regions with elevated average sequence conservation that typically extend for up to hundreds of contiguous bases [10–14]. Such regions encompass canonical coding exons, as well as so-called “conserved noncoding sequences” that presumably result from purifying selection, and are thereby indicative of functional importance [15,16].

Recently, comparative genomic sequence of unprecedented depth has been generated by sequencing of multiple mammalian and other vertebrate genomes orthologous to 1% of the human genome defined by the ENCODE regions [17,18]. Several alignment techniques have been applied to construct multiple sequence alignments within ENCODE regions [18]. These alignments have in turn been subjected to analysis with existing sequence conservation detection algorithms, including phastCons[10], GERP [11], and BinCons [13]. The conserved regions identified by these analyses show statistically significant overlap with experimentally identified coding and noncoding functional elements. However, the majority of experimentally characterized noncoding functional elements fall outside of currently delineated conserved regions, and, conversely, most conserved regions were located outside of experimentally detected elements [18]. The fact that many functional elements reside in noncoding regions that do not exhibit uniformly high conservation is perhaps not surprising given that binding sites for transcriptional factors that mediate many biological processes are quite plastic evolutionarily [19]. Conversely, many individual nucleotides located outside of well-defined conserved regions exhibit sequence conservation across multiple species. Such conservation may be due to mere chance or, for a certain fraction of these nucleotides, may reflect their importance for fitness and hence function. The aforementioned observations emphasize the need for higher resolution methods for analysis of evolutionary conservation within functional elements and generally across the genome.

Here we develop an approach for analyzing sequence conservation at the individual base-pair level, with an aim toward correlating conservation with human genetic variation and with functional genomic annotations. We present a new probabilistic conservation score, SCONE (Sequence Conservation Evaluation). SCONE provides conservation scores for individual nucleotide positions, and can be applied to predict continuous sequence regions with an elevated level of conservation.

We apply SCONE to the study of annotated functional elements and human sequence polymorphism. We focus on the statistical distribution of position-specific conservation scores rather than on the bulk overlap between conserved regions and functional features. It is clear from the outset that the power to detect conservation at the single base-pair resolution is limited, even when comparing multiple species [20]. We surmount this obstacle by deriving considerable statistical power from combined analysis of numerous individual nucleotide positions from many genomic regions. While this analysis does not allow us to detect individual functional positions accurately, we can show that, collectively, a subset of noncontiguous individual positions are important. A key advantage of the analysis of the distribution of position-specific scores is that it is unbiased with respect to the pattern of conservation along a given sequence region. SCONE thus has the potential to analyze putative functional elements in which the conservation signal is not homogeneous or manifested by exon-like contiguous conserved stretches.

We report herein on the relationship between sequence conservation, functional sequence elements, and human allelic variation, as revealed by single-nucleotide conservation analysis.


Opmerking van de uitgever Springer Nature blijft neutraal met betrekking tot jurisdictieclaims in gepubliceerde kaarten en institutionele voorkeuren.

Extended Data Fig. 1 Highly efficient base editing by A3A-BE4max or hyA3A-BE4max in mouse embryos.

(een, b) Genotyping of F0 generation pups by A3A-BE4max and hyA3A-BE4max. The frequencies of WT and mutant alleles were determined by analyzing HTS using BE-analyzer. The percentage on the right represents the frequency of the indicated mutant allele with the corresponding mutation-induced amino acid conversion shown in parentheses. The frequency of the wild-type allele was omitted. Wt, wild-type.

Extended Data Fig. 2 Off-target analysis and germline transmission of the founders derived from hyA3A-BE4max injection.

(een) HTS was performed with mouse tails to determine editing efficiencies at 15 potential off-target sites in three Dmd mutant F0 mice (#BD03, #BD04 and #BD07). Mismatched nucleotide letters are indicated in lowercase. Data are means ± SD (n = 3 mice).(B) HTS alignments of mutant sequences from F1 generated by mating founder #BD12(♀) with Wt (♂). The column on the right indicates frequencies of mutant alleles. Wt, wild-type.Statistical source data are provided in Source Data Extended Data Fig. 2.

Extended Data Fig. 3 Comparison of base editing efficiency and protein levels by CBEs and hyCBEs in HEK293T cells.

(een)Comparison of base editing efficiency induced by A3A-BE4max or hyeA3A-BE4max in HEK293T cells. The average mutation percentage derived from three independent experiments of A3A-BE4max and hyeA3A-BE4max at the same site is listed. Some of the data (hyeA3A-BE4max) are the same as presented in Fig. 4a. Statistical source data are provided in Source Extended Data Fig. 3. (B) The protein levels of BE4max, hyBE4max, A3A-BE4max, hyA3A-BE4max, eA3A-BE4max and hyeA3A-BE4max were determined by Western blotting in HEK293T cells 3 days after transfection of similar amounts of plasmid DNA. Specific antibodies against Cas9 (top) or GAPDH (bottom) were used. Western blotting images are representative of three independent experiments. Unprocessed blots are shown in Source Data Extended Data Fig. 3.

Extended Data Fig. 4 Comparison of base editing product purity induced by variant base editors in HEK293T cells.

(een) Comparison of base editing products induced by BE4max vs hyBE4max. HTS data were analyzed and the ratio of each type of nucleotides was listed on each target position. Data are means ± SD (n = 3 independent experiments). (B) Comparison of base editing products induced by A3A-BE4max vs hyA3A-BE4max. HTS data were analyzed and the ratio of each type of nucleotides was listed on each target position. Data are means ± SD (n = 3 independent experiments) (C) Comparison of base editing product induced by eA3A-BE4max vs hyeA3A-BE4max. HTS data were analyzed and the ratio of each type of nucleotides was listed on each target position. The individual data points are shown as black (C > T), light green (C > A) and light red (C > G) dots. Data are means ± SD (n = 3 independent experiments). Statistical source data are provided in Source Data Extended Data Fig. 4.

Extended Data Fig. 5 Whole genome sequencing of Dmd F0 (#DD11) and wild-type (Wt) mice.

(een) Summary of genome sequencing analysis. WGS for a Dmd mutant mouse (#DD11) and a wild type mouse (Wt) were performed. A total of 82,573 and 62,359 SNPs were identified for #DD11 and Wt, respectively. After filtering out dbSNP (naturally occurring variants in the SNP database), 20,387 SNPs were obtained in the #DD11 genome. Then the sequences at the remaining SNP sites were compared with all on-/off-target sequences (20 bp). (B) Summary of on-/off-target site information. A total of 175,058 sites, including 1 on-target site and 20 374 2,869 22,335 and 148,569 off-target sites with 3, 4, 5, 6, or 7 mismatch/es, respectively, were analyzed. (C) Summary of the whole-genome sequencing. (NS) Summary of off-target analysis. After comparing the sequences at the remaining SNP sites with the 175,058 on-/off-target sequences (20 bp), the C-to-T substitution was only detected within the on-target sequencing in #DD11. (e) Validation the off-target candidate site determined in (d) using targeted deep sequencing of genomic DNA isolated from various #DD11 organs (heart, liver, lung and tail). Mismatched nucleotides and PAM sequences are shown in red and in blue, respectively. Data represent mean from two independent experiments. Statistical source data are provided in Source Data Extended Data Fig. 5.


Bekijk de video: DNA replication and RNA transcription and translation. Khan Academy (November 2021).