How do I use vcftools to filter by variant ID OR variant position?

Question

[jwhite007]

I'm having a problem obtaining a recoded vcf which contains variants with their IDs in a "snp file" OR their positions in a "position file." The following gives only variants with their IDs in the "snp file" AND their positions in the "position file."

vcftools --vcf <input vcf file> --snps <variant IDs file> --positions <variant positions file> --recode --out <output vcf file>

Is there any way that I can obtain that which I desire without having to do each separately, concatenating, and then getting rid of duplicates?

Thanks in advance.

[ATpoint]

Please use the format bar to indicate code. I did it for you this time.

Answer 1

[Pierre Lindenbaum]

not vcftools

using vcffilterjdk: http://lindenb.github.io/jvarkit/VcfFilterJdk.html

echo "rs12345" > ids.txt 
echo -e "1\t916834" > pos.txt
java -jar dist/vcffilterjdk.jar --body -f biostar347173.code src/test/resources/test_vcf01.vcf

(...)
1	916834	.	G	A	589	PASS	AA=A;AC=5;AF=0.625;AN=8;ANN=A|upstream_gene_variant|MODIFIER|C1orf170|ENSG00000187642|transcript|ENST00000479361|retained_intron||n.-4813C>T|||||4813|,A|intron_variant|MODIFIER|C1orf170|ENSG00000187642|transcript|ENST00000433179|protein_coding|1/4|c.54-281C>T||||||,A|intron_variant|MODIFIER|C1orf170|ENSG00000187642|transcript|ENST00000341290|protein_coding|1/4|c.-7-281C>T||||||	GT:AD:DP:DPF:GQ:GQX	./.:.:.:.:.:.	./.:.:.:.:.:.	0/1:18,17:35:1:99:42	0/1:19,14:33:0:99:35	1/1:0,41:41:1:99:36	0/1:22,18:40:4:99:40
1	935222	rs12345	C	A	380	PASS	AA=C;AC=7;AF=0.7;AN=10;ANN=A|missense_variant|MODERATE|HES4|ENSG00000188290|transcript|ENST00000428771|protein_coding|1/3|c.132G>T|p.Arg44Ser|331/1040|132/744|44/247||,A|upstream_gene_variant|MODIFIER|RP11-54O7.17|ENSG00000272512|transcript|ENST00000606034|lincRNA||n.-1791G>T|||||1791|,A|intron_variant|MODIFIER|HES4|ENSG00000188290|transcript|ENST00000304952|protein_coding|1/3|c.108+24G>T||||||,A|intron_variant|MODIFIER|HES4|ENSG00000188290|transcript|ENST00000484667|protein_coding|1/2|c.108+24G>T||||||,A|non_coding_transcript_exon_variant|MODIFIER|HES4|ENSG00000188290|transcript|ENST00000481869|retained_intron|1/2|n.255G>T||||||;CSQR=1|ENSR00000528874|regulatory_region_variant	GT:AD:DP:DPF:GQ:GQX	./.:.:.:.:.:.	1/1:0,22:22:8:63:26	0/1:15,14:29:0:99:36	0/1:16,16:32:5:99:35	1/1:0,25:25:7:72:26	0/1:21,15:36:2:99:35

view raw README.md hosted with ❤ by GitHub

	final Map<String,Set<Integer>> contig2pos = new HashMap<>();
	final Set<String> rsSet = new HashSet<>();
	public Object apply(final VariantContext vc) {
	/** first time, load the data */
	if(rsSet.isEmpty()) {
	try {
	rsSet.addAll(IOUtil.slurpLines(new java.io.File("ids.txt")));
	IOUtil.slurpLines(new java.io.File("pos.txt")).stream().forEach(L->{
	final int tab = L.indexOf('\t');
	final String c= L.substring(0,tab);
	final int p = Integer.parseInt(L.substring(tab+1));
	Set<Integer> pos = contig2pos.get(c);
	if(pos==null) {
	pos = new HashSet<>();
	contig2pos.put(c,pos);
	}
	pos.add(p);
	});
	}
	catch(Exception err)
	{
	throw new RuntimeException(err);
	}
	}
	if(vc.hasID() && this.rsSet.contains(vc.getID())) return true;
	Set<Integer> pos = contig2pos.get(vc.getContig());
	if(pos==null) return false;
	return pos.contains(vc.getStart());
	}

view raw biostar347173.code hosted with ❤ by GitHub

[jwhite007]

Pierre,

I just realized that you are the creator of jvarkit. Thank you for your hard work in getting that code going and making it freely available and open source.

Cheers, James

[jwhite007]

Pierre,

I'd definitely like this in an easy-to-use tool like vcftools or bcftools, but what you've provided seems to work. Thank you for posting this. I will interrogate further later, however...

vcftools --vcf <input vcf file> --snps <variant IDs file> --recode --out <output vcf file>

gives 1695 variants

vcftools --vcf <input vcf file> --positions <variant positions file> --recode --out <output vcf file>

gives 483 variants

vcftools --vcf <input vcf file> --snps <variant IDs file> --positions <variant positions file> --recode --out <output vcf file>

gives 60 variants

java -jar dist/vcffilterjdk.jar --body -f biostar347173.code <input vcf file> > <output vcf file>

gives 2118 variants

1695 + 483 - 60 = 2118 = Yay!!

[Pierre Lindenbaum]

If an answer was helpful you should upvote it, if the answer resolved your question you should mark it as accepted.