Mercurial > repos > devteam > quality_filter
annotate quality_filter.py @ 3:2f4393f9f912 draft
planemo upload commit 00a7926c285bc4a339bd7deebf40b28f39c7d947
| author | devteam |
|---|---|
| date | Tue, 21 Jul 2015 15:40:35 -0400 |
| parents | 6c6f15373f96 |
| children |
| rev | line source |
|---|---|
| 0 | 1 #!/usr/bin/env python |
| 2 #Guruprasad Ananda | |
| 3 """ | |
| 4 Filter based on nucleotide quality (PHRED score). | |
| 5 | |
| 6 usage: %prog input out_file primary_species mask_species score mask_char mask_region mask_region_length | |
| 7 """ | |
| 8 | |
| 9 | |
| 10 from __future__ import division | |
| 11 from galaxy import eggs | |
| 12 import pkg_resources | |
| 13 pkg_resources.require( "lrucache" ) | |
| 14 import numpy | |
| 15 | |
| 16 import sys | |
| 17 import os, os.path | |
| 18 from UserDict import DictMixin | |
| 19 from bx.binned_array import FileBinnedArray | |
| 20 from bx.bitset import * | |
| 21 from bx.bitset_builders import * | |
| 22 from bx.cookbook import doc_optparse | |
|
3
2f4393f9f912
planemo upload commit 00a7926c285bc4a339bd7deebf40b28f39c7d947
devteam
parents:
2
diff
changeset
|
23 from galaxy.tools.exception_handling import * |
| 0 | 24 import bx.align.maf |
| 25 | |
| 26 class FileBinnedArrayDir( DictMixin ): | |
| 27 """ | |
| 28 Adapter that makes a directory of FileBinnedArray files look like | |
| 29 a regular dict of BinnedArray objects. | |
| 30 """ | |
| 31 def __init__( self, dir ): | |
| 32 self.dir = dir | |
| 33 self.cache = dict() | |
| 34 def __getitem__( self, key ): | |
| 35 value = None | |
| 36 if key in self.cache: | |
| 37 value = self.cache[key] | |
| 38 else: | |
| 39 fname = os.path.join( self.dir, "%s.qa.bqv" % key ) | |
| 40 if os.path.exists( fname ): | |
| 41 value = FileBinnedArray( open( fname ) ) | |
| 42 self.cache[key] = value | |
| 43 if value is None: | |
| 44 raise KeyError( "File does not exist: " + fname ) | |
| 45 return value | |
| 46 | |
| 47 def stop_err(msg): | |
| 48 sys.stderr.write(msg) | |
| 49 sys.exit() | |
| 50 | |
| 51 def load_scores_ba_dir( dir ): | |
| 52 """ | |
| 53 Return a dict-like object (keyed by chromosome) that returns | |
| 54 FileBinnedArray objects created from "key.ba" files in `dir` | |
| 55 """ | |
| 56 return FileBinnedArrayDir( dir ) | |
| 57 | |
| 58 def bitwise_and ( string1, string2, maskch ): | |
| 59 result = [] | |
| 60 for i, ch in enumerate(string1): | |
| 61 try: | |
| 62 ch = int(ch) | |
| 63 except: | |
| 64 pass | |
| 65 if string2[i] == '-': | |
| 66 ch = 1 | |
| 67 if ch and string2[i]: | |
| 68 result.append(string2[i]) | |
| 69 else: | |
| 70 result.append(maskch) | |
| 71 return ''.join(result) | |
| 72 | |
| 73 def main(): | |
| 74 # Parsing Command Line here | |
| 75 options, args = doc_optparse.parse( __doc__ ) | |
| 76 | |
| 77 try: | |
| 78 #chr_col_1, start_col_1, end_col_1, strand_col_1 = parse_cols_arg( options.cols ) | |
| 79 inp_file, out_file, pri_species, mask_species, qual_cutoff, mask_chr, mask_region, mask_length, loc_file = args | |
| 80 qual_cutoff = int(qual_cutoff) | |
| 81 mask_chr = int(mask_chr) | |
| 82 mask_region = int(mask_region) | |
| 83 if mask_region != 3: | |
| 84 mask_length = int(mask_length) | |
| 85 else: | |
| 86 mask_length_r = int(mask_length.split(',')[0]) | |
| 87 mask_length_l = int(mask_length.split(',')[1]) | |
| 88 except: | |
| 89 stop_err( "Data issue, click the pencil icon in the history item to correct the metadata attributes of the input dataset." ) | |
| 90 | |
| 91 if pri_species == 'None': | |
| 92 stop_err( "No primary species selected, try again by selecting at least one primary species." ) | |
| 93 if mask_species == 'None': | |
| 94 stop_err( "No mask species selected, try again by selecting at least one species to mask." ) | |
| 95 | |
| 96 mask_chr_count = 0 | |
| 97 mask_chr_dict = {0:'#', 1:'$', 2:'^', 3:'*', 4:'?', 5:'N'} | |
| 98 mask_reg_dict = {0:'Current pos', 1:'Current+Downstream', 2:'Current+Upstream', 3:'Current+Both sides'} | |
| 99 | |
| 100 #ensure dbkey is present in the twobit loc file | |
| 101 try: | |
| 102 pspecies_all = pri_species.split(',') | |
| 103 pspecies_all2 = pri_species.split(',') | |
| 104 pspecies = [] | |
| 105 filepaths = [] | |
| 106 for line in open(loc_file): | |
| 107 if pspecies_all2 == []: | |
| 108 break | |
| 109 if line[0:1] == "#": | |
| 110 continue | |
| 111 fields = line.split('\t') | |
| 112 try: | |
| 113 build = fields[0] | |
| 114 for i, dbkey in enumerate(pspecies_all2): | |
| 115 if dbkey == build: | |
| 116 pspecies.append(build) | |
| 117 filepaths.append(fields[1]) | |
| 118 del pspecies_all2[i] | |
| 119 else: | |
| 120 continue | |
| 121 except: | |
| 122 pass | |
| 123 except Exception, exc: | |
| 124 stop_err( 'Initialization errorL %s' % str( exc ) ) | |
| 125 | |
| 126 if len(pspecies) == 0: | |
| 127 stop_err( "Quality scores are not available for the following genome builds: %s" % ( pspecies_all2 ) ) | |
| 128 if len(pspecies) < len(pspecies_all): | |
| 129 print "Quality scores are not available for the following genome builds: %s" % (pspecies_all2) | |
| 130 | |
| 131 scores_by_chrom = [] | |
| 132 #Get scores for all the primary species | |
| 133 for file in filepaths: | |
| 134 scores_by_chrom.append(load_scores_ba_dir( file.strip() )) | |
| 135 | |
| 136 try: | |
| 137 maf_reader = bx.align.maf.Reader( open(inp_file, 'r') ) | |
| 138 maf_writer = bx.align.maf.Writer( open(out_file,'w') ) | |
| 139 except Exception, e: | |
| 140 stop_err( "Your MAF file appears to be malformed: %s" % str( e ) ) | |
| 141 | |
| 142 maf_count = 0 | |
| 143 for block in maf_reader: | |
| 144 status_strings = [] | |
| 145 for seq in range (len(block.components)): | |
| 146 src = block.components[seq].src | |
| 147 dbkey = src.split('.')[0] | |
| 148 chr = src.split('.')[1] | |
| 149 if not (dbkey in pspecies): | |
| 150 continue | |
| 151 else: #enter if the species is a primary species | |
| 152 index = pspecies.index(dbkey) | |
| 153 sequence = block.components[seq].text | |
| 154 s_start = block.components[seq].start | |
| 155 size = len(sequence) #this includes the gaps too | |
| 156 status_str = '1'*size | |
| 157 status_list = list(status_str) | |
| 158 if status_strings == []: | |
| 159 status_strings.append(status_str) | |
| 160 ind = 0 | |
| 161 s_end = block.components[seq].end | |
| 162 #Get scores for the entire sequence | |
| 163 try: | |
| 164 scores = scores_by_chrom[index][chr][s_start:s_end] | |
| 165 except: | |
| 166 continue | |
| 167 pos = 0 | |
| 168 while pos < (s_end-s_start): | |
| 169 if sequence[ind] == '-': #No score for GAPS | |
| 170 ind += 1 | |
| 171 continue | |
| 172 score = scores[pos] | |
| 173 if score < qual_cutoff: | |
| 174 score = 0 | |
| 175 | |
| 176 if not(score): | |
| 177 if mask_region == 0: #Mask Corresponding position only | |
| 178 status_list[ind] = '0' | |
| 179 ind += 1 | |
| 180 pos += 1 | |
| 181 elif mask_region == 1: #Mask Corresponding position + downstream neighbors | |
| 182 for n in range(mask_length+1): | |
| 183 try: | |
| 184 status_list[ind+n] = '0' | |
| 185 except: | |
| 186 pass | |
| 187 ind = ind + mask_length + 1 | |
| 188 pos = pos + mask_length + 1 | |
| 189 elif mask_region == 2: #Mask Corresponding position + upstream neighbors | |
| 190 for n in range(mask_length+1): | |
| 191 try: | |
| 192 status_list[ind-n] = '0' | |
| 193 except: | |
| 194 pass | |
| 195 ind += 1 | |
| 196 pos += 1 | |
| 197 elif mask_region == 3: #Mask Corresponding position + neighbors on both sides | |
| 198 for n in range(-mask_length_l, mask_length_r+1): | |
| 199 try: | |
| 200 status_list[ind+n] = '0' | |
| 201 except: | |
| 202 pass | |
| 203 ind = ind + mask_length_r + 1 | |
| 204 pos = pos + mask_length_r + 1 | |
| 205 else: | |
| 206 pos += 1 | |
| 207 ind += 1 | |
| 208 | |
| 209 status_strings.append(''.join(status_list)) | |
| 210 | |
| 211 if status_strings == []: #this block has no primary species | |
| 212 continue | |
| 213 output_status_str = status_strings[0] | |
| 214 for stat in status_strings[1:]: | |
| 215 try: | |
| 216 output_status_str = bitwise_and (status_strings[0], stat, '0') | |
| 217 except Exception, e: | |
| 218 break | |
| 219 | |
| 220 for seq in range (len(block.components)): | |
| 221 src = block.components[seq].src | |
| 222 dbkey = src.split('.')[0] | |
| 223 if dbkey not in mask_species.split(','): | |
| 224 continue | |
| 225 sequence = block.components[seq].text | |
| 226 sequence = bitwise_and (output_status_str, sequence, mask_chr_dict[mask_chr]) | |
| 227 block.components[seq].text = sequence | |
| 228 mask_chr_count += output_status_str.count('0') | |
| 229 maf_writer.write(block) | |
| 230 maf_count += 1 | |
| 231 | |
| 232 maf_reader.close() | |
| 233 maf_writer.close() | |
| 234 print "No. of blocks = %d; No. of masked nucleotides = %s; Mask character = %s; Mask region = %s; Cutoff used = %d" % (maf_count, mask_chr_count, mask_chr_dict[mask_chr], mask_reg_dict[mask_region], qual_cutoff) | |
| 235 | |
| 236 | |
| 237 if __name__ == "__main__": | |
| 238 main() |