update output plotting of chrom continuity

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

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__":

vcf_to_sfs.py

@@ -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)

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)