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

better gestion of memory while parsing vcf

master
tforest 2022-02-15 18:02:06 +01:00
parent a439f0cef6
commit 8d4d36a1f9
2 changed files with 17 additions and 13 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
customgraphics.py
View File

@ -237,8 +237,8 @@ def plot_whole_karyotype(recent_variants, mem_clean = False, step = 1, show = Tr
nb_iter = len(recent_variants) -1 nb_iter = len(recent_variants) -1
if show : if show :
iter_start = min_chr_id + 1 iter_start = min_chr_id + 1
if not step : if step == "auto" :
step = round(len(recent_variants[list(recent_variants.keys())[min_chr_id]]) / step) step = round(len(recent_variants[list(recent_variants.keys())[min_chr_id]]) / 1000)
if stacked: if stacked:
nb_subplots = nb_iter - min_chr_id nb_subplots = nb_iter - min_chr_id
subplot_init = True subplot_init = True
@ -262,6 +262,7 @@ def plot_whole_karyotype(recent_variants, mem_clean = False, step = 1, show = Tr
print("Cleaned mem. in", str(datetime.timedelta(seconds=end - start))) print("Cleaned mem. in", str(datetime.timedelta(seconds=end - start)))
else: else:
# if show is enable, use a step # if show is enable, use a step
if step == "auto":
step = round(len(recent_variants[list(recent_variants.keys())[chr]]) / 1000) step = round(len(recent_variants[list(recent_variants.keys())[chr]]) / 1000)
vcf_utils.customgraphics.plot_chrom_continuity(recent_variants, chr_id = chr, show = False, returned = False, step = step, subplot_id = chr) vcf_utils.customgraphics.plot_chrom_continuity(recent_variants, chr_id = chr, show = False, returned = False, step = step, subplot_id = chr)
# last case # last case

21
vcf_utils.py
View File

@ -43,6 +43,7 @@ def parse_vcf(vcf_file, phased=False, stop_at=None, chr_starts_with="*"):
# # every snp line, not comment or header # # every snp line, not comment or header
if not line.startswith("##") and not line.startswith("#"): if not line.startswith("##") and not line.startswith("#"):
FIELDS = line.split("\t") FIELDS = line.split("\t")
# when line is parsed, delete it to save some memory
CHROM = FIELDS[0] CHROM = FIELDS[0]
POS = int(FIELDS[1]) POS = int(FIELDS[1])
if stop_at: if stop_at:
@ -52,8 +53,8 @@ def parse_vcf(vcf_file, phased=False, stop_at=None, chr_starts_with="*"):
REF = FIELDS[3].split(",") REF = FIELDS[3].split(",")
# ALT is col 5 of VCF # ALT is col 5 of VCF
ALT = FIELDS[4].split(",") ALT = FIELDS[4].split(",")
FORMAT = line.split("\t")[8:9] FORMAT = FIELDS[8:9]
SAMPLES = line.split("\t")[9:] SAMPLES = FIELDS[9:]
QUALITY = float(FIELDS[5]) QUALITY = float(FIELDS[5])
INFO = FIELDS[7] INFO = FIELDS[7]
INFOS = {} INFOS = {}
@ -68,7 +69,7 @@ def parse_vcf(vcf_file, phased=False, stop_at=None, chr_starts_with="*"):
# 1 : missing # 1 : missing
# 2 : deletion among REF # 2 : deletion among REF
# 3 : deletion among ALT # 3 : deletion among ALT
if "./.:." in line \ if "./.:." in SAMPLES \
or len(ALT[0]) > 1 \ or len(ALT[0]) > 1 \
or len(REF[0]) > 1: or len(REF[0]) > 1:
# sites that are not kept # sites that are not kept
@ -104,7 +105,7 @@ def parse_vcf(vcf_file, phased=False, stop_at=None, chr_starts_with="*"):
entries = { entries = {
'POS':POS, 'POS':POS,
'CHR':CHROM, 'CHR':CHROM,
'FIELDS':FIELDS, #'FIELDS':FIELDS,
'REF':REF, 'REF':REF,
'ALT':ALT, 'ALT':ALT,
'FORMAT':FORMAT, 'FORMAT':FORMAT,
@ -173,20 +174,22 @@ def build_polymorph_coverage_matrix(entries, noGenotype, diploid=True, na_omit =
def genotyping_continuity_plot(vcf_entries, def genotyping_continuity_plot(vcf_entries,
verbose=False, verbose=False,
step = 1): step = 1):
last_pos = int(sorted(list(vcf_entries.keys()))[-1]) genotyped_pos = sorted(list(vcf_entries.keys()))
last_pos = genotyped_pos[-1]
x = 0 x = 0
y = 1 y = 1
coords = [[], []] coords = [[], []]
print(last_pos, "sites to scan") print("Chr. len. =", last_pos, "bp \t ; nb. SNPs =", len(genotyped_pos[::step]))
for k, pos in enumerate(range(0, last_pos, step)): for k, pos in enumerate(genotyped_pos[::step]):
if verbose: if verbose:
progress = round(k/int(last_pos))*100 progress = round(k/int(last_pos))*100
if progress % 10 == 0: if progress % 10 == 0:
print(progress, "%") print(progress, "%")
# if pos is genotyped # if pos is genotyped
if k in vcf_entries: # if k in vcf_entries:
# y=k*step
y+=1*step y+=1*step
x+=1*step x=pos
coords[0].append(x) coords[0].append(x)
coords[1].append(y) coords[1].append(y)
return coords return coords