From a7e78958b29128b12655db4526d322ed7f4af2ff Mon Sep 17 00:00:00 2001 From: tforest Date: Wed, 16 Feb 2022 15:42:34 +0100 Subject: [PATCH] update output plotting of chrom continuity --- customgraphics.py | 26 ++++++----- sfs_tools.py | 4 +- vcf_to_sfs.py | 110 ---------------------------------------------- vcf_utils.py | 3 -- 4 files changed, 18 insertions(+), 125 deletions(-) delete mode 100755 vcf_to_sfs.py diff --git a/customgraphics.py b/customgraphics.py index 727414f..08bcd3f 100644 --- a/customgraphics.py +++ b/customgraphics.py @@ -153,7 +153,11 @@ def plot_matrix(mat, legend=None, color_scale_type="YlGn", cbarlabel = "qt", tit def plot(x, y, outfile = None, outfolder = None, ylab=None, xlab=None, title=None, label = None, show=True, nb_subplots = None, subplot_init = False, - subplot_id = None): + subplot_id = None, output = None, dpi = 300, width = 15, height = 15, plot_init = True): + + # before fig is generated, set its dimensions + if plot_init: + plt.figure(figsize=(width, height)) if subplot_init: # define a certain amount of subplots fig, axs = plt.subplots(nb_subplots) @@ -181,10 +185,10 @@ def plot(x, y, outfile = None, outfolder = None, ylab=None, xlab=None, if title: plt.title(title) if outfile: - plt.savefig(outfile) - else: - if show == True: - plt.show() + plt.savefig(outfile, dpi = dpi) + if show == True: + plt.show() + def scatter(x, y, ylab=None, xlab=None, title=None): plt.scatter(x, y) @@ -208,7 +212,7 @@ def barplot(x, y, ylab=None, xlab=None, title=None): def plot_chrom_continuity(vcf_entries, chr_id, x=None, y=None, outfile = None, outfolder = None, returned=False, show=True, label=True, step=1, nb_subplots = None, - subplot_init = False, subplot_id = None, title = None): + subplot_init = False, subplot_id = None, title = None, plot_init = False): chr_name = list(vcf_entries.keys())[chr_id] if label: label = chr_name @@ -226,10 +230,10 @@ def plot_chrom_continuity(vcf_entries, chr_id, x=None, y=None, outfile = None, xlab = "pos. in ref.", title = title, outfile = outfile, outfolder = outfolder, show=show, label=label, - nb_subplots = nb_subplots, subplot_init = subplot_init, subplot_id = subplot_id) + nb_subplots = nb_subplots, subplot_init = subplot_init, subplot_id = subplot_id, plot_init = plot_init) def plot_whole_karyotype(recent_variants, mem_clean = False, step = 1, show = True, min_chr_id = 0, - max_chr_id = None, stacked = False, title = None): + max_chr_id = None, stacked = False, title = None, outfile = None): coords = [] if max_chr_id : nb_iter = max_chr_id @@ -246,7 +250,7 @@ def plot_whole_karyotype(recent_variants, mem_clean = False, step = 1, show = Tr nb_subplots = None subplot_init = False vcf_utils.customgraphics.plot_chrom_continuity(recent_variants, chr_id = min_chr_id, show = False, returned = False, step = step, - nb_subplots = nb_subplots, subplot_init = subplot_init, subplot_id = min_chr_id) + nb_subplots = nb_subplots, subplot_init = subplot_init, subplot_id = min_chr_id, plot_init = True) else : iter_start = 0 for chr in range(iter_start, nb_iter): @@ -267,7 +271,9 @@ def plot_whole_karyotype(recent_variants, mem_clean = False, step = 1, show = Tr vcf_utils.customgraphics.plot_chrom_continuity(recent_variants, chr_id = chr, show = False, returned = False, step = step, subplot_id = chr) # last case if show == True: - vcf_utils.customgraphics.plot_chrom_continuity(recent_variants, chr_id = nb_iter, show = True, returned = False, step = step, subplot_id = nb_iter, title = title) + vcf_utils.customgraphics.plot_chrom_continuity(recent_variants, chr_id = nb_iter, show = True, returned = False, step = step, subplot_id = nb_iter, + title = title, + outfile = outfile, plot_init = False) # maybe add a clean of recent_variants in extreme cases, before building the plots if show == False: return coords diff --git a/sfs_tools.py b/sfs_tools.py index d06fec4..139f260 100755 --- a/sfs_tools.py +++ b/sfs_tools.py @@ -100,7 +100,7 @@ def sfs_from_vcf(n, vcf_file, folded = True, diploid = True, phased = False, ver if verbose: print("SFS=", SFS_values) print("Pluriallelic sites =", count_pluriall) - return SFS_values + return SFS_values, count_pluriall def barplot_sfs(sfs, folded=True, title = "Barplot"): sfs_val = [] @@ -117,7 +117,7 @@ def barplot_sfs(sfs, folded=True, title = "Barplot"): #build the plot title = title+" [folded="+str(folded)+"]" plt.title(title) - plt.bar(sfs.keys(), sfs_val) + plt.bar([i+1 for i in sfs.keys()], sfs_val) plt.show() if __name__ == "__main__": diff --git a/vcf_to_sfs.py b/vcf_to_sfs.py deleted file mode 100755 index 4b21b39..0000000 --- a/vcf_to_sfs.py +++ /dev/null @@ -1,110 +0,0 @@ -#!/usr/bin/env python3 - -""" -FOREST Thomas (thomas.forest@college-de-france.fr) - -Caution : At the moment for gzipped files only. - -ARGS --------- - -standalone usage : vcf_to_sfs.py VCF.gz nb_indiv - -""" - -import gzip -import sys - -def sfs_from_vcf(n, vcf_file, folded = True, diploid = True, phased = False, verbose = False): - - """ Returns an SFS from a VCF file. - - Parameters - ---------- - n : int - Nb of individuals in sample. - vcf_file : str - SNPs in VCF file format. - - Used to generate a Site Frequency Spectrum (SFS) from a VCF. - - Returns - ------- - dict - Site Frequency Spectrum (SFS) - - - """ - - if diploid and not folded: - n *= 2 - # initiate SFS_values with a zeros dict - SFS_values = dict.fromkeys(range(n),0) - # store nb polyallellic sites - polyall = 0 - with gzip.open(vcf_file, "rb") as inputgz: - line = inputgz.readline() - genotypes = [] - print("Parsing VCF", vcf_file, "... Please wait...") - while line: - # decode gzipped binary lines - line = line.decode('utf-8').strip() - # every snp line, not comment or header - if not line.startswith("##") and not line.startswith("#"): - FIELDS = line.split("\t") - # REF is col 4 of VCF - REF = FIELDS[3].split(",") - # ALT is col 5 of VCF - ALT = FIELDS[4].split(",") - FORMAT = line.split("\t")[8:9] - SAMPLES = line.split("\t")[9:] - snp_genotypes = [] - allele_counts = {} - allele_counts_list = [] - # SKIP the SNP if : - # 1 : missing - # 2 : deletion among REF - # 3 : deletion among ALT - if "./.:." in line \ - or len(ALT[0]) > 1 \ - or len(REF[0]) > 1: - line = inputgz.readline() - continue - for sample in SAMPLES: - if not phased: - # for UNPHASED data - smpl_genotype = [int(a) for a in sample.split(':')[0].split('/') if a != '.'] - else: - # for PHASED - smpl_genotype = [int(a) for a in sample.split(':')[0].split('|') if a != '.'] - nb_alleles = set(smpl_genotype) - snp_genotypes += smpl_genotype - # skip if all individuals have the same genotype - if len(set(snp_genotypes)) == 1: - line = inputgz.readline() - continue - for k in set(snp_genotypes): - allele_counts[snp_genotypes.count(k)] = k - allele_counts_list.append(snp_genotypes.count(k)) - if folded and len(ALT) >= 2: - polyall += 1 - else: - SFS_values[min(allele_counts_list)-1] += 1 - line = inputgz.readline() - if verbose: - print(SFS_values) - return SFS_values, polyall - -if __name__ == "__main__": - - if len(sys.argv) != 3: - print("Need 2 args") - exit(0) - - # PARAM : vcf_file - vcf_file = sys.argv[1] - # PARAM : Nb of indiv - n = int(sys.argv[2]) - - sfs, nb_polyall = sfs_from_vcf(n, vcf_file, folded = True, diploid = True, phased = False) - print(sfs) diff --git a/vcf_utils.py b/vcf_utils.py index 2894a5a..ccf2b46 100755 --- a/vcf_utils.py +++ b/vcf_utils.py @@ -185,9 +185,6 @@ def genotyping_continuity_plot(vcf_entries, progress = round(k/int(last_pos))*100 if progress % 10 == 0: print(progress, "%") - # if pos is genotyped - # if k in vcf_entries: - # y=k*step y+=1*step x=pos coords[0].append(x)