Removing homozygous reference genotypes from multi-sample vcf file

Question: Removing homozygous reference genotypes from multi-sample vcf file

11 months ago by

seta • 1.3k

Sweden

seta • 1.3k wrote:

Hi everybody,

I want to remove all sites that are homozygous reference (0/0) per each row in the vcf file and get the output with vcf format. So, I simply used zgrep -v "0[/|]0" < file1.gz.vcf > output.vcf. However, it sounds that many variants is removed. Just for make sure, please kindly let me know if I did right?

Thanks

genotype removing vcf • 806 views

ADD COMMENT • link •

modified 11 months ago by finswimmer ♦ 13k • written 11 months ago by seta • 1.3k

11 months ago by

Pierre Lindenbaum ♦ 128k

France/Nantes/Institut du Thorax - INSERM UMR1087

Pierre Lindenbaum ♦ 128k wrote:

~~bcftools view --min-ac 1 input.vcf~~

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

java -jar src/jvarkit-git/dist/vcffilterjdk.jar -e 'return variant.getGenotypes().stream().noneMatch(G->G.isHomRef() || G.isNoCall());'  input.vcf.gz

ADD COMMENT • link modified 11 months ago • written 11 months ago by Pierre Lindenbaum ♦ 128k

Sorry, I forgot to mention the vcf file is multisample with about 1200 individuals. So, the above command didn't work. How I can get just hom-var and heterozygous genotypes per each row (variant)?

ADD REPLY • link written 11 months ago by seta • 1.3k

I updated my answer.

ADD REPLY • link written 11 months ago by Pierre Lindenbaum ♦ 128k

11 months ago by

finswimmer ♦ 13k

Germany

finswimmer ♦ 13k wrote:

I strongly recommend using existing tools for filtering vcf files.

With bcftools it works like this:

bcftools view -i 'GT[*]="alt"' input.vcf

This will give you any site where at least one sample has an alternate allel.

ADD COMMENT • link modified 11 months ago • written 11 months ago by finswimmer ♦ 13k

as far as I undersand OP don't want ANY HOM_REF sample in the row.

ADD REPLY • link written 11 months ago by Pierre Lindenbaum ♦ 128k

Hmm, yes might be.

In this case:

bcftools view -e 'GT[*]="RR"' input.vcf

ADD REPLY • link written 11 months ago by finswimmer ♦ 13k

Right, I just hetero and homo-var NOT homo-ref per row. I before used the bcftools view -e 'GT[*]="RR"' file.vcf, but sounds worked wrong as it just kept 530 variants from 2134618 variants. Our admin must install Pierre's tool on the cluster to test it.

ADD REPLY • link modified 11 months ago • written 11 months ago by seta • 1.3k

For me this does not sound surprisingly.

If you have no sample within 1200 samples, which is hom REF on a given position, this means, that what is called ref, is very rare. Of course this is possible, but I would guess it is rare.

ADD REPLY • link modified 11 months ago • written 11 months ago by finswimmer ♦ 13k

Our admin must install Pierre's tool on the cluster to test it

As long as you have java installed on the cluster you should be able to run Pierre's program from your own directory.

ADD REPLY • link modified 11 months ago • written 11 months ago by genomax ♦ 83k

11 months ago by

Dave Carlson • 320

Stony Brook University, NY

Dave Carlson • 320 wrote:

I like to use SnpSift for pulling out specific genotype combinations. In your case, I think you could do the following:

zcat file1.gz.vcf | java -jar SnpSift.jar filter "isHom( GEN[0] )  & isRef( GEN[0] ) " --inverse > output.vcf

This will pull out all lines in the vcf that are not homozygous reference. If you're trying to do this on a specific sample, just adjust "GEN[0]" to reference the sample you're looking for (e.g., "GEN[1]" for the second sample, "GEN[2]" for the third sample, and so on).

ADD COMMENT • link written 11 months ago by Dave Carlson • 320

Please log in to add an answer.

Similar posts • Search »