This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Grep Position List

Hi,

I have VCF file and list of position file. I want to grep those records from VCF which have position corresponding to position present in position file. For this purpose I used this command.

grep -wfF Position.txt myVCF > insec.vcf

here is the list of position.

`POS
10248
10321
10327
10492
10583
12783
13116
13118
13273
13302
14354
14542
14907
14930
15029
15118
15190
15208
15211
15274
`

Here is the result

#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  LP6005441-DNA_G11
chr1    10492   rs55998931      C       T       182.77  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=2.356;ClippingRankSum=0;DB;DP=4
chr1    10583   rs58108140      G       A       170.77  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=-3.006;ClippingRankSum=0;DB;DP=
chr1    12783   rs62635284      G       A       522.77  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=4.292;ClippingRankSum=0;DB;DP=2
chr1    13116   rs62635286      T       G       125.77  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=-0.301;ClippingRankSum=0;DB;DP=
chr1    13118   rs62028691      A       G       125.77  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=-1.009;ClippingRankSum=0;DB;DP=
chr1    13273   rs531730856     G       C       675.77  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=-2.577;ClippingRankSum=0;DB;DP=
chr1    13302   rs180734498     C       T       135.77  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=1.88;ClippingRankSum=0;DB;DP=41
chr1    13417   rs777038595     C       CGAGA   606.73  PASS   AC=1;AF=0.5;AN=2;BaseQRankSum=-2.053;ClippingRankSum=0;DB;DP=
chr1    13896   rs201696125     C       A       195.77  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=-1.19;ClippingRankSum=0;DB;DP=5
chr1    14699   rs62635298      C       G       160.77  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=1.278;ClippingRankSum=0;DB;DP=3
chr1    14907   rs6682375       A       G       922.77  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=-0.114;ClippingRankSum=0;DB;DP=

But it is giving some extra records as well like record containing 13415,13896. These positions are not present in list of positions. What shall I do to overcome this problem?

grep vcf
• 7,414 views • link
updated by Maxime Lamontagne Maxime Lamontagne • written by BAGeno BAGeno

Hi,

Use same command (without -w) by including <tab>YOUR_POSITION<tab> in 'Position.txt' file. This might work.

Example to add tabs: 

sed -i 's/^/\t/' Position.txt

sed -i 's/$/\t/' Position.txt
• 0 views • link
written by EagleEye EagleEye

It worked :) Thank you

• 0 views • link
written by BAGeno BAGeno

3 answers

Grep is not all that well suited for filtering VCF files, or in general, column-oriented files where you'd want to apply the filter on just one column. Perhaps you could extend the pattern to include the leading and trailing tabs.

See bcftools for a whole range of options for filtering that initially may look a little more complicated yet will serve you well as soon as you want to apply more diverse filters:

https://samtools.github.io/bcftools/bcftools.html

• 0 views • link
written by Istvan Albert Istvan Albert

I had also used bcftools isec but it gave empty intersection files. I want to compare two VCF. These files my have same position but different rsIDs on that position.

• 0 views • link
written by BAGeno BAGeno

On the theme of filtering column-oriented files, if you are prepared to code stuff up yourself, consider awk within the shell, or modules like csv within python etc.

• 0 views • link
written by Joe Joe

What about vcftools:

--positions

Include or exclude a set of sites on the basis of a list of positions in a file. Each line of the input file should contain a (tab-separated) chromosome and position. The file can have comment lines that start with a "#", they will be ignored.

• 0 views • link
written by dschika dschika

I tried this command but it did not work.

 vcftools --vcf myNew.vcf --positions Chr.pos.tsv --out Filtered_File
• 0 views • link
written by BAGeno BAGeno

Try adding:

--recode
• 0 views • link
written by dschika dschika

This is giving VCF file with just header.

• 0 views • link
written by BAGeno BAGeno

Hm...could you please double check:

How does Chr.pos.tsv look like? Have you included chromosomes? Are chromosome and position tab-separated?

How does myNew.vcf look like? Is there more than just the header-section? Are there chromosome-position pairs matching your Chr.pos.tsv file?

• 0 views • link
written by dschika dschika

This is what myNew.vcf look like.

##bcftools_filterCommand=filter -i DP>=20 myNew.raw30.vcf
##bcftools_filterCommand=filter -i GQ>=30 New.DP20.vcf
##liftOverProgram=CrossMap(https://sourceforge.net/projects/crossmap/)
##liftOverTime=April01,2017
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  LP6005441-DNA_G11
chr1    10492   rs55998931      C       T       182.77  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=2.356;ClippingRankSum=0;DB;DP=42;ExcessHet=3.0103;FS=0;MLEAC=1;
chr1    10583   rs58108140      G       A       170.77  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=-3.006;ClippingRankSum=0;DB;DP=37;ExcessHet=3.0103;FS=0;MLEAC=1
chr1    12783   rs62635284      G       A       522.77  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=4.292;ClippingRankSum=0;DB;DP=29;ExcessHet=3.0103;FS=0;MLEAC=1;
chr1    13116   rs62635286      T       G       125.77  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=-0.301;ClippingRankSum=0;DB;DP=28;ExcessHet=3.0103;FS=8.671;MLE
chr1    13118   rs62028691      A       G       125.77  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=-1.009;ClippingRankSum=0;DB;DP=27;ExcessHet=3.0103;FS=5.521;MLE
chr1    13273   rs531730856     G       C       675.77  PASS    AC=1;AF=0.5;AN=2;BaseQRankSum=-2.577;ClippingRankSum=0;DB;DP=46;ExcessHet=3.0103;FS=6.947;MLE

This is what Chr.pos.tsv look like.

##reference=file:///groups/reich/reference-genomes/hs37d5/unzipped/hs37d5.fa
##bcftools_viewVersion=1.2+htslib-1.2.1
##bcftools_viewCommand=view -h /n/data1/hms/genetics/reich/1000Genomes/cteam_remap/A-samples/LP6005441-DNA_G11/annotvcf/LP6005441-DNA_G11.annotated.vcf.bgz
##liftOverProgram=CrossMap(https://sourceforge.net/projects/crossmap/)
##liftOverFile=/data1/Sarah/database/hg19ToHg38.over.chain.gz
##new_reference_genome=/data1/Sarah/database/new_hg38.fa
##liftOverTime=May22,2017
#CHROM  POS
1       10248
1       10321
1       10327
1       10492
1       10583
1       12783
1       13116
1       13118
1       13273
• 0 views • link
written by BAGeno BAGeno

Well, "chr1" is not equal "1" ...

• 0 views • link
written by dschika dschika

Why not use R with the merge command?

data1 <- read.table("VCF File", header=F, sep="\t")

data2 <- read.table("Position File", header=F, sep="\t") 
merge.results <- merge(data2, data1, by.x="V1", by.y="V2", all=F)

write.table(merge.results, "VCF Subset", col.names=F, row.names=F, quote=F, sep="\t")
• 0 views • link
updated by WouterDeCoster WouterDeCoster • written by Maxime Lamontagne Maxime Lamontagne

I added code markup to your post for increased readability. You can do this by selecting the text and clicking the 101010 button. When you compose or edit a post that button is in your toolbar, see image below:

101010 Button

• 0 views • link
written by WouterDeCoster WouterDeCoster

Log in to answer this question.

More posts like this
  • PGT only available for some variants in GATK .vcf
    written by michael.flower.14 • michael.flower.14
    I've got a vcf file someone else prepared using GATK. I'm interested in the phasing information in the PGT tag e.g. 0|1. This information seems …
  • Extract the DP from VCF file along with the chromosome postion and alteration
    written by manojkumarbioinfo manojkumarbioinfo
    Hi, I would like to extract the DP from my VCF file, along with the chromosome position and alteration Example VCF file CHROM POS ID …
  • print only columns with data from every line
    written by HL • HL
    Hi, I have a vcf file where is about 60 000 columns. Here is example of the first three lines: #CHROM POS ID REF ALT …
  • Convert missing genotype ./. to homozygous 0/0 in vcf
    written by User000 User000
    Dear all, Is there a way to change missing genotype `./.` to homozygous `0/0` in a vcf file? Fo example: Input: chr1 69849 . G …
  • Strange Allele Frequency in vcf file
    written by GP_PG • GP_PG
    Hi everyone. I performed a snp calling using GATK and following all the suggestions in the BestPractices sections. I performed single snp calling for two …
  • gatk v4.14.0 does not recognize matematical operators
    written by dago dago
    Dear all I have problem with filtering the SNPs called with gatk v4.1.4.0. I have a vcf file that looks like this: contig00001 8227 . …
  • Problem with allele frequency calculation in vcf
    written by BAGeno BAGeno
    Hi, I am encountering problem in vcf related to allele frequency. Here is my sample file. chr1 113834946 rs2476601 A G,T 938986 PASS AC=21,1;AF=0.983315,0.000538213;AN=22;BaseQRankSum=1.13;ClippingRankSum=0.217;DP=33382;ExcessHet=4.0054;FS=0;InbreedingCoeff=-0.0158;MLEAC=1897,1;MLEAF=0.984,0.0005187;MQ=25.78;MQRankSum=0.253;POSITIVE_TRAIN_SITE;QD=28.3;ReadPosRankSum=0.638;SOR=0.773;VQSLOD=5.18;culprit=QD;VQSRMODE=SNP;NS=929;ExcHet=1,0;AA_chimp=G;AA_ensembl=G GT:AD:DP:GQ:PGT:PID:PL …
  • Diploid calls on chrY?
    written by dec986 dec986
    Hello, I recently got my genome sequenced. chrY 59030922 . A G 1244.77 PASS AC=1;AF=0.500;AN=2;BaseQRankSum=-6.580e-01;ClippingRankSum=0.00;DP=41;ExcessHet=3.0103;FS=9.230;MLEAC=1;MLEAF=0.500;MQ=37.33;MQRankSum=-4.563e+00;QD=30.36;ReadPosRankSum=-7.570e-01;SOR=0.076;VQSLOD=0.452;culprit=FS GT:AD:DP:GQ:PL 0/1:9,32:41:99:1273,0,180 chrY 59032665 . A G 418.77 PASS AC=1;AF=0.500;AN=2;BaseQRankSum=-4.182e+00;ClippingRankSum=0.00;DP=52;ExcessHet=3.0103;FS=0.000;MLEAC=1;MLEAF=0.500;MQ=45.59;MQRankSum=-6.082e+00;QD=8.21;ReadPosRankSum=-1.528e+00;SOR=0.756;VQSLOD=2.78;culprit=DP …
  • Is this GT a mistake?
    written by aostern • aostern
    I started with a VCF that had this variant: #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT [SAMPLE] 12 7045891 . ACAGCAGCAGCAGCAGCAGCAG ACAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAG,ACAGCAG,ACAGCAGCAGCAGCAGCAG,ACAGCAGCAGCAGCAGCAGCAGCAG,A,ACAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAGCAG 7038.49 …
  • Mapping errors in SEGDUP regions
    written by arne • arne
    Hi, Are there any best practices for variant calling how to avoid mapping errors in segmental duplication regions? I have a computer science background and …