Question: Extract samples with alternative genotype from a multi-sample VCF
0
gravatar for tacrolimus
3 months ago by
tacrolimus80
tacrolimus80 wrote:

Dear all,

I have a multisample vcf and want to extract the IDs of those with rare impactful SNPs with at least one alternative allele at the a loci within the region of interest.

I am running the following on a normalised multi-sample vcf which has been annotated with VEP:

bcftools view -r REGION OF INTEREST | filter_vep -filter 'IMPACT is MODERATE or IMPACT is HIGH' |  filter_vep -filter MAX_AF <0.01 or not MAX_AF'' | bcftools view -i 'GT="alt"' | bcftools query -f "'[%CHROM:%PO %SAMPLE %GT\n]' >output.txt

I work in an airlocked environment with patient data so can't share the VCF I am using but it is a standard multi-sample vcf created by merging using bcftools.

I get an output with a variant and ID per line irrespective of the genotype e.g.:

chr1:123:456 sampleID123 ./.
chr1:123:456 sampleID234 ./.
chr1:123:456 sampleID156 0/1.
chr1:123:456 sampleID123 ./.

etc

I would like to just get the sample IDs with at least one alt genotype.

Many thanks

vep bcftools vcf • 177 views
ADD COMMENTlink modified 3 months ago by abascalfederico1.1k • written 3 months ago by tacrolimus80

chr1:123:456 sampleID156 0/1.

I don't understand your problem. In your example above , there is a "sample IDs with at least one alt genotype."

ADD REPLYlink written 3 months ago by Pierre Lindenbaum134k

yes the 0/1 variant is the one I want but the output also includes those with two ./. variants. I only want those with at least 0/1.

ADD REPLYlink written 12 weeks ago by tacrolimus80
0
gravatar for abascalfederico
3 months ago by
abascalfederico1.1k
Spain
abascalfederico1.1k wrote:

If you just want those lines with a "1" in the third column, you can add an awk command:

bcftools view -r REGION OF INTEREST | filter_vep -filter 'IMPACT is MODERATE or IMPACT is HIGH' |  filter_vep -filter MAX_AF <0.01 or not MAX_AF'' | bcftools view -i 'GT="alt"' | bcftools query -f "'[%CHROM:%PO %SAMPLE %GT\n]' | awk '$3 ~ /1/ { print }' >output.txt
ADD COMMENTlink written 3 months ago by abascalfederico1.1k

True, I could also just grep '0/1' '1/1' etc but I was wondering what the error with the code above was? Thanks for taking the time to answer.

ADD REPLYlink written 12 weeks ago by tacrolimus80
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 2178 users visited in the last hour
_