richard
/
geinf
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Cleaned

master
tforest 2021-10-25 11:24:27 +02:00
parent 77991f07ad
commit 6799dd9a53
1 changed files with 0 additions and 17 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

17
vcf_to_sfs.py
View File

@ -28,8 +28,6 @@ if diploid and not folded:
# initiate SFS_values with a zeros dict # initiate SFS_values with a zeros dict
SFS_values = dict.fromkeys(range(n),0) SFS_values = dict.fromkeys(range(n),0)
polycount = 0
with gzip.open(sys.argv[1], "rb") as inputgz: with gzip.open(sys.argv[1], "rb") as inputgz:
line = inputgz.readline() line = inputgz.readline()
genotypes = [] 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 != '.'] smpl_genotype = [int(a) for a in sample.split(':')[0].split('|') if a != '.']
nb_alleles = set(smpl_genotype) nb_alleles = set(smpl_genotype)
snp_genotypes += 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 # skip if all individuals have the same genotype
if len(set(snp_genotypes)) == 1: if len(set(snp_genotypes)) == 1:
line = inputgz.readline() 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[snp_genotypes.count(k)] = k
allele_counts_list.append(snp_genotypes.count(k)) allele_counts_list.append(snp_genotypes.count(k))
if folded : if folded :
#allele_counts_list = list(allele_counts.keys())
##print("ALC", allele_counts_list, "POLY", polyallelic, ALT)
for al in range(polyallelic-1): for al in range(polyallelic-1):
SFS_values[min(allele_counts_list)-1] += 1/len(ALT) SFS_values[min(allele_counts_list)-1] += 1/len(ALT)
allele_counts_list.remove(min(allele_counts_list)) 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() line = inputgz.readline()
print(SFS_values) print(SFS_values)
print(polycount) print(polycount)