Question: calculate Per variant Heterozygosity from VCF file
1
gravatar for SOHAIL
12 months ago by
SOHAIL230
Beijing Institute of Genomics, CAS.
SOHAIL230 wrote:

Hi everybody,

Is there any way to calculate the per variant heterozygosity (i.e. number of 1/0 or 0/1 genotypes observed at given variant site for set of individuals in VCF file) from VCF file?

I knew per individual heterozygosity can be calculated by --het tag from VCFtools

Thanks!

-sohail

ngs vcf • 1.9k views
ADD COMMENTlink modified 12 months ago by Kevin Blighe35k • written 12 months ago by SOHAIL230
6
gravatar for Kevin Blighe
12 months ago by
Kevin Blighe35k
Republic of Ireland
Kevin Blighe35k wrote:

Assuming that you're confident that your VCF is normalised and that you genuinely just want to count occurrences of heterozygous calls per line, then the following will work for either phased (0|1 or 1|0) or unphased (0/1 or 1/0) genotypes, or a mixture of these.

Here, I'm actually accessing a BCF and it's chr1 from 1000 Genomes Phase III data. So, the genotypes are phased. For just un-phased, use gsub(/0\/1|1\/0/,""); for just phased, use gsub(/0\|1|1\|0/,""). the gsub function in AWK conveniently returns the number of matched patterns.

bcftools view chr1.1kg.phase3.v5.bcf | awk -F"\t" 'BEGIN {print "CHR\tPOS\tID\tREF\tALT\tHetCount"} !/^#/ {print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t" gsub(/0\|1|1\|0|0\/1|1\/0/,"")}'
CHR POS     ID          REF                     ALT HetCount
1   10177   .           A                       AC  1490
1   10235   .           T                       TA  6
1   10352   rs145072688 T                       TA  2025
1   10616   rs376342519 CCGCCGTTGCAAAGGCGCGCCG  C   35
1   10642   .           G                       A   21
1   11008   .           C                       G   403
1   11012   .           C                       G   403
1   11063   .           T                       G   15
1   13110   .           G                       A   134
1   13116   rs201725126 T                       G   414
1   13118   rs200579949 A                       G   414
1   13273   .           G                       C   444
1   13284   .           G                       A   7
1   13380   .           C                       G   41
1   13483   .           G                       C   10
1   13494   .           A                       G   7
1   13550   .           G                       A   17
1   14464   .           A                       T   428
1   14599   .           T                       A   711

---------------------------------------------------

Edit 17th September, 2018:

A modified version for counting on bi-allelic records:

bcftools view test.vcf.gz | \
awk -F"\t" '{line=$0} BEGIN {
        print "CHR\tPOS\tID\tREF\tALT\tAltHetCount\tAltHomCount\tRefHomCount"
    } !/^#/ {
        if (gsub(/,/, ",", $5)==0) {
            print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t" gsub(/0\|1|1\|0|0\/1|1\/0/,"") "\t" gsub(/1\/1|1\|1/,"") "\t" gsub(/0\/0|0\|0/,"")
        } else if (gsub(/,/, ",", $5)==1) {
            print $1"\t"$2"\t"$3"\t"$4"\t"$5"\t" gsub(/1\/0|0\/1|1\|0|0\|1|1\|2|2\|1|1\/2|2\/1/,"")","gsub(/2\/0|0\/2|2\|0|0\|2|1\|2|2\|1|1\/2|2\/1/,"",line) "\t" gsub(/1\/1|1\|1/,"")","gsub(/2\/2|2\|2/,"") "\t" gsub(/0\/0|0\|0/,"")
        }
    }'

CHR POS         ID          REF ALT  AltHetCount    AltHomCount RefHomCount
16  56380495    rs78560610  G   A    1              0           10
16  56380502    rs3811361   G   A    3              0           8
16  56380570    rs3811360   A   G    6              0           5
16  56381544    rs11297695  CT  C    6              2           3
16  56381677    .           A   C    1              0           10
16  56381699    rs11640263  C   T    7              2           2
16  56388821    rs148650019 A   G    1              0           10
16  56389119    rs143245214 T   TTA  3              3           5
16  56389827    rs5817056   CA  C    5              1           5
16  56390685    .           CTT C,CT 4,9            0,0         1
16  56390968    rs34097380  GA  G    4              2           5
16  56396430    rs3748401   C   A    1              0           10
16  56396486    rs4924      C   T    5              4           2
16  56397153    rs2617846   G   A    2              5           3
16  56397174    rs2617847   C   G    2              4           3
16  56921684    .           T   TA,A 3,8            1,1         1
ADD COMMENTlink modified 4 months ago • written 12 months ago by Kevin Blighe35k
1

It worked perfectly fine.. thanks for the help.. @Kevin

ADD REPLYlink written 12 months ago by SOHAIL230

You're welcome, Sohail

ADD REPLYlink written 12 months ago by Kevin Blighe35k

This is almost perfect. The problem are sites like this:

0/1:1,1:2:26:0|1:27975_T_C:39,0,26

This site gives me two heterozygous counts instead of 1, I am trying to solve it, adding to gsub something like this but is not working gsub(/^0\|1|^1\|0|^0\/1|^1\/0/,"")

Any Idea?

ADD REPLYlink written 6 weeks ago by shinken12380

Hey, that will require an extra if statement. The code above is designed for multi-allelic records that are encoded as 1/2, 2/2, 1/0, etc. I do not exactly have time to edit this, so, my recommendation would be to just split your multiallelic records via bcftools norm -m-any and then run the above AWK code.

How did you produce your file?

ADD REPLYlink written 6 weeks ago by Kevin Blighe35k

I made the SNP calling with GATK.

ADD REPLYlink written 6 weeks ago by shinken12380

I am not surprised (seriously) - this is an issue that re-occurs every now and then. GATK's programs can produce VCFs in a format that is difficult to use with other scripts/programs.

ADD REPLYlink written 6 weeks ago by Kevin Blighe35k
1

Yep, thats True, I this case, because my file is not phased, the solution was to use your script only with the nonphases counts. like this gsub(/0\/1|1\/0/,"").

Thank you for your script It help me a lot.

ADD REPLYlink written 6 weeks ago by shinken12380
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: 2327 users visited in the last hour