Küsimus:
Kuidas eemaldada (mitte triviaalsed) duplikaadid VCF-failist?
terdon
2019-02-27 21:56:46 UTC
view on stackexchange narkive permalink

See on seotud küsimusega, mille esitasin siin. Mõelgem vcf-failile, mis sisaldab duplikaatvariante, kuid kus duplikaadid ei ole lihtsalt ühesuguses tähistuses sama asi, vaid üks on teise alamhulk. Näiteks:

  ## 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. ÜLE ANDMA    . GT 0 / 1chr12 529516. C T. ÜLE ANDMA    . GT 0/1  

Need kaks varianti on tegelikult samad. Nende tulemuseks on täpselt sama genotüüp. AACAC-i muutmine asendis 529514 AATAC-ks tähendab lihtsalt positsiooni C muutmist asendisse T asendis 529516.

Kas on olemas mõni tööriist, mis suudaks sellised duplikaadid tuvastada ja need eemaldada? Proovisin vcfuniq vcflib -st, kuid see ei tundu seda duplikaadina ära tundvat. Ma arvan, et see vaatleb ainult 4 esimest välja ja arvestab, et need duplikaadid muudavad täpselt samade väärtustega variante neljandal neljandal väljal:

  $ ./bin/vcfuniq test.vcf ## fileformat = VCFv4.1 ## reference = foo ## FORMAT = <ID = GT, Number = 1, Type = String, Description = "Genotype" > ## contig = <ID = chr12> # CHROM POS ID REF ALT QUAL FILTER14 INFO FORMAT S. AACAC AATAC. ÜLE ANDMA    . GT 0 / 1chr12 529516. C T. ÜLE ANDMA    . GT 0/1  

Nagu selgitatud lingitud küsimuses, peab EBI vcf_validator seda kehtetuks. Ja tegelikult pole mõtet neid duplikaate igal juhul omada, nii et kas saan kuidagi neid tuvastada ja eemaldada? Eelistatavalt olemasolev tööriist, kuid olen avatud ka skriptimislahendustele.


Sellised juhtumid muudavad seda keerulisemaks:

  ## fileformat = VCFv4 .1 ## reference = foo ## FORMAT = <ID = GT, arv = 1, tüüp = string, kirjeldus = "genotüüp" >
## 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. ÜLE ANDMA    . GT 0/1  

Kahjuks ei jõua see Danieli nutika skripti lähenemisviisiga:

  $ kassi test2.vcf | foo.py ## fileformat = VCFv4.1 ## reference = foo ## FORMAT = <ID = GT, arv = 1, tüüp = string, kirjeldus = "genotüüp" > ## contig = <ID = chr12AiliyghyaltID = chr12> # CHROM POS ID REF ALT KVALITEETFILTRI INFO VORMING Proov1chr12 529514 529514 AACAC AAT, AATAC 0,00. . GT 0 / 1chr12 529516 529516 C T. ÜLE ANDMA    . GT 0/1  
Kaks vastused:
terdon
2019-02-27 22:10:34 UTC
view on stackexchange narkive permalink

Selgub, et bcftools saab seda teha (testitud bcftools-1.8-ga), kui annate sellele viitegenoomi testimiseks:

  $ bcftools norm -d pole -f hg19.fa test.vcf ## fileformat = VCFv4.1 ## FILTER = <ID = PASS, Description = "Kõik filtrid on läbitud" > ## reference = foo ## FORMAT = <ID = GT, Number = 1 , Type = String, Description = "Genotype" > ## contig = <ID = chr12> ## bcftools_normVersion = 1.8 + htslib-1.8 ## bcftools_normCommand = norm -d pole -f hg19.fa test.vcf; Kuupäev = kolmapäev 27. veebruar 16:08:44 2019 # CHROM POS ID REF ALT QUAL FILTER INFO FORMAT Sample1chr12 529516. C T. ÜLE ANDMA    . GT 0 / 1Jooned kokku / jagatud / ümber joondatud / vahele jäetud: 2/0/1/0  

Mitme alleelilise variandi keerukama juhtumi jaoks küsimuse teises näites VCF saate seda kaks korda bcftoolsi kaudu käivitada. Kui olete mitme alleeli variandi vasakule joondamiseks ja jagamiseks ning duplikaatide eemaldamiseks kasutanud norm :

  $ bcftools norm -m -any -NO z - O v -o - ~ / test2.vcf | bcftools norm -d pole -f hg19.faRead kokku / jagatud / ümber joondatud / vahele jäetud: 2/1/0/0 ## fileformat = VCFv4.1 ## FILTER = <ID = PASS, Description = "Kõik filtrid on läbitud" > ## viide = foo ## FORMAT = <ID = GT, arv = 1, tüüp = string, kirjeldus = "genotüüp" > ## contig = <ID = chr12> ## bcftools_normVersion = 1.8 + htslib-1.8 ## bcfanyols_norm -NO z-O v -o - test2.vcf; Kuupäev = kolmapäev 27. veebruar 18:18:32 2019 ## bcftools_normCommand = norm -d pole -f hg19.fa -; Kuupäev = kolmapäev 27. veebruar 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 / 0Jooned kokku / jagatud / ümber joondatud / vahele jäetud: 03/02/0  
Suurepärane teave tagataskus. Tore lahendus!
Soovin, et saaksin sellele veel +1 lisada.
Daniel Standage
2019-02-27 22:16:00 UTC
view on stackexchange narkive permalink

Ma ei ole VCF-i ekspert (vähesed võivad öelda, et nad on!), kuid olen viimastel aastatel töötanud palju VCF-i andmetega, mõlemad tööriistad VCF-i tarbimiseks ja tootmiseks. Ma pole kunagi näinud selliselt kodeeritud variante ja see näib olevat mittekanooniline. Tavaliselt:

  • ühe nukleotiidi variandid (SNV) on kodeeritud ühe alusega REF alleelina ja ühe alusega ALT alleelina.
  • sisestuste või kustutamise korral on REF- ja ALT-alleelidest lühem üks alus, alus eelneb sisestatud / kustutatud järjestusele. Seega on REF- ja ALT-alleelide esimene alus alati sama.
  • Harvadel juhtudel, kui kaks või enam järjestikust asendust moodustavad multinukleotiidvariandi (MNV), on REF-i ja ALT-alleelide pikkus sama.

Ühepikkuste mitme bp-stringide kasutamine SNV-de kodeerimiseks on tarbetu ja, nagu te juba märkisite, problemaatiline. See paneb mind mõtlema, et see on viga või VCF-i tekitanud variandi ennustaja "funktsioon".

Sel juhul kirjutaksin väikese skripti, mis kontrolliks variante, kus REF ja ALT alleelid on sama pikk. Kui baas on REF-i ja ALT-i jaoks ühes ja samas asukohas sama, langetage see ja kohandage asukohta vastavalt sellele.

Allolev skript teisendab need funky SNV-d kanooniliseks esituseks ja töötab ka MNV-del. Seejärel peaksid duplikaatide eemaldamiseks töötama standardsed tööriistad.

  #! / Usr / bin / env python3def canonicalize (voogesitus): voosisesele reale: kui mitte line.startswith ('#'): väärtused = line.split ('\ t') pos = int (väärtused [1]) ref, alt = väärtused [3: 5] kui len (ref) > 1 ja len ( ref) == len (alt): # Mitu bp lõppu n, (r, a) jaoks kärpida loendis (zip (viide [:: - 1], alt [:: - 1])): kui r! = a: revoffset = -1 * n murda # Mitu bp eest ära kärpida
n, (r, a) loendamiseks (zip (ref, alt)): kui r! = a: nihe = n väärtused [1] = str (pos + nihe) väärtused [3] = ref [nihe: tagasimakse] väärtused [4] = alt [offset: revoffset] katkendjoon = '\ t'.join (väärtused) annavad lineif __nimi__ ==' __main__ ': impordi sys reale kanoniseerimisel (sys.stdin): print (rida, lõpp = '')  

UPDATE : edasisel järelemõtlemisel on teie loetletud keerulisem näide tegelikult mõttekas. Viites on meil järjestus AACAC ja alternatiivsed alleelid esindavad selles kahte variatsiooni: kahe viimase bp kustutamine (esimesel juhul) ja keskmise C punktmutatsioon kuni T (mõlemal juhul). Tavaliselt eelneb indeli määratlusele ainult üks bp, seega oleksin selle keerulise variandi kodeerinud järgmiselt: ref = ACAC alt = AT, ATAC .

Nii et SNV on "implicit by" / "kodeeritud keerulises variandis" / "kodeeritud, kuid see pole rangelt duplikaat. Mind huvitab, kas VCF-i valideerija kaebab ka nende juhtumite üle?

See on tõepoolest ebatüüpiline, kuid see küsimus esitati seetõttu, et tegelikult kohtasin seda looduses. Ma ei genereerinud vcf-d, see saadeti mulle, kuid päise põhjal võib öelda, et selle on tootnud `freebayes` ja kahe eraldi faili ühendamine. Nii et see oli tõenäoliselt ühinemise artefakt. Kahjuks pean tegelema klientide antud VCF-failidega, nii et kuigi saan nõuda, et need vastaksid standarditele, on see [_does_ vastavuses] (https://bioinformatics.stackexchange.com/q/7105/298 ) (AFAIK), nii et vajasin viisi selle parandamiseks.
Kena! Kahjuks (ja vabandust, ma oleksin pidanud selle selgeks tegema) sisaldab minu tegelik fail mitme alleeli variante, kus ainult üks oli dupe (vt uuendatud küsimust). Teie lähenemisviis ei taba neid, kuid bcftools siiski (vt minu vastust). See sobib suurepäraselt vallalistele.
Kutt, see on mingi räpane VCF.
Tere tulemast minu maailma :(
Ikka lahe lahendus @DanielStandage!


See küsimus ja vastus tõlgiti automaatselt inglise keelest.Algne sisu on saadaval stackexchange-is, mida täname cc by-sa 4.0-litsentsi eest, mille all seda levitatakse.
Loading...