This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Slicing a vcf file to just a few thousand SNPs with PyVCF.

Hi, I'm trying to reduce a vcf file from 1000 Genomes to just a set of SNPs (around 40k). I have the list of dbSNPs IDs stored on a pandas Series, and I'm trying to retrieve just those genes with a "if record in series:" , but it's not working. Any suggestion on how to call for just certaing SNPs?

Thank you.

python vcf 1000genomes
• 3,592 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by Pedro Morell • Pedro Morell

1 answer

using gatk : https://software.broadinstitute.org/gatk/documentation/tooldocs/current/org_broadinstitute_gatk_tools_walkers_variantutils_SelectVariants.php

--keepIDs / -IDs

List of variant IDs to select
If a file containing a list of IDs is provided to this argument, the tool will only select variants whose ID field is present in this list of IDs. The matching is done by exact string matching. The expected file format is simply plain text with one ID per line
• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

Thanks! I'm trying to use this method, but I'm having an Error with the input. My file looks like this:

rs369720564
rs369711657
rs369708999
rs369700777
rs369687712

So, according to the description (plain text, a single ID per row) it should be working. Any advices?

• 0 views • link
updated by Pierre Lindenbaum Pierre Lindenbaum • written by Pedro Morell • Pedro Morell

ERROR MESSAGE: Invalid argument value 'HCL.txt' at position 4.

Also tryed with a toy list copied from the IDs in the vcf and got the same message, so the format is the issue.

• 0 views • link
written by Pedro Morell • Pedro Morell 
java -jar GenomeAnalysisTK.jar -T SelectVariants --variant input.vcf --keepIDs snps.txt --out output.vcf
• 0 views • link
written by Pedro Morell • Pedro Morell

first, GATK cannot work without the option -R ref.fa

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

second, this is not the cmd line you used as there is no such argument value 'HCL.txt' at position 4.

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

I actually copy/pasted it, but decided to change the names in favor of clarity. I tryed with the reference genome, and now I'm getting this error instead: 

ERROR MESSAGE: Input files /home/hjorvik/Escritorio/TFM/Datos/1KG/ALL.chr6.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf and reference have incompatible contigs. Please see https://software.broadinstitute.org/gatk/documentation/article?id=63for more information. Error details: No overlapping contigs found.

This is the cmd line I used, this time with its original names:

java -jar GenomeAnalysisTK.jar -R GRCh38_full_analysis_set_plus_decoy_hla.fa -T SelectVariants --variant ALL.chr6.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf -IDs HCL.txt --out chrm6_trimed.vcf
• 0 views • link
written by Pedro Morell • Pedro Morell

check #chrom column in your VCF and match with that from fasta file headers (chromosome/contig names)

• 0 views • link
written by cpad0112 cpad0112

I done this, but still have the same problem. I've realized that, even if I changed my chrom column, the contigs in the header are still labeled as 1,2,3,4. How can I change them using gawk?

• 0 views • link
written by Pedro Morell • Pedro Morell

Log in to answer this question.

More posts like this