modules.plot

  1import gzip
  2import matplotlib.pyplot as plt
  3import numpy as np
  4
  5'''
  6plot.py, handles logic of the 'plot' subcommand
  7'''
  8
  9class AnalyseFile:
 10
 11
 12    def __init__(self, path: str):
 13        self.path = path
 14
 15        # stores the sum of each tail type, indices:
 16        # [0] = no_tail, [1] = polyU, [2] = polyA, 
 17        # [3] = mixed_AU, [4] = mixed_GC, [5] = other
 18        self.type_count_arr = np.zeros(6, dtype=np.uint32)
 19
 20        # stores the number of rows in input
 21        self.row_count = 0
 22
 23    # === File handling methods ===
 24
 25    def _get_header(self):
 26        '''
 27        Internal method, used for reading the header of input text file.
 28
 29        Args:
 30            handle: file handle generated by one of the open(file, 'r') variants, eg. gzip.open(file, 'rt')
 31
 32        Returns:
 33            list[str, ...]: header as a \t separated list of strings
 34        '''
 35        return self.handle.readline().strip().split('\t')
 36
 37    def open_gzip(self):
 38        '''
 39        Opens a gzipped file for reading as text and stores its 1st line (header) separately
 40        '''
 41        self.handle = gzip.open(self.path, "rt")
 42        self.header = self._get_header()
 43
 44    def open_plain(self):
 45        '''
 46        Opens a plain text file for reading and stores its 1st line (header) separately
 47        '''
 48        self.handle = open(self.path, 'r')
 49        self.header = self._get_header()
 50
 51    def close(self):
 52        '''
 53        For closing the file after reading
 54        '''
 55        self.handle.close()
 56
 57    def _line_process(self, line):
 58        '''
 59        Used for processing the input lines into manageable data structures (only lists for now)
 60        '''
 61        return line.strip().split('\t')
 62
 63    def __iter__(self):
 64        '''
 65        Return object iterator
 66        '''
 67        return self
 68
 69    def __next__(self):
 70        '''
 71        Define how to get the next element, i.e. line
 72        '''
 73        line = self.handle.readline()
 74        self.line = self._line_process(self.handle.readline())
 75        if line:
 76            self.line = self._line_process(line)
 77            self.row_count += 1     # tick the row counter
 78            self._count_type()
 79        else:
 80            raise StopIteration
 81
 82    # === Statistic calculation methods ===
 83
 84    def _count_type(self) -> None:
 85        '''
 86        Tick the appropriate tail type counter by 1. Value are stored as
 87        np.uint32 in a shape (6,) np array. Indice-value pairs are as follows:
 88        [0] = no_tail, [1] = polyU, [2] = polyA,
 89        [3] = mixed_AU, [4] = mixed_GC, [5] = other
 90        '''
 91        # 13th column stores the tail type
 92        type = self.line[13]
 93
 94        match type:
 95            case "no_tail":
 96                self.type_count_arr[0] += 1
 97            case "polyU":
 98                self.type_count_arr[1] += 1
 99            case "polyA":
100                self.type_count_arr[2] += 1
101            case "mixed_AU":
102                self.type_count_arr[3] += 1
103            case "mixed_GC":
104                self.type_count_arr[4] += 1
105            case "other":
106                self.type_count_arr[5] += 1
107
108    def tail_perc(self) -> np.ndarray:
109        '''
110        Calculate the percentage of reads with each tail type
111
112        Returns:
113            np.ndarray: shape (6,) float64 array with values between 0 and 100
114                representing percentages
115        '''
116        return self.type_count_arr / self.row_count * 100
117
118def plot(filenames: list[str, ...], values: np.ndarray) -> None:
119    '''
120    Make plots from calculated statistics.
121    This function is the main() od plotting. In the future each specific plot
122    should have their own function, and plot() should ONLY be used for
123    coordinating them
124
125    Args:
126        filenames (list[str, ...]): list of filenames to be used as xticks,
127            probably will be changed to "experiment_names" in the next minor
128            refactor
129        values (np.ndarray): 2D array of float64s, rows represent tail types,
130            columns represent filenames/experimental setups
131    
132    Returns:
133        None
134    '''
135
136    # tail types taken into account
137    types = ("no_tail", "polyU", "polyA", "mixed_AU", "mixed_GC", "other")
138
139    # Percentage arrays are stored under their respective tail type keys
140    data = {tail_type: perc_arr for tail_type, perc_arr in zip(types, values)}
141
142    # plot bar width
143
144    width = 0.5
145
146    # === Stacked barplot ===
147    # === log Y scale ===
148
149    bottom = np.zeros(len(filenames))
150
151    for label, value in data.items():
152        plt.bar(filenames, value, width, label=label, bottom=bottom)
153        bottom += value
154
155    plt.title("Tail type percentage")
156    plt.semilogy()
157    plt.xticks(filenames, rotation = 45)
158
159    plt.legend(loc='upper left', bbox_to_anchor=(1.0, 1.0))
160    plt.tight_layout()  # helps avoid clipping
161 
162    plt.savefig("plottest_log.png", dpi=300)
163
164    plt.clf()
165
166    # === linear Y scale ===
167
168    bottom = np.zeros(len(filenames))
169
170    for label, value in data.items():
171        plt.bar(filenames, value, width, label=label, bottom=bottom)
172        bottom += value
173
174    plt.title("Tail type percentage")
175    plt.xticks(filenames, rotation = 45)
176
177    plt.legend(loc='upper left', bbox_to_anchor=(1.0, 1.0))
178    plt.tight_layout()  # helps avoid clipping
179
180    plt.savefig("plottest.png", dpi=300)
181
182    return None
183
184def ploter(infiles: list[str, ...]) -> None:
185    '''
186    Plot figures for input files generated by Slipper analyse. For now only
187    tail type percent per input file is available
188
189    Args:
190        infiles (list[str, ...]): list of filenames or paths to tsv.gz files
191            generated by Slipper analyse
192
193    Returns:
194        None
195    '''
196
197    files = {i: (path, AnalyseFile(path)) for i, path in zip(range(len(infiles)), infiles)}
198    filenames = [filename[0].split('/')[-1] for filename in files.values()]
199
200    for file in files.values():
201
202        # remember file = (path, AnaluseFile(path))
203        # open gziped file for reading
204        file[1].open_gzip()
205
206        # counting rows and tail_type occurences is handled by
207        # AnalyseFile.__next__() method
208        for _ in file[1]:
209            pass
210        
211        file[1].close()
212
213    # AnalyseFile.tail_perc() calculates the relative frequency of each tail
214    # type per input file
215    percs = np.vstack([file[1].tail_perc() for file in files.values()]).transpose()
216
217    # run the main plotting function
218    plot(filenames=filenames, values=percs)