modules.analyse

  1import io, gzip
  2from dataclasses import dataclass, astuple
  3
  4'''
  5analyse.py, handles the analysis of 'extract' output
  6'''
  7
  8'''
  9What this needs to do:
 10    * calculate tail length             DONE
 11    * count nucleotides                 DONE
 12    * assign strand                     DONE
 13    * assign tail type:                 DONE
 14        all G or C -> G/C               DONE
 15        Empty (len = 0) -> no_tail      DONE
 16        all T -> polyU                  DONE
 17        all A -> polyA                  DONE
 18        A or T -> mixed_AU              DONE
 19        otherwise -> other              DONE
 20    * Count terminal Us (Ts)            DONE
 21'''
 22
 23class ExtractFile:
 24    def __init__(self, path: str):
 25        self.path = path
 26
 27    #
 28    # File handling methods
 29    #
 30
 31    def _get_header(self):
 32        '''
 33        Internal method, used for reading the header of input text file.
 34
 35        Args:
 36            handle: file handle generated by one of the open(file, 'r') variants, eg. gzip.open(file, 'rt')
 37
 38        Returns:
 39            list[str, ...]: header as a \t separated list of strings
 40        '''
 41        return self.handle.readline().strip().split('\t')
 42
 43    def open_gzip(self):
 44        '''
 45        Opens a gzipped file for reading as text and stores its 1st line (header) separately
 46        '''
 47        self.handle = gzip.open(self.path, "rt")
 48        self.header = self._get_header()
 49
 50    def open_plain(self):
 51        '''
 52        Opens a plain text file for reading and stores its 1st line (header) separately
 53        '''
 54        self.handle = open(self.path, 'r')
 55        self.header = self._get_header()
 56
 57    def close(self):
 58        '''
 59        For closing the file after reading
 60        '''
 61        self.handle.close()
 62
 63    def _line_process(self, line):
 64        '''
 65        Used for processing the input lines into manageable data structures (only lists for now)
 66        '''
 67        return line.strip().split('\t')
 68
 69    def __iter__(self):
 70        '''
 71        Return object iterator
 72        '''
 73        return self
 74
 75    def __next__(self):
 76        '''
 77        Define how to get the next element, i.e. line
 78        '''
 79        line = self.handle.readline()
 80        self.line = self._line_process(self.handle.readline())
 81        if line:
 82            self.line = self._line_process(line)
 83        else:
 84            raise StopIteration
 85
 86    #
 87    # File processing/analysis methods
 88    #
 89
 90    def tail_len(self):
 91        '''
 92        Calculate tail length (CLIP3 field)
 93        '''
 94        if self.line[5] != "NA":
 95            return len(self.line[5])
 96        else:
 97            return 0
 98
 99    def count_bases(self):
100        '''
101        Count bases for 3' tail sequence
102        '''
103        tail_seq = self.line[5]
104
105        if tail_seq != "NA":
106            return tuple([tail_seq.count('A'),
107                          tail_seq.count('C'),
108                          tail_seq.count('G'),
109                          tail_seq.count('T')])
110        else:
111            return tuple([0, 0, 0, 0])
112
113    def get_strand(self):
114        '''
115        Get mapping strand from SAM flags. Flags are stored as `int:12-bit`
116        format, eg. 123:000001111011.
117
118        Outputs the Entrez strand convention, i.e. + = 1, - = -1
119        '''
120        flags = bin(int(self.line[1].split(':')[0]))
121
122        # 5th bit is "REVERSE - SEQ is reverse complemented"
123        return -1 if flags[-5] else 1
124
125    def tail_type(self) -> str:
126        '''
127        Assign a tail type based on it's base composition
128
129        Args:
130            base_counts (tuple): (A, C, G, T)
131
132        Returns:
133            str: tail type, one of: no_tail, polyU, polyA, mixed_AU, mixed_GC,
134                other
135        '''
136        base_counts = self.count_bases()
137        total = sum(base_counts)
138        if total == 0:
139            return "no_tail"
140
141        elif base_counts[3] == total:
142            return "polyU"
143
144        elif base_counts[0] == total:
145            return "polyA"
146
147        elif base_counts[0] + base_counts[3] == total:
148            return "mixed_AU"
149
150        elif base_counts[2] + base_counts[3] == total:
151            return "mixed_GC"
152
153        else:
154            return "other"
155
156    def count_Us(self):
157        tail = self.line[5]
158        u_count = 0
159
160        for i in range(1, len(tail) + 1):
161            if tail[-i] == 'T':
162                u_count += 1
163            else:
164                break
165
166        return u_count
167
168    def header_format(self):
169        '''
170        Prepare header for printing to file
171        '''
172        outhdr = '\t'.join(self.header + ["TAIL_LEN", 
173                                          "A", "C", "G", "U",
174                                          "STRAND", "TYPE",
175                                          "U_COUNT"]) + '\n'
176
177        return outhdr
178
179    def line_format(self):
180        '''
181        Prepare line for writing to file
182        '''
183        return '\t'.join(self.line + [str(self.tail_len()),
184                                      '\t'.join([str(i) for i in self.count_bases()]),
185                                      str(self.get_strand()),
186                                      self.tail_type(),
187                                      str(self.count_Us())]) + '\n'
188
189def analyser(infile: str, out_name: str, briefmode: bool) -> None:
190
191    f = ExtractFile(infile)
192    f.open_gzip()
193
194    if briefmode:
195        with gzip.open(f"{out_name}.tsv.gz", "wt") as fo:
196
197            fo.write(f.header_format())
198
199            for line in f:
200                if f.tail_type() != "no_tail":
201                    fo.write(f.line_format())
202    else:
203        with gzip.open(f"{out_name}.tsv.gz", "wt") as fo:
204
205            fo.write(f.header_format())
206
207            for line in f:
208                fo.write(f.line_format())
209
210    f.close()
211
212    return None