Removing homozygous reference genotypes from multi-sample vcf file

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

5 weeks ago by

seta • 1.2k

Sweden

seta • 1.2k 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 • 206 views

ADD COMMENT • link •

modified 5 weeks ago by finswimmer ♦ 11k • written 5 weeks ago by seta • 1.2k

5 weeks ago by

Pierre Lindenbaum ♦ 121k

France/Nantes/Institut du Thorax - INSERM UMR1087

Pierre Lindenbaum ♦ 121k 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 5 weeks ago • written 5 weeks ago by Pierre Lindenbaum ♦ 121k

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 5 weeks ago by seta • 1.2k

I updated my answer.

ADD REPLY • link written 5 weeks ago by Pierre Lindenbaum ♦ 121k

5 weeks ago by

finswimmer ♦ 11k

Germany

finswimmer ♦ 11k 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 5 weeks ago • written 5 weeks ago by finswimmer ♦ 11k

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

ADD REPLY • link written 5 weeks ago by Pierre Lindenbaum ♦ 121k

Hmm, yes might be.

In this case:

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

ADD REPLY • link written 5 weeks ago by finswimmer ♦ 11k

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 5 weeks ago • written 5 weeks ago by seta • 1.2k

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 5 weeks ago • written 5 weeks ago by finswimmer ♦ 11k

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 5 weeks ago • written 5 weeks ago by genomax ♦ 70k

5 weeks ago by

Dave Carlson • 90

Stony Brook University, NY

Dave Carlson • 90 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 5 weeks ago by Dave Carlson • 90

Please log in to add an answer.

Similar posts • Search »