diff --git a/vcf_to_sfs.py b/vcf_to_sfs.py index fcc49c0..34ca424 100755 --- a/vcf_to_sfs.py +++ b/vcf_to_sfs.py @@ -28,8 +28,6 @@ if diploid and not folded: # initiate SFS_values with a zeros dict SFS_values = dict.fromkeys(range(n),0) -polycount = 0 - with gzip.open(sys.argv[1], "rb") as inputgz: line = inputgz.readline() genotypes = [] @@ -66,12 +64,6 @@ with gzip.open(sys.argv[1], "rb") as inputgz: smpl_genotype = [int(a) for a in sample.split(':')[0].split('|') if a != '.'] nb_alleles = set(smpl_genotype) snp_genotypes += smpl_genotype - # if set(snp_genotypes) > 2: - # polyallelic = set(snp_genotypes) - # else: - # polyallelic = False - polyallelic = len(ALT) - ##print(REF, ALT, snp_genotypes) # skip if all individuals have the same genotype if len(set(snp_genotypes)) == 1: line = inputgz.readline() @@ -80,18 +72,9 @@ with gzip.open(sys.argv[1], "rb") as inputgz: allele_counts[snp_genotypes.count(k)] = k allele_counts_list.append(snp_genotypes.count(k)) if folded : - #allele_counts_list = list(allele_counts.keys()) - ##print("ALC", allele_counts_list, "POLY", polyallelic, ALT) for al in range(polyallelic-1): SFS_values[min(allele_counts_list)-1] += 1/len(ALT) allele_counts_list.remove(min(allele_counts_list)) - # if len(ALT) == 1: - # SFS_values[min(allele_counts_list)-1] += 1 - # else: - # for al in range(polyallelic-1): - # SFS_values[min(allele_counts_list)-1] += 1/len(ALT) - # allele_counts_list.remove(min(allele_counts_list)) - # polycount += 1 line = inputgz.readline() print(SFS_values) print(polycount)