Küsimus:
VCF-ist saitide valimine, mille alt AD> 10
Matt Bashton
2017-07-03 19:41:40 UTC
view on stackexchange narkive permalink

Mul on HaplotypeCalleri abil loodud ulatuslik teisendkõne funktsiooniga --output_mode EMIT_ALL_SITES . Olen huvitatud kõigi saitide leidmisest (olenemata genotüübi kõne heterosügootsest või homosügootsest), kus on vähemalt üks alternatiivsete alleelide AD väärtus (alleelisügavus) on suurem kui 10, I . e . toetab rohkem kui 10 lugemist. Ideaalis soovin ka tagasi rohkem kui lihtsalt esimest alternatiivset alleeli. Pange tähele, et ma ei taha VCF-i tagumisi ridu, kui näeme ainult ref-alleeli AD-arvu.

Nii et allpool toodud näites VCF-i koodilõigu soovin valida read: 6,7, 8,12,13 ja 14, millel on GT: AD väärtused 1/1: 1,988: 989 0/1: 116,92 0/1: 220,234 0/1: 62,611 1/1: 0,109 . 

  #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 12908_DIAG3 187446740. T. Lõpmatus. AN = 2; DP = 1095; MQ = 60,00 GT: AD: DP 0/0: 1095: 10953 187446741. C. Lõpmatus. AN = 2; DP = 1117; MQ = 60,00 GT: AD: DP 0/0: 1117: 11173 187446752. A. Lõpmatus. AN = 2; DP = 1297; MQ = 60,00 GT: AD: DP 0/0: 1297: 12973 187446763. C. Lõpmatus. AN = 2; DP = 1494; MQ = 60,00 GT: AD: DP 0/0: 1494: 14943 187451574. C. Lõpmatus. 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

Faili dropbox-link

Kaalusin algselt GATKi SelectVariantsi kasutamist, kuid ma pole kindel, et JEXL-il on võimalus välja valida, mida ma konkreetselt soovin, välja arvatud tekk AD> 10, mis annab mulle nii viite- kui ka alleelialused AD-ga> 10. Võib-olla Kas on olemas biowk lahendus või midagi muud keerukamat koos coreutilsiga, mis võiks edukalt tagastada saidid, mille alt AD arv on> 10?

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

kasutades vcfilterjs

ja järgmist skripti: 

  function accept (vc) {var i, j; jaoks (i = 0; i< vc.getNSamples (); ++ i) {var genotüüp = vc.getGenotype (i); kui (! genotüüp.hasAD ()) jätkub; var ad = genotüüp.getAD (); / * silmus AD kohal alates '1' = esimene ALT * / for (j = 1; j< ad.length; ++ j) {if (ad [j] >10) tagastab true }} return false; } accept (variant);

kasutamine: 

  java -jar dist / vcffilterjs.jar -f script.js Test.vcf
Ma arvan, et see teie loodud süsteem näib olevat väga paindlik. Ma oleksin väga huvitatud Githubis veel mõningate ülevaadete ja dokumenteeritud näidete nägemisest, sest olen kindel, et see suudab lahendada palju keerukate variantide filtreerimise probleeme.
@MatthewBashton aitäh! ja alates tänasest on palju kiirem alternatiivne versioon (mitte javascript, vaid 100% Java-põhine) http://lindenb.github.io/jvarkit/VcfFilterJdk.html
@DanielStandage mitte klambrid, vaid lambdas, kasutades ** vcffilterjdk **: `java -jar dist / vcffilterjdk.jar -e 'tagasivariand.getGenotypes (). Stream (). Filter (G-> G.hasAD () && java.util .Arrays.stream (G.getAD ()). Skip (1) .filter (AD-> AD> 10) .findAny (). IsPresent ()). FindAny (). IsPresent (); " sisend.vcf `
Tore, et pean ka jdk-versiooni käima panema. Mul pole probleeme teie traksidega ega taandestiiliga ¯ \\ _ (ツ) _ / ¯
Cotton Seed
2017-07-10 00:15:02 UTC
view on stackexchange narkive permalink

Seda saate teha jaotises Tervitus: 

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

See töötab suvalise arvu proovidega ja säilitab variandid, kus vähemalt ühel proovil on alternatiivse alleelitoega genotüüp, üle 10 lugemise.

Kontrollimaks, kas saime oodatud 6 varianti: 

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

count tagastab näidiste arvu (1) ja variantide arvu (6).

Vaadake alustamise lehte või õpetusi, kui soovite proovida!

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

See töötab nüüd Bcftools v1.5 arendusversiooniga (toime 4f134df). Täname Petr Danecekit funktsiooni lisamise eest. Eeldan, et see funktsioon jõuab Bcftools järgmise versiooni juurde: 

  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

Väljund ilma päiseta (filtreerimistööde demonstreerimiseks muutsin teist rida kolmealleeliseks: 

  3 187451609 rs1880101 AG 39794 PASS AC = 2; AF = 1; AN = 2; BaseQRankSum = 1,859; ClippingRankSum = 0; DB; DP = 995; ExcessHet = 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: 1988: 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


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