From 598fb213514bd6dec55219d66e7030a70d15a16a Mon Sep 17 00:00:00 2001
From: tforest <thomas.forest@college-de-france.fr>
Date: Mon, 18 Oct 2021 16:17:38 +0200
Subject: [PATCH] first commit

---
 vcf_to_sfs.py | 41 +++++++++++++++++++++++++++++++++++++++++
 1 file changed, 41 insertions(+)
 create mode 100755 vcf_to_sfs.py

diff --git a/vcf_to_sfs.py b/vcf_to_sfs.py
new file mode 100755
index 0000000..686c8f6
--- /dev/null
+++ b/vcf_to_sfs.py
@@ -0,0 +1,41 @@
+#!/usr/bin/env python3
+
+"""
+Caution : At the moment for gzipped files only.
+
+
+"""
+
+import gzip
+import sys
+
+# default folded SFS
+folded = True
+
+with gzip.open(sys.argv[1], "rb") as inputgz:
+    line = inputgz.readline()
+    genotypes = []
+    SFS_values = {}
+    while line:
+        line = line.decode('utf-8').strip()
+        if not line.startswith("##") and not line.startswith("#"):
+            FORMAT = line.split("\t")[8:9]
+            SAMPLES = line.split("\t")[9:]
+            snp_genotypes = []
+            for sample in SAMPLES:
+                # for UNPHASED data
+                smpl_genotype = [int(a) for a in sample.split(':')[0].split('/') if a != '.']
+                #if not folded:
+                    
+                print(smpl_genotype)
+                nb_alleles = len(set(smpl_genotype))
+                snp_genotypes.append(nb_alleles)
+            print(snp_genotypes)
+            nb_derived_allele = len([val for val in snp_genotypes if val != 0])
+            print("nb derived allele", nb_derived_allele)
+            if nb_derived_allele not in SFS_values.keys():
+                SFS_values[nb_derived_allele] = 1
+            else:
+                SFS_values[nb_derived_allele] += 1
+        line = inputgz.readline()
+        print(SFS_values)