Question: Hi,

I want to calculate statistics of missing data per each site in my vcf file.

Using `vcftools --missing-site` gives wrong stats for several sites.

Is there is any other way to calculate it?

Thank you!

**I have 36 samples and here is an example of the `vcftools --missing-site` output for four positions below**

```
chr1 10616 rs376342519 CCGCCGTTGCAAAGGCGCGCCG C 194.49 PASS AC=3;AF=0.083;AN=36;DB;DP=132;ExcessHet=0.1902;FS=0;InbreedingCoeff=0.4307;MLEAC=5;MLEAF=0.139;MQ=38.12;MQRankSum=-0.842;QD=24.31;SOR=1.329;VQSLOD=-6.944;culprit=MQ GT:AD:DP:GQ 0/0:2,0:2:6 ./.:0,0:0:. ./.:0,0:0:. ./.:0,0:0:. 0/0:5,0:5:15 0/0:2,0:2:3 ./.:0,0:0:. 1/1:0,3:3:10 ./.:0,0:0:. ./.:0,0:0:. 0/0:9,0:9:24 0/0:11,0:11:27 ./.:0,0:0:. ./.:0,0:0:. 0/0:9,0:9:15 0/0:19,0:19:54 0/0:4,0:4:12 0/0:9,0:9:2 ./.:0,0:0:. 0/0:9,0:9:9 0/0:8,0:8:18 0/0:8,0:8:24 0/0:6,0:6:15 ./.:0,0:0:. ./.:0,0:0:. 0/0:2,0:2:3 0/0:6,0:6:6 ./.:0,0:0:. ./.:0,0:0:. ./.:0,0:0:. 0/0:8,0:8:24 ./.:0,0:0:. 0/1:3,2:5:76 ./.:0,0:0:. ./.:0,0:0:. ./.:2,0:2:.
chr2 120680534 rs74654922 A ATGTGTGTG 12887.7 PASS AC=7;AF=0.175;AN=40;DB;DP=751;ExcessHet=3.0103;FS=0;InbreedingCoeff=0.5714;MLEAC=12;MLEAF=0.3;MQ=58.52;NEGATIVE_TRAIN_SITE;QD=29.16;SOR=0.883;VQSLOD=-1.646;culprit=MQ GT:AD:DP:GQ 1|1:0,10:12:30 0|0:0,1:17:48 .|.:0,2:16:. 0|0:0,1:21:60 0/0:0,1:6:16 ./.:38,0:38:. 0|0:0,0:24:72 ./.:21,0:21:. 0/0:0,0:13:39 .|.:0,1:23:. 0/0:0,2:11:48 0|0:0,1:22:63 .|.:0,0:18:. ./.:20,0:20:. ./.:15,0:15:. 1/0:0,12:30:99 0|0:0,0:12:36 1|1:0,19:20:57 .|.:0,1:25:. ./.:22,0:22:. 0|0:0,0:19:57 0|0:0,0:13:39 .|.:0,1:16:. 0|0:0,0:26:78 ./.:26,0:26:. 0|0:0,0:29:87 0|0:0,1:29:84 1|1:0,19:21:57 .|.:0,0:23:. 0|0:0,2:23:63 ./.:14,0:14:. 0|0:0,0:20:60 ./.:11,0:11:. 0|0:0,0:15:45 ./.:32,0:32:. ./.:35,0:35:.
chr22 50808269 . AG A 105.55 PASS AC=2;AF=0.059;AN=34;DP=1217;ExcessHet=0.1296;FS=0;InbreedingCoeff=0.1633;MLEAC=7;MLEAF=0.206;MQ=36.17;NEGATIVE_TRAIN_SITE;QD=30.35;SOR=2.584;VQSLOD=-3.491;culprit=MQ GT:AD:DP:GQ ./.:22,0:22:. 0/0:38,0:38:0 ./.:24,0:24:. ./.:20,0:20:. ./.:26,0:26:. ./.:43,0:43:. ./.:26,0:26:. ./.:19,0:19:. 0/0:30,0:30:0 ./.:26,0:26:. ./.:27,0:27:. 0/0:43,0:43:0 0/0:22,0:22:0 ./.:31,0:31:. ./.:28,0:28:. 0/0:34,0:34:0 ./.:26,0:26:. 1/1:0,2:7:7 0/0:51,0:51:0 0/0:33,0:33:0 ./.:21,0:21:. ./.:9,0:9:. 0/0:22,0:22:0 0/0:28,0:28:0 0/0:41,0:41:0 0/0:61,0:61:0 ./.:54,0:54:. 0/0:43,0:43:0 ./.:36,0:36:. 0/0:52,0:52:0 ./.:17,0:17:. 0/0:51,0:51:0 ./.:19,0:19:. 0/0:63,0:63:0 ./.:34,0:34:. 0/0:43,0:43:0
chr18 67740193 . C CTTTTTTTTTTTTTT 8423.06 PASS AC=9;AF=0.196;AN=46;DP=616;ExcessHet=0;FS=0;InbreedingCoeff=0.618;MLEAC=10;MLEAF=0.217;MQ=57.96;NEGATIVE_TRAIN_SITE;QD=26.16;SOR=0.956;VQSLOD=-2.121;culprit=MQ GT:AD:DP:GQ .|.:0,4:15:. .|.:0,5:15:. 1/1:0,2:5:42 .|.:0,2:14:. 1/1:0,5:6:0 0/0:29,0:29:0 0/0:15,0:15:6 0/0:0,1:6:14 0/0:13,0:13:0 .|.:0,4:12:. 0|0:0,6:26:60 .|.:0,3:11:. .|.:0,1:6:. 0/0:0,1:4:20 1|1:0,20:25:60 0/0:18,0:18:0 .|.:0,4:15:. .|.:0,4:18:. 0/0:29,0:29:12 0/0:0,5:15:90 0/0:0,2:6:18 0/0:0,1:6:38 .|.:0,4:15:. 0|0:0,1:17:46 0/0:0,3:9:44 .|.:0,6:12:. 0|0:0,1:18:51 0/0:21,0:21:6 .|.:0,8:16:. .|.:0,5:17:. .|.:0,4:14:. 1/1:0,1:1:0 0|0:0,2:27:75 1/0:0,1:5:38 0/0:24,0:24:15 0/0:16,0:16:9
```

**And here is the output of `vcftools --missing-site`**

```
CHR POS N_DATA N_GENOTYPE_FILTERED N_MISS F_MISS
chr1 10616 72 0 36 0.5
chr2 120680534 66 0 26 0.393939
chr22 50808269 72 0 38 0.527778
chr18 67740193 59 0 13 0.220339
```

**When in reality it should be this**

```
chr1 10616 72 0 36 0.5
chr2 120680534 72 0 32 0.44
chr22 50808269 72 0 38 0.527778
chr18 67740193 72 0 26 0.361
```

Answer: This python code would work for this example:

import sys

infile = sys.argv[1]
try:
ofile = sys.argv[2]
except IndexError:
ofile = "/dev/stdout"

with open(infile, 'r') as fin:
with open(ofile, 'w') as ofile:
ofile.write("\t".join(["CHR", "POS", "N_DATA", "N_GENOTYPE_FILTERED", "N_MISS", "F_MISS\n"]))

for line in fin:

if line.startswith("#"):
continue

line = line.strip().split()
chrom = line[0]
pos = line[1]
missing = 0

for smp in line[9:]:
missing = missing + smp.split(":")[0].replace("|", "/").split("/").count(".")

total = 2 * len(line[9:])

oline = "\t".join([chrom, pos, str(total), "0", str(missing), str(round(missing / total, 6)), "\n"])
ofile.write(oline)

Assuming it's in a file called. e.g. `missing_sites.py`, it can be invoked by:

python missing_sites.py input.vcf

or

python missing_sites.py input.vcf outfile