Pitanje:
Alat za pretraživanje za izračunavanje stope pogreške faze / sklopke
Wouter De Coster
2019-03-13 02:53:51 UTC
view on stackexchange narkive permalink

Tražim alat koji, s obzirom na vcf datoteku istine i testnu vcf datoteku, izračunava stopu pogreške faze / prekidača. Izvršio sam fazno oblikovanje vcf-a pomoću WhatsHap i želim usporediti ishod s nekim osnovnim istinskim faznim vcf-om koji imam. Ne mogu pronaći alat i ne želim ga sam napisati. Čini se kao prilično uobičajena stvar ... pa bih očekivao da postoji, ali moj je google-fu danas slab.

Primjer testnih linija: 

  #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT samplechr10 100588. G A. PASS AF = 0,18; AC = 1; NS = 150; AN = 2; EAS_AF = 0,14; EUR_AF = 0,39; AFR_AF = 0,08; AMR_AF = 0,31; SAS_AF = 0,2; VT = SNP; DP = 11907 GT: PS 1 | 0 : 48005chr10 102385. A G. PASS AF = 0,07; AC = 1; NS = 150; AN = 2; EAS_AF = 0,02; EUR_AF = 0,0; AFR_AF = 0,1; AMR_AF = 0,19; SAS_AF = 0,01; VT = SNP; DP = 22665 GT: PS 0 | 1 : 48005chr10 105170. G A. PASS AF = 0,03; AC = 1; NS = 150; AN = 2; EAS_AF = 0,02; EUR_AF = 0,0; AFR_AF = 0,0; AMR_AF = 0,19; SAS_AF = 0,0; VT = SNP; DP = 17140 GT: PS 0 | 1 : 48005chr10 105365. G T. PASS AF = 0,19; AC = 1; NS = 150; AN = 2; EAS_AF = 0,14; EUR_AF = 0,39; AFR_AF = 0,09; AMR_AF = 0,31; SAS_AF = 0,2; VT = SNP; DP = 21845 GT: PS 1 | 0 : 48005chr10 106057. T C. PASS AF = 0,13; AC = 1; NS = 150; AN = 2; EAS_AF = 0,02; EUR_AF = 0,0; AFR_AF = 0,32; AMR_AF = 0,19; SAS_AF = 0,01; VT = SNP; DP = 22559 GT: PS 0 | 1 : 48005

Primjer linija istine 

  #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT samplechr10 100554. C CT, CTT 30. . GT: AD 1 | 2: 0,1,1 chr10 100588. G A 30. . GT: AD 0 | 1: 1,1chr10 102385. A G 30. . GT: AD 1 | 0: 1,1chr10 102636. T C 30. . GT: AD 1 | 1: 0,2chr10 102757. T C 30. . GT: AD 0 | 1: 1,1
chr10 104815. OUU G, OZ 30. . GT: AD 2 | 1: 0,1,1 chr10 105170. G A 30. . GT: AD 1 | 0: 1,1chr10 105365. G T 30. . GT: AD 0 | 1: 1,1chr10 106057. T C 30. . GT: AD 1 | 0: 1,1chr10 106110. C G 30. . GT: AD 1 | 1: 0,2chr10 106261. A G 30. . GT: AD 1 | 0: 1,1chr10 108612. T G 30. . GT: AD 1 | 0: 1,1chr10 108646. A U 30. . GT: AD 0 | 1: 1,1chr10 108834. G A 30. . GT: AD 1 | 0: 1,1chr10 110840. C T 30. . GT: AD 1 | 1: 0,2chr10 111743. A U 30. . GT: AD 1 | 0: 1,1chr10 112016. T A 30. . GT: AD 1 | 0: 1,1chr10 112262. C T 30. . GT: AD 1 | 1: 0,2chr10 113006. C T 30. . GT: AD 0 | 1: 1,1chr10 113031. G A 30. . GT: AD 0 | 1: 1,1chr10 113136. G A 30. . GT: AD 1 | 0: 1,1chr10 113359. A T 30. . GT: AD 1 | 1: 0,2chr10 113583. C G 30. . GT: AD 0 | 1: 1,1chr10 113997. G C 30. . GT: AD 1 | 1: 0,2
Možete li nam dati nekoliko redaka svake datoteke za testiranje? Želite li samo usporediti polje genotipa za istu varijantu u dvije vcf datoteke?
Uredio sam svoje pitanje tako da uključuje primjere redaka. Datoteka istine ima fazne faze (kreirana je na osnovi sintetske diploidne metode https://www.nature.com/articles/s41592-018-0054-7 i https://github.com/lh3/CHM-eval / stablo / master / dip-call).
Jedan odgovor:
terdon
2019-03-13 18:32:41 UTC
view on stackexchange narkive permalink

Pod pretpostavkom da su vaše datoteke onakve kakve prikazujete, a polje genotipa je uvijek prvi unos uzorka informacijskog polja (10.), onda to možete učiniti izravno jednostavnom skriptom awk: 

   $ awk -F "\ t" -vOFS = "\ t" 'POČINITE {ispis "#CHROM", "POS", "ID", "REF", " ALT "," Istina "," Test "} {if (/ ^ [^ #] /) {var = $ 1 $ 2 $ 4 $ 5; gt = substr (10,1,3 USD); ako je (NR == FNR) {a [var] = gt; sljedeći;} if (var u && a [var]! = gt) {print $ 1, $ 2, $ 3, $ 4, $  5, gt, a [var]}}}'uth.vcf test.vcf #CHROM POS ID REF ALT Test istinechr10 100588. G A 1 | 0 0 | 1chr10 102385. A G 0 | 1 1 | 0chr10 105170. G A 0 | 1 1 | 0chr10 105365. G T 1 | 0 0 | 1chr10 106057. TC 0 | 1 1 | 0

Evo skripte podijeljene u zasebne retke radi jasnosti (i dalje možete kopirati zalijepiti izravno u terminal): 

  awk -F "\ t" -vOFS = "\ t" 'POČINITE {ispis "#CHROM", "POS", "ID", "REF", "ALT", "Istina", "Test"} {if (/ ^ [^ #] /) {var =  $ 1 $  2  $ 4 $  5; gt = substr ( 10,1,3 USD); ako je (NR == FNR) {a [var] = gt; Sljedeći; } if (var u && a [var]! = gt) {print $ 1, $ 2, $ 3, $ 4, $  5, gt, [var]}}}'uth.vcf test.vcf

Osnovna ideja je:

  • awk -F "\ t" -vOFS = "\ t" : postavite ulaz i izlaz separatora polja na karticu.
  • BEGIN {print "#CHROM ...} : ispis zaglavlja.
  • if (/ ^ [^ #] /) {: preskočite retke zaglavlja ulaznih vcf datoteka
  • var = $ 1 $ 2 $ 4 $ 5 : dobijte "potpis" za ovu varijantu , spajanjem polja kromosoma, položaja, ref i alt.
  • gt = substr ($ 10,1,3); : dobivanje vrijednosti genotipa.
  • if (NR == FNR) {: ako je ovo prva datoteka /
  • a [var] = gt; sljedeća : spremite genotip za ovu varijantu i preskočite na sljedeći redak.
  • if (var u && a [var]! = gt) {: ovo dio će se pokrenuti samo za drugu datoteku zbog "sljedećeg" u prethodnom bloku. Ako je ova varijanta bila u skupu istine i ako njezin genotip nije isti kao u testnom skupu, ispišite varijantu i dva genotipa.

Ako trebate stope pogrešaka, možete tako proširiti skriptu: 

  awk -F "\ t" -vOFS = "\ t" 'BEGIN {print "#CHROM", "POS", "ID" , "REF", "ALT", "Istina", "Test"} {if (/ ^ [^ #] /) {var =  $ 1 $  2  $ 4 $  5; gt = substr ( 10,1,3 USD); ako je (NR == FNR) {a [var] = gt; Sljedeći; } tot ++; if (var u && a [var]! = gt) {bad ++; ispiši $ 1, $ 2, $ 3, $ 4, $  5, gt, a [var]}}} END {printf "Točno:% s od% s (% s%) \ n", tot-bad, tot , (ukupno-loše) * 100 / ukupno; printf "Pogrešno:% s od% s (% s%) \ n", loše, malo, loše * 100 / tot;}'uth.vcf test.vcf

Na vašem primjeru skup podataka koji se ispisuje: 

  #CHROM POS ID REF ALT Truth Testchr10 100588. G A 1 | 0 0 | 1chr10 102385. A G 0 | 1 1 | 0chr10 105170. G A 0 | 1 1 | 0chr10 105365. G T 1 | 0 0 | 1chr10 106057. T C 0 | 1 1 | 0Točno: 0 od 5 (0%) Pogrešno: 5 od 5 (100%)
U testnom VCF, "PS" znači "fazni skup". Prekidači između skupova faza ne računaju se kao pogreške. Općenito, trebamo odgovarajuću skriptu za ocjenu. Ne postoji uobičajena skripta vjerojatno zato što se ulazi često puno mijenjaju.
@user172818 oh vau, hvala! Prije nisam sreo `PS`. Znači li to da je moj pojednostavljeni pristup u osnovi beskoristan?
Mislim da će biti korisno kada "PS" nije uključen. Na primjer, Hi-C faziranje treba dati fazu cijelog kromosoma. Ne treba nam "PS". Vaša je skripta primjenjiva u tom slučaju.
Da, dok skup istine ima "lijevi" i "desni" haplotip (očinski ili majčin), testni skup je podijeljen u haplotipske / fazne blokove. Unutar bloka možemo shvatiti koji su SNP-ovi na istom haplotipu (ili nisu), ali ne među blokovima. Moj primjer nije savršen, ali vjerujem da je ovdje 100% ispravno postavljen.


ⓘ
Ova pitanja su automatski prevedena s engleskog jezika.Izvorni sadržaj dostupan je na stackexchange-u, što zahvaljujemo na cc by-sa 4.0 licenci pod kojom se distribuira.
Loading...