this:

  bedtools intersect -a <myvcf>.vcf.gz -b <myinterval>.bed -wa | \ java -Xmx10g -jar snpSift.jar filtr --set <myrsid>.txt "ID v SET [0]"  

lze nahradit jedním GATK SelectVariants https: // softwarem .broadinstitute.org / gatk / documentation / tooldocs / current / org_broadinstitute_gatk_tools_walkers_variantutils_SelectVariants.php

  java -jar gatk.jar -T SelectVariants -R ref.fa \ -L myinterval.bed --keepIDs myrsid.txt -V myvcf.vcf.gz -o out.vcf  

a navíc můžete také pracovat paralelně rozdělením souboru v posteli (např. chromozomy) a sloučit výsledky.

Obecně s BEDOPS:

  $ vcf2bed < < (gunzip -c snps.vcf) | bedops -e 1 - myRegions.bed > answer.bed  

Nebo:

  $ vcf2bed < < (gunzip -c snps.vcf) | bedmap --echo --echo-map-id --delim '\ t' myRegions.bed - > answer.bed  

atd.

BEDOPS podporuje vstupní toky ; použijte je, kde můžete.

To znamená, že pokud se chystáte často dotazovat na data VCF, převeďte je jednou a použijte je v následných operacích.

  $ vcf2bed < < (gunzip -c snps.vcf) > snps.bed  

Potom:

  $ bedops -e 1 snps.bed myRegions.bed > odpověď. bed  

atd.

Pokud máte výpočetní klastr, můžete dále paralelizovat práci rozdělením SNP na chromozomy pomocí bedextract :

  $ pro chr v `bedextract --list-chr snps.bed`; do echo $ {chr}; bedextract $ {chr} snps.bed > snps. $ {chr} .bed; hotovo  

Poté odešlete dávkovou úlohu pro každý chromozom:

  $ bedmap --chrom $ {chr} --echo --echo-map-id --delim '\ t' myRegions.bed snps. $ {chr} .bed > odpověď. $ {chr} .bed  

Použití bedextract a přidání - chrom $ {chr} k vašemu volání bedmap přinese velké zlepšení rychlosti oproti jiným sadám nástrojů, protože tyto možnosti se zaměřují na zajímavý chromozom a nic jiného.

Těžil jsem s příkazem a dokázal jsem dokončit provádění invertováním pořadí filtrování intervalů rsid a bed. Příkaz je následující:

  gzcat <myvcf>.vcf.gz | \ java -Xmx10g -jar filtr snpSift.jar --set <myrsid>.txt "ID v sadě [0]" | \ java -Xmx10g -jar snpSift.jar intervaly <mybed>.bed  

Možná je moje podmnožina rsid menší než moje BED intervaly, a proto je proces rychlejší a efektivnější.

Každopádně všem děkuji

Tady je další způsob, jak to udělat v roce 2019:

  $ bcftools view -T regions.bed -i 'ID=@<myrsid>.txt' input.vcf > output.vcf  

Věci lze zrychlit, pokud zkomprimujeme a indexujeme soubor input.vcf . bcftools pak bude mít téměř okamžitý náhodný přístup k pozici.

   $ bgzip -c input.vcf > vstup. vcf.gz $  tabix input.vcf.gz $ bcftools zobrazit -T regiony.bed -i 'ID=@<myrsid>.txt' input.vcf.gz > výstup.vcf  

Interně bcftools pracuje s bcf . Kdykoli je vstup a výstup vcf , převede vcf>bcf> vcf . Při práci s velkými soubory to vytváří velkou režii. První převod vcf na bcf může snížit dobu běhu pro další následné analýzy dramaticky. Podle mých zkušeností 2-3krát rychlejší.

   $ bcftools view -Ob input.vcf > input.bcf $  bcftools index input.bcf $ bcftools view -Ob -T regions.bed -i 'ID=@<myrsid>.txt' input.bcf > output.bcf