vcf files - extract SNPs

Question

[paraskevopou]

Hi there all! I have 2 vcf files with filtered SNPs and now I want to extract a. The genes that are shared between the 2 vcf files b. The SNPs that are different on the shared genes between these two vcf files. for example 1.vcf

TRINITY_DN17985_c2_g1    539            0      G     A
TRINITY_DN17985_c2_g1   576         0   G   A 
TRINITY_DN17985_c2_g1   710         0   A   T
TRINITY_DN17985_c2_g1   722         0   T   C
TRINITY_DN17985_c2_g1   822         0   A   G
TRINITY_DN17985_c2_g1   1608    0   T   G
TRINITY_DN17985_c2_g1   2353    0   T   A
TRINITY_DN17985_c2_g1   4371    0   T   A
TRINITY_DN17985_c2_g1   6200    0   C   T
TRINITY_DN17985_c2_g1   7444    0   G   T
TRINITY_DN17985_c2_g1   7456    0   A   T
TRINITY_DN17985_c2_g1   2434    0   G   T
TRINITY_DN17985_c2_g1   4190    0   G   A
TRINITY_DN17985_c2_g1   4275    0   A   G
TRINITY_DN15509_c1_g1   626         0   A   C
TRINITY_DN15509_c1_g1   1647    0   G   A
TRINITY_DN15425_c1_g4   741         0   G   A
TRINITY_DN12988_c0_g1   1553    0   T   A
TRINITY_DN14521_c2_g1   826         0   T   A

2.vcf

TRINITY_DN17985_c2_g1   539         0      G     T
TRINITY_DN17985_c2_g1   576         0   G   T 
TRINITY_DN17985_c2_g1   710         0   A   G
TRINITY_DN17985_c2_g1   722         0   T   C
TRINITY_DN17985_c2_g1   822         0   A   G
TRINITY_DN17985_c2_g1   1608    0   T   G
TRINITY_DN17985_c2_g1   2353    0   T   A
TRINITY_DN17985_c2_g1   4371    0   T   A
TRINITY_DN17985_c2_g1   6200    0   C   T
TRINITY_DN17985_c2_g1   7444    0   G   T
TRINITY_DN17985_c2_g1   7456    0   A   T
TRINITY_DN16743_c0_g1   5739    0   G   A
TRINITY_DN16743_c0_g1   5764    0   T   A
TRINITY_DN16743_c0_g2   199         0   A   T
TRINITY_DN16295_c0_g1   982         0   A   G

and I want to extract only the following for which the ALT allele is different between 1.vcf and 2.vcf files.

TRINITY_DN17985_c2_g1   539         0      G     A
TRINITY_DN17985_c2_g1   576         0   G   A 
TRINITY_DN17985_c2_g1   710         0   A   T

TRINITY_DN17985_c2_g1   539         0      G     T
TRINITY_DN17985_c2_g1   576         0   G   T 
TRINITY_DN17985_c2_g1   710         0   A   G

I know that I can use bcftools isec for getting the SNPs that are shared between the two or the SNPs that are unique in the first and second file respectively. But how I can get first the genes that are present in both files?

Thanks a lot in advance for any help! Sofia

[Matteo Schiavinato]

I suppose this is not the format that the post should have. The VCF format doesn't look like a VCF format. Can you please wrap it around the code flag? (little button with 101010)

[Kevin Blighe]

I have just done that (wrapped it in code format)

Answer 1

[Michael Kosicki]

If the first column is the "genes", then:

bedtools intersect -wa -wb -a 1.vcf -b 2.vcf > shared.vcf
grep -v '^#' shared.vcf | cut -f1 -d' ' shared.vcf | sort | uniq  # list of unique genes

Then, assuming it's the same reference (ie ref is the same on the same position) and that both vcf files actually have 10 columns:

grep -v '^#' shared.vcf | awk '$5!=$15' > diff.vcf   # differing positions, no header

Of course, depending on what you want to do with that file afterwards, you may want to restore the header and remove columns 11-20 (or add them as new rows etc).