Spørsmål:
Velge nettsteder fra VCF som har et alt AD> 10
Matt Bashton
2017-07-03 19:41:40 UTC
view on stackexchange narkive permalink

Jeg har høydybde varianter som er opprettet ved hjelp av HaplotypeCaller med --output_mode EMIT_ALL_SITES Jeg er interessert i å finne alle nettsteder (uansett genotype kaller heterozygot eller homozygot) hvor minst en av alternative alleler har en AD -verdi (Allelic Depth) større enn 10, I . e . støttes av mer enn 10 lesninger. Også ideelt vil jeg ha tilbake mer enn bare den første alternative allelen. Merk at jeg ikke vil ha tilbake linjer med VCF hvis vi bare ser et AD-antall for ref-allelen.

Så i eksemplet på VCF-utdraget nedenfor ønsker jeg å velge linjer: 6,7, 8,12,13 og 14, som har GT: AD-verdier 1/1: 1,988: 989 0/1: 116,92 0/1: 220,234 henholdsvis 0/1: 62,611 1/1: 0,109 . 

  #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 12908_DIAG3 187446740. T. Evighet    . AN = 2; DP = 1095; MQ = 60,00 GT: AD: DP 0/0: 1095: 10953 187446741. C. Evighet    . AN = 2; DP = 1117; MQ = 60,00 GT: AD: DP 0/0: 1117: 11173 187446752. A. Evighet    . AN = 2; DP = 1297; MQ = 60,00 GT: AD: DP 0/0: 1297: 12973 187446763. C. Evighet    . AN = 2; DP = 1494; MQ = 60,00 GT: AD: DP 0/0: 1494: 14943 187451574. C. Evighet    . AN = 2; DP = 1493; MQ = 60,00 GT: AD: DP 0/0: 1493: 14933 187451609 rs1880101 A G 39794,03. AC = 2; AF = 1,00; AN = 2; BaseQRankSum = 1,859; ClippingRankSum = 0,000; DB; DP = 995; ExcessHet = 3,0103; FS = 0,000; MLEAC = 2; MLEAF = 1,00; MQ = 60,00; MQRankSum = 0,000; QD = 24,56; ReadPosRankSum = 0,406; SOR = 8,234 GT: AD: DP: GQ: PL 1/1: 1,988: 989: 99: 39808,2949,04 1803279. T G 0. AC = 0; AF = 0,00; AN = 2; BaseQRankSum = -6,652; ClippingRankSum = 0,000; DP = 245; ExcessHet = 3,0103; FS = 89,753; MLEAC = 0; MLEAF = 0,00; MQ = 59,97; MQRankSum = 0,000; ReadPosRankSum = -2.523; SOR = 6.357 GT: AD: DP: GQ: PL 0/0: 211,23: 234: 99: 0,364,6739
4 1803307 rs2305183 T C 2486.60. AC = 1; AF = 0,500; AN = 2; BaseQRankSum = -5,049; ClippingRankSum = 0,000; DB; DP = 215; ExcessHet = 3,0103; FS = 1,110; MLEAC = 1; MLEAF = 0,500; MQ = 59,97; MQRankSum = 0,000 ; QD = 11,95; ReadPosRankSum = -0,045; SOR = 0,809 GT: AD: DP: GQ: PL 0/1: 116,92: 208: 99: 2494,0,36734 1803671. C A 0. AC = 0; AF = 0,00; AN = 2; BaseQRankSum = -0,880; ClippingRankSum = 0,000; DP = 450; ExcessHet = 3,0103; FS = 0,000; MLEAC = 0; MLEAF = 0,00; MQ = 60,00; MQRankSum = 0,000; ReadPosRankSum = -0,953; SOR = 0,572 GT: AD: DP: GQ: PL 0/0: 445,2: 447: 99: 0,1272,159584 1803681. T C 0. AC = 0; AF = 0,00; AN = 2; BaseQRankSum = -1,654; ClippingRankSum = 0,000; DP = 483; ExcessHet = 3,0103; FS = 0,000; MLEAC = 0; MLEAF = 0,00; MQ = 60,00; MQRankSum = 0,000; ReadPosRankSum = -0,422; SOR = 0,664 GT: AD: DP: GQ: PL 0/0: 479,2: 481: 99: 0,1408,185384 1803703. A G 0. AC = 0; AF = 0,00; AN = 2; BaseQRankSum = -1,704; ClippingRankSum = 0,000; DP = 458; ExcessHet = 3,0103; FS = 0,000; MLEAC = 0; MLEAF = 0,00; MQ = 60,00; MQRankSum = 0,000; ReadPosRankSum = 0.299; SOR = 0.497 GT: AD: DP: GQ: PL 0/0: 454,2: 456: 99: 0,1325,180954 1803704 rs2234909 TC 6676,60. AC = 1; AF = 0,500; AN = 2; BaseQRankSum = -2,605; ClippingRankSum = 0,000; DB; DP = 456; ExcessHet = 3,0103; FS = 1,753; MLEAC = 1; MLEAF = 0,500; MQ = 60,00; MQRankSum = 0,000 ; QD = 14,71; ReadPosRankSum = 0,324; SOR = 0,849 GT: AD: DP: GQ: PL 0/1: 220,234: 454: 99: 6684,0,63664 1803824 rs2305184 CG 2030,60. AC = 1; AF = 0,500; AN = 2; BaseQRankSum = 8,083; ClippingRankSum = 0,000; DB; DP = 124; ExcessHet = 3,0103; FS = 6,128; MLEAC = 1; MLEAF = 0,500; MQ = 60,00; MQRankSum = 0,000; QD = 16,51; ReadPosRankSum = 0,180; SOR = 0,096 GT: AD: DP: GQ: PL 0/1: 62,61: 123: 99: 2038,0,17664 1805296 rs3135883 GA 3876.03. AC = 2; AF = 1,00; AN = 2; DB; DP = 110; ExcessHet = 3,0103; FS = 0,000; MLEAC = 2; MLEAF = 1,00; MQ = 60,00; QD = 29,22; SOR = 9,401 GT: AD: DP : GQ: PL 1/1: 0,109: 109: 99: 3890,326,0

En dropbox-lenke for fil

Jeg hadde i utgangspunktet vurdert å bruke GATKs SelectVariants, men jeg er ikke sikker på at JEXL har muligheten til å velge ut hva jeg vil spesifikt annet enn et teppe AD> 10 som vil gi meg både ref- og alt-alleler med AD> 10. Kanskje det er en bioawk-løsning eller noe mer forseggjort med coreutils som vellykket kan returnere nettsteder med et alt AD-antall> 10?

krysspost: http://gatkforums.broadinstitute.org/gatk/discussion/comment/40003#Comment_40003
Tre svar:
Pierre
2017-07-03 20:41:14 UTC
view on stackexchange narkive permalink

ved hjelp av vcfilterjs

og følgende skript: 

  function accept (vc) {var i, j; for (i = 0; i< vc.getNSamples (); ++ i) {var genotype = vc.getGenotype (i); hvis (! genotype.hasAD ()) fortsetter; var ad = genotype.getAD (); / * loop over AD som starter fra '1' = første ALT * / for (j = 1; j< ad.length; ++ j) {if (ad [j] >10) return true; }} returner falsk; } godta (variant);

bruk: 

  java -jar dist / vcffilterjs.jar -f script.js Test.vcf
Jeg tror dette systemet du har laget ser ut til å være veldig fleksibelt. Jeg vil være veldig interessert i å se noen flere gjennomganger og dokumenterte eksempler på Github fordi jeg er sikker på at det kan løse mange komplekse filtreringsproblemer med variasjoner.
@MatthewBashton takk! og det er en mye raskere alternativ versjon siden i dag (ikke javascript, men 100% java basert) http://lindenb.github.io/jvarkit/VcfFilterJdk.html
@DanielStandage ikke seler, bare lambdas bruker ** vcffilterjdk **: `java -jar dist / vcffilterjdk.jar -e 'returvariant.getGenotypes (). Stream (). Filter (G-> G.hasAD () && java.util .Arrays.stream (G.getAD ()). Hopp over (1) .filter (AD-> AD> 10) .findAny (). IsPresent ()). FindAny (). IsPresent (); ' input.vcf `
Hyggelig jeg må også gi jdk-versjonen en sjanse. Jeg har ingen problemer med seler eller innrykkstil ¯ \\ _ (ツ) _ / ¯
Cotton Seed
2017-07-10 00:15:02 UTC
view on stackexchange narkive permalink

Du kan gjøre dette i Hail: 

  fra haglimport * hc = HailContext () (hc.import_vcf ('test.vcf') .filter_variants_expr ( 'gs.eksisterer (g = > g.ad [1:]. eksisterer (d = > d > 10))') .export_vcf ('filtered.vcf'))

Dette fungerer med et hvilket som helst antall prøver og vil beholde varianter der minst en prøve har en genotype med en alternativ allelstøtte med mer enn 10 lesninger. 

  >>> hc.import_vcf ('filtered.vcf'). count () (1L, 6L)

count returnerer antall prøver (1) og antall varianter (6).

Ta en titt på start-siden eller tutorials hvis du vil prøve det!

winni2k
2017-07-04 21:24:14 UTC
view on stackexchange narkive permalink

Dette fungerer nå med utviklingsversjonen av Bcftools v1.5 (commit 4f134df). Takk til Petr Danecek for å legge til funksjonen. Jeg forventer at denne funksjonen vil komme seg inn i neste utgivelse av Bcftools: 

  git clone git: //github.com/samtools/htslib.gitgit clone git: //github.com/samtools /bcftools.git(cd bcftools; make) bgzip Test.vcf./bcftools/bcftools index Test.vcf.gz./bcftools/bcftools filter -i 'AD [1-] > 10' Test.vcf.gz

Output uten header (jeg har endret den andre linjen for å være tri-allel for å demonstrere filtreringsarbeidet): 

  3 187451609 rs1880101 AG 39794 PASS AC = 2; AF = 1; AN = 2; BaseQRankSum = 1.859; ClippingRankSum = 0; DB; DP = 995; OverskuddHet = 3.0103; FS = 0; MLEAC = 2; MLEAF = 1; MQ = 60; MQRankSum = 0; QD = 24.56; ReadPosRankSum = 0.406; SOR = 8.234 GT: AD: DP: GQ: PL 1/1: 1.988: 989: 99: 39808,2949,04 1803279. TG 0 PASS AC = 0; AF = 0; AN = 2; BaseQRankSum = -6.652; ClippingRankSum = 0; DP = 245; ExcessHet = 3.0103; FS = 89.753; MLEAC = 0; MLEAF = 0; MQ = 59.97; MQRankSum = 0; ReadPosRankSum = -2.523; SOR = 6.357 GT: AD: DP: GQ: PL 0/0/0: 211,3,34: 234: 99: 0,364,67394 1803307 rs2305183 TC 2486,6 PASS AC = 1; AF = 0,5 ; AN = 2; BaseQRankSum = -5.049; ClippingRankSum = 0; DB; DP = 215; ExcessHet = 3.0103; FS = 1.11; MLEAC = 1; MLEAF = 0.5; MQ = 59.97; MQRankSum = 0; QD = 11.95; ReadPosRankSum = -0,045; SOR = 0,809 GT: AD: DP: GQ: PL 0/1: 116,92: 208: 99: 2494,0,36734 1803704 rs2234909 TC 6676,6 PASS AC = 1; AF = 0,5; AN = 2; BaseQRankSum = -2,605; ClippingRankSum = 0; DB; DP = 456; ExcessHet = 3,0103; FS = 1,753; MLEAC = 1; MLEAF = 0,5; MQ = 60; MQRankSum = 0; QD = 14,71; ReadPosRankSum = 0,324; SOR = 0,849 GT : AD: DP: GQ: PL 0/1: 220,234: 454: 99: 6684,0,63664 1803824 rs2305184 CG 2030,6 PASS AC = 1; AF = 0,5; AN = 2; BaseQRankSum = 8,083; ClippingRankSum = 0; DB; DP = 124; ExcessHet = 3.0103; FS = 6.128; MLEAC = 1; MLEAF = 0.5; MQ = 60; MQRankSum = 0; QD = 16.51; ReadPosRankSum = 0.18; SOR = 0.096 GT: AD: DP: GQ: PL 0 / 1: 62,61: 123: 99: 2038,0,1766
4 1805296 rs3135883 GA 3876.03 PASS AC = 2; AF = 1; AN = 2; DB; DP = 110; ExcessHet = 3.0103; FS = 0; MLEAC = 2; MLEAF = 1; MQ = 60; QD = 29.22; SOR = 9.401 GT: AD: DP: GQ: PL 1/1: 0,109: 109: 99: 3890,326,0


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