Spørsmål:
Hvordan kan jeg fjerne (ikke-trivielle) duplikater fra en VCF-fil?
terdon
2019-02-27 21:56:46 UTC
view on stackexchange narkive permalink

Dette er relatert til spørsmålet jeg stilte her. Tenk på en vcf-fil som inneholder duplikatvarianter, men der duplikatene ikke bare er det samme i samme notasjon, men i stedet er den ene en delmengde av den andre. For eksempel: 

  ## fileformat = VCFv4.1 ## reference = foo ## FORMAT = <ID = GT, Number = 1, Type = String, Description = "Genotype" > ## contig = <ID = chr12> # CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample1chr12 529514. AACAC AATAC. SENDE    . GT 0 / 1chr12 529516. C T. SENDE    . GT 0/1

Disse de to variantene er faktisk de samme. De resulterer i nøyaktig samme genotype. Å endre AACAC til AATAC i posisjon 529514 betyr bare å endre C til T i posisjon 529516.

Er det noe verktøy som kan oppdage slike duplikater og fjerne dem? Jeg prøvde vcfuniq fra vcflib , men det ser ikke ut til å gjenkjenne dette som et duplikat. Jeg tror det bare ser på de første 4 feltene og anser bare duplikater av disse variantene med nøyaktig de samme verdiene i de første 4 feltene: 

  $ ./bin/vcfuniq test.vcf ## fileformat = VCFv4.1 ## referanse = foo ## FORMAT = <ID = GT, Number = 1, Type = String, Description = "Genotype" > ## contig = <ID = chr12> # CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample1 Sample1. AACAC AATAC. SENDE    . GT 0 / 1chr12 529516. C T. SENDE    . GT 0/1

Men som forklart i spørsmål knyttet til, anser EBIs vcf_validator dette som ugyldig. Og det er egentlig ikke fornuftig å ha disse duplikatene i alle fall, så er det noen måte jeg kan oppdage og fjerne dem? Gjerne et eksisterende verktøy, men jeg er også åpen for skriptløsninger.

Dette kompliseres ytterligere av tilfeller som denne: 

  ## fileformat = VCFv4 .1 ## referanse = foo ## FORMAT = <ID = GT, Number = 1, Type = String, Description = "Genotype" >
## contig = <ID = chr12> # CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample1chr12 529514 529514 AACAC AAT, AATAC 0.00. . GT 0 / 1chr12 529516 529516 C T. SENDE    . GT 0/1

Dessverre vil denne ikke bli fanget av tilnærmingen i Daniels smarte skript: 

  $ cat test2.vcf | foo.py ## fileformat = VCFv4.1 ## reference = foo ## FORMAT = <ID = GT, Number = 1, Type = String, Description = "Genotype" > ## contig = <ID = chr12> # CHROM POS ID REF ALT KVALTFILTERINFORMASJONSFORM Eksempel1chr12 529514 529514 AACAC AAT, AATAC 0,00. . GT 0 / 1chr12 529516 529516 C T. SENDE    . GT 0/1
To svar:
terdon
2019-02-27 22:10:34 UTC
view on stackexchange narkive permalink

Det viser seg at bcftools kan gjøre dette (testet på bcftools-1.8), hvis du gir det referansegenomet å teste mot: 

  $ bcftools norm -d ingen -f hg19.fa test.vcf ## fileformat = VCFv4.1 ## FILTER = <ID = PASS, Beskrivelse = "Alle filtre bestått" > ## referanse = foo ## FORMAT = <ID = GT, Number = 1 , Type = Streng, beskrivelse = "Genotype" > ## contig = <ID = chr12> ## bcftools_normVersion = 1.8 + htslib-1.8 ## bcftools_normCommand = norm -d ingen -f hg19.fa test.vcf; Dato = ons 27. februar 16:08:44 2019 # CHROM POS ID REF ALT KVALTFILTERINFORMASJONSFORMAT Eksempel1chr12 529516. C T. SENDE    . GT 0 / 1Linjer totalt / delt / omjustert / hoppet over: 2/0/1/0

For det mer komplekse tilfellet av multiallelisk variant i det andre VCF-eksemplet fra spørsmålet, du kan kjøre den gjennom bcftools to ganger. Når du bruker norm for å venstrejustere og dele multi-alleliske varianter, og deretter igjen for å fjerne duplikatene: 

  $ bcftools norm -m -any -NO z - O v -o - ~ / test2.vcf | bcftools norm -d ingen -f hg19.faLines totalt / delt / omlagt / hoppet over: 2/1/0/0 ## fileformat = VCFv4.1 ## FILTER = <ID = PASS, Beskrivelse = "Alle filtre passerte" > ## referanse = foo ## FORMAT = <ID = GT, nummer = 1, Type = streng, beskrivelse = "Genotype" > ## contig = <ID = chr12> ## bcftools_normVersion = 1.8 + htslib-1.8 ## bcftools_norm -NO z -O v -o - test2.vcf; Date = Wed Feb 27 18:18:32 2019 ## bcftools_normCommand = norm -d none -f hg19.fa -; Dato = Ons 27. februar 18:18:32 2019 # CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample1chr12 529516 529514 CAC T 0. . GT 0 / 1chr12 529516 529514 C T 0. . GT 0/0 Linjer totalt / delt / omlagt / hoppet over: 3/0/2/0
Stor informasjon å ha i baklommen. Fin løsning!
Skulle ønske jeg kunne +1 dette igjen.
Daniel Standage
2019-02-27 22:16:00 UTC
view on stackexchange narkive permalink

Jeg er ingen ekspert med VCF (få kan si at de er det!), men jeg har jobbet mye med VCF-data de siste årene, begge verktøy for å konsumere og produsere VCF. Jeg har aldri sett varianter kodet på denne måten, og det ser ut til å være ikke-kanonisk. Vanligvis:

  • Single nucleotide variants (SNVs) er kodet med en enkelt base som REF-allelet og en enkelt base som ALT-allelet.
  • I tilfelle innsettinger eller slettinger, den kortere av REF- og ALT-allelene vil være en enkelt base, basen går foran den innsatte / slettede sekvensen. Dermed er den første basen av REF- og ALT-allelene alltid den samme.
  • I det sjeldnere tilfellet av to eller flere påfølgende substitusjoner som danner en multinukleotidvariant (MNV), vil REF- og ALT-allelene ha samme lengde.

Å bruke multi-bp-strenger av samme lengde for å kode SNV er unødvendig og, som du har påpekt, problematisk. Dette får meg til å tro at det er en feil eller en "funksjon" av variantprediktoren som produserte VCF.

I dette tilfellet vil jeg skrive et lite skript som vil se etter varianter der REF og ALT alleler har samme lengde. Hvis basen er den samme for REF og ALT i en hvilken som helst posisjon, kan du slippe den og justere posisjonen deretter.

Skriptet nedenfor konverterer disse funky SNV-ene til den kanoniske representasjonen, og vil også fungere på MNV-er. Standardverktøy bør da arbeide for å fjerne duplikatene. 

  #! / Usr / bin / env python3def kanonikalisere (instream): for linje i instream: hvis ikke line.startswith ('#'): verdier = line.split ('\ t') pos = int (verdier [1]) ref, alt = verdier [3: 5] hvis len (ref) > 1 og len ( ref) == len (alt): # Hvor mange bp å trimme av enden for n, (r, a) i oppsummering (zip (ref [:: - 1], alt [:: - 1])): hvis r! = a: revoffset = -1 * n break # Hvor mange bp for å trimme av fronten
for n, (r, a) i enumerate (zip (ref, alt)): hvis r! = a: offset = n-verdier [1] = str (pos + offset) verdier [3] = ref [offset: revoffset] verdier [4] = alt [offset: revoffset] break line = '\ t'.join (values) yield lineif __name__ ==' __main__ ': import sys for line in canonicalize (sys.stdin): print (line, end = '')

OPPDATERING : Ved ytterligere refleksjon er det mer kompliserte eksemplet du listet opp, fornuftig. I referansen har vi sekvensen AACAC , og de alternative allelene representerer to variasjoner på dette: sletting av de siste to bp (i første omgang), og en punktmutasjon av den midterste C til T (i begge tilfeller). Vanligvis er det bare en enkelt bp som går foran definisjonen på en indel, så jeg ville ha kodet denne komplekse varianten som ref = ACAC alt = AT, ATAC .

Så SNV er "underforstått av" / "kodet i" / "overflødig med" den komplekse varianten, men det er ikke strengt tatt en duplikat. Jeg er nysgjerrig på om VCF-validatoren også klager over disse sakene?

Det er virkelig atypisk, men dette spørsmålet ble bedt om fordi jeg faktisk opplevde dette i naturen. Jeg genererte ikke vcf, den ble sendt til meg, men overskriften antyder at den ble produsert av `freebayes` og sammenslåing av to separate filer. Så dette var sannsynligvis en gjenstand for sammenslåingen. Dessverre må jeg håndtere VCF-filer som blir gitt til meg av klienter, så selv om jeg kan insistere på at de samsvarer med standardene, er dette [_does_ conform] (https://bioinformatics.stackexchange.com/q/7105/298 ) (AFAIK), så jeg trengte en måte å fikse det på.
Den var bra! Dessverre (og beklager, jeg burde ha gjort dette klart) inneholder den faktiske filen jeg hadde multi-alleliske varianter der bare en var en dupe (se oppdatert spørsmål). Din tilnærming vil ikke fange dem, men bcftools gjør det (se svaret mitt). Dette fungerer bra for singler skjønt.
Dude, det er noe knasende VCF.
Velkommen til min verden :(
Fortsatt en kul løsning @DanielStandage!


ⓘ
Denne spørsmålet ble automatisk oversatt fra engelsk.Det opprinnelige innholdet er tilgjengelig på stackexchange, som vi takker for cc by-sa 4.0-lisensen den distribueres under.
Loading...