Mercurial > repos > charles-bernard > alfa
comparison ALFA/ALFA.py @ 18:a1e2ab10b317 draft
Uploaded
| author | charles-bernard |
|---|---|
| date | Tue, 11 Oct 2016 09:18:48 -0400 |
| parents | |
| children |
comparison
equal
deleted
inserted
replaced
| 17:e3d439570972 | 18:a1e2ab10b317 |
|---|---|
| 1 #!/usr/bin/python | |
| 2 #-*- coding: utf-8 -*- | |
| 3 | |
| 4 __author__ = 'noel & bahin' | |
| 5 ''' <decription> ''' | |
| 6 | |
| 7 import argparse | |
| 8 import os | |
| 9 import numpy | |
| 10 import sys | |
| 11 import subprocess | |
| 12 import matplotlib.pyplot as plt | |
| 13 import matplotlib.cm as cmx | |
| 14 import matplotlib.colors as colors | |
| 15 import matplotlib.patheffects as PathEffects | |
| 16 import re | |
| 17 from matplotlib.backends.backend_pdf import PdfPages | |
| 18 # To correctly embbed the texts when saving plots in svg format | |
| 19 import matplotlib | |
| 20 matplotlib.rcParams['svg.fonttype'] = 'none' | |
| 21 | |
| 22 ########################################################################## | |
| 23 # FUNCTIONS # | |
| 24 ########################################################################## | |
| 25 | |
| 26 def init_dict(d, key, init): | |
| 27 if key not in d: | |
| 28 d[key] = init | |
| 29 | |
| 30 def get_chromosome_lengths(args): | |
| 31 """ | |
| 32 Parse the file containing the chromosomes lengths. | |
| 33 If no length file is provided, browse the annotation file (GTF) to estimate the chromosome sizes ( | |
| 34 """ | |
| 35 lengths={} | |
| 36 gtf_chrom_names=set() | |
| 37 force_get_lengths = False | |
| 38 # If the user provided the chromosome length file | |
| 39 if args.chr_len: | |
| 40 with open(args.chr_len, 'r') as chr_len_file: | |
| 41 for line in chr_len_file: | |
| 42 lengths[line.split('\t')[0]] = int(line.rstrip().split('\t')[1]) | |
| 43 with open(args.annotation,'r') as gtf_file: | |
| 44 for line in gtf_file: | |
| 45 if not line.startswith('#'): | |
| 46 chrom = line.split('\t')[0] | |
| 47 if chrom not in gtf_chrom_names: | |
| 48 gtf_chrom_names.add(chrom) | |
| 49 for chrom in lengths.keys(): | |
| 50 if chrom not in gtf_chrom_names: | |
| 51 print "Warning: at least one chromosome name ('"+chrom+"') of the file '"+args.chr_len+"'does not match any chromosome name if GTF and was ignored." | |
| 52 #del lengths[chrom] | |
| 53 break | |
| 54 for chrom in gtf_chrom_names: | |
| 55 if force_get_lengths: break | |
| 56 if chrom not in lengths.keys(): | |
| 57 print "WARNING: chromosome name '"+chrom+"' was found in gtf and does not match any chromosome name provided in",args.chr_len+". " | |
| 58 print "\t=> The chromosome lenghts will be approximated using annotations in the GTF file." | |
| 59 continue_value ="" | |
| 60 while continue_value not in {"yes","y","no","n"}: | |
| 61 continue_value = raw_input("\tDo you want to continue ('yes' or 'y')?\n\tElse write 'no' or 'n' to exit the script and check your file of lengths.\n") | |
| 62 if continue_value == "no" or continue_value == "n": | |
| 63 sys.exit("Exiting") | |
| 64 elif continue_value == "yes" or continue_value == "y": | |
| 65 force_get_lengths = True | |
| 66 break | |
| 67 print "Error: use 'yes/y/no/n' only" | |
| 68 if not force_get_lengths: | |
| 69 return lengths | |
| 70 # Otherwise, (or if at least one chromosome was missing in chromosome lengths file) we consider the end of the last annotation of the chromosome in the GTF file as the chromosome length | |
| 71 with open(args.annotation, 'r') as gtf_file: | |
| 72 for line in gtf_file: | |
| 73 if not line.startswith('#'): | |
| 74 chrom = line.split('\t')[0] | |
| 75 end = int(line.split('\t')[4]) | |
| 76 init_dict(lengths, chrom, 0) | |
| 77 lengths[chrom] = max(lengths[chrom], end) | |
| 78 if force_get_lengths: | |
| 79 print "The chromosome lenghts have been approximated using the last annotations in the GTF file." | |
| 80 return lengths | |
| 81 | |
| 82 def write_feature_on_index(feat,chrom, start, stop, sign, stranded_genome_index, unstranded_genome_index=None): | |
| 83 grouped_by_biotype_features = [] | |
| 84 for biotype,categs in feat.iteritems(): | |
| 85 categ_list=[] | |
| 86 for cat in set(categs): | |
| 87 categ_list.append(cat) | |
| 88 grouped_by_biotype_features.append(":".join((str(biotype),",".join(categ_list)))) | |
| 89 stranded_genome_index.write('\t'.join((chrom, start, stop, sign,''))+'\t'.join(grouped_by_biotype_features)+'\n') | |
| 90 if unstranded_genome_index : | |
| 91 unstranded_genome_index.write('\t'.join((chrom, start, stop, '.',''))+'\t'.join(grouped_by_biotype_features)+'\n') | |
| 92 | |
| 93 | |
| 94 def add_info(cpt, feat_values, start, stop, chrom=None, unstranded_genome_index=None, stranded_genome_index = None , biotype_prios=None, coverage=1, categ_prios=None): | |
| 95 """ | |
| 96 From an annotated genomic interval (start/stop positions and associated feature : one or more category(ies)/biotype pair(s) ) | |
| 97 - If a file is provided: write the information at the end of the index file being generated; | |
| 98 - else : browse the features and update the counts of categories/biotypes found in the genome. | |
| 99 """ | |
| 100 ## Writing in the file is it was provided | |
| 101 if stranded_genome_index : | |
| 102 unstranded_features = None | |
| 103 # If only one strand has a feature, this feature will directly be written on the unstranded index | |
| 104 if feat_values[0] == {} : | |
| 105 # A interval with no feature corresponds to a region annotated only on the reverse strand : update 'antisense' counts | |
| 106 stranded_genome_index.write('\t'.join((chrom, start, stop, '+','antisense\n'))) | |
| 107 write_feature_on_index(feat_values[1], chrom ,start, stop, '-', stranded_genome_index, unstranded_genome_index) | |
| 108 else : | |
| 109 if feat_values[1] == {} : | |
| 110 write_feature_on_index(feat_values[0], chrom ,start, stop, '+', stranded_genome_index, unstranded_genome_index) | |
| 111 stranded_genome_index.write('\t'.join((chrom, start, stop, '-','antisense\n'))) | |
| 112 | |
| 113 # Else, the unstranded index should contain the union of plus and minus features | |
| 114 else : | |
| 115 write_feature_on_index(feat_values[0], chrom ,start, stop, '+', stranded_genome_index) | |
| 116 write_feature_on_index(feat_values[1], chrom ,start, stop, '-', stranded_genome_index) | |
| 117 unstranded_feat = dict(feat_values[0], **feat_values[1]) | |
| 118 for name in set(feat_values[0]) & set(feat_values[1]): | |
| 119 unstranded_feat[name]+=feat_values[0][name] | |
| 120 write_feature_on_index(unstranded_feat, chrom ,start, stop, '.', unstranded_genome_index) | |
| 121 | |
| 122 ## Increasing category counter(s) | |
| 123 else : | |
| 124 # Default behavior if no biotype priorities : category with the highest priority for each found biotype has the same weight (1/n_biotypes) | |
| 125 if not biotype_prios: | |
| 126 nb_feat = len(feat_values) | |
| 127 # For every categ(s)/biotype pairs | |
| 128 for feat in feat_values: | |
| 129 cur_prio = 0 | |
| 130 selected_categ = [] | |
| 131 # Separate categorie(s) and biotype | |
| 132 try: | |
| 133 biot,cats = feat.split(":") | |
| 134 # Error if feature equal "antisense" : update the 'antisense/antisense' counts | |
| 135 except ValueError : | |
| 136 try : | |
| 137 cpt[(feat,feat)] += (int(stop) - int(start)) * coverage / float(nb_feat) | |
| 138 except : | |
| 139 cpt[(feat,feat)] = (int(stop) - int(start)) * coverage / float(nb_feat) | |
| 140 return | |
| 141 # Browse the categories and get only the one(s) with highest priority | |
| 142 for cat in cats.split(','): | |
| 143 try: prio = prios[cat] | |
| 144 except: | |
| 145 #TODO Find a way to add unknown categories | |
| 146 if cat not in unknown_feature: | |
| 147 print >> sys.stderr, "Warning: Unknown categorie %s found and ignored\n...\r" %cat, | |
| 148 unknown_feature.append(cat) | |
| 149 continue | |
| 150 if prio > cur_prio : | |
| 151 cur_prio = prio | |
| 152 selected_categ = [cat] | |
| 153 if prio == cur_prio : | |
| 154 if cat != selected_categ : | |
| 155 try: | |
| 156 if cat not in selected_categ : | |
| 157 selected_categ.append(cat) | |
| 158 except TypeError : | |
| 159 selected_categ = [selected_categ,cat] | |
| 160 # Increase each counts by the coverage divided by the number of categories and biotypes | |
| 161 nb_cats = len(selected_categ) | |
| 162 for cat in selected_categ : | |
| 163 try: | |
| 164 cpt[(cat,biot)] += (int(stop) - int(start)) * coverage / (float(nb_feat * nb_cats)) | |
| 165 except KeyError: | |
| 166 cpt[(cat,biot)] = (int(stop) - int(start)) * coverage / (float(nb_feat * nb_cats)) | |
| 167 #else : | |
| 168 #cpt[(cats,biot)] = (int(stop) - int(start)) / float(nb_feat) * coverage | |
| 169 # Else, apply biotype selection according to the priority set | |
| 170 else: | |
| 171 #TODO Add an option to pass biotype priorities and handle it | |
| 172 pass | |
| 173 | |
| 174 def print_chrom(features_dict, chrom, stranded_index_file, unstranded_index_file, cpt_genome): | |
| 175 with open(unstranded_index_file,'a') as findex, open(stranded_index_file,'a') as fstrandedindex: | |
| 176 # Initialization of variables : start position of the first interval and associated features for +/- strands | |
| 177 start = "" | |
| 178 for pos in sorted(features_dict['+'].keys()): | |
| 179 if start != "": | |
| 180 add_info(cpt_genome, [features_plus,features_minus], str(start), str(pos), chrom, stranded_genome_index = fstrandedindex, unstranded_genome_index = findex) | |
| 181 start = pos | |
| 182 features_plus = features_dict['+'][pos] | |
| 183 features_minus = features_dict['-'][pos] | |
| 184 | |
| 185 | |
| 186 def create_genome_index(annotation, unstranded_genome_index, stranded_genome_index,cpt_genome,prios,biotypes,chrom_sizes): | |
| 187 ''' Create an index of the genome annotations and save it in a file''' | |
| 188 print '\n### Generating genome indexes\n', | |
| 189 sys.stdout.flush() | |
| 190 # Initializations | |
| 191 intervals_dict = 0 | |
| 192 max_value = -1 | |
| 193 prev_chrom = '' | |
| 194 reverse_strand = {'+':'-','-':'+'} | |
| 195 i = 0 # Line counter | |
| 196 # Write the chromosomes lengths as comment lines before the genome index | |
| 197 with open(unstranded_genome_index,'w') as findex, open(stranded_genome_index,'w') as fstrandedindex: | |
| 198 for key,value in chrom_sizes.items(): | |
| 199 findex.write("#%s\t%s\n" %(key, value)) | |
| 200 fstrandedindex.write("#%s\t%s\n" %(key, value)) | |
| 201 # Running through the GTF file and writing into genome index files | |
| 202 with open(annotation, 'r') as gtf_file: | |
| 203 for line in gtf_file: | |
| 204 # Print messages after X processed lines | |
| 205 i += 1 | |
| 206 if i % 100000 == 0: | |
| 207 print >> sys.stderr, '\r%s line processed...' %str(i) | |
| 208 print '\r \r. . .', | |
| 209 sys.stdout.flush() | |
| 210 elif i % 20000 == 0: | |
| 211 print '\r \r. . .', | |
| 212 sys.stdout.flush() | |
| 213 elif i % 2000 == 0: | |
| 214 print '.', | |
| 215 sys.stdout.flush() | |
| 216 # Processing lines except comment ones | |
| 217 if not line.startswith('#'): | |
| 218 # Getting the line infos | |
| 219 line_split=line[:-1].split('\t') | |
| 220 chrom = line_split[0] | |
| 221 cat=line_split[2] | |
| 222 start = int(line_split[3]) - 1 | |
| 223 stop = int(line_split[4]) | |
| 224 strand = line_split[6] | |
| 225 antisense = reverse_strand[strand] | |
| 226 biotype=line_split[8].split('biotype')[1].split(';')[0].strip('" ') | |
| 227 feat = [(cat,biotype)] | |
| 228 # Registering chromosome info in the genome index file if this is a new chromosome or a annotation not overlapping previously recorded features | |
| 229 if start > max_value or chrom != prev_chrom: | |
| 230 # Write the previous features | |
| 231 if intervals_dict != 0: | |
| 232 if chrom != prev_chrom : | |
| 233 print_chrom(intervals_dict, prev_chrom, stranded_genome_index, unstranded_genome_index, cpt_genome) | |
| 234 print "\rChromosome '" + prev_chrom + "' registered." | |
| 235 else: | |
| 236 print_chrom(intervals_dict, chrom, stranded_genome_index, unstranded_genome_index, cpt_genome) | |
| 237 prev_chrom = chrom | |
| 238 # (Re)Initializing the chromosome lists | |
| 239 intervals_dict = {strand:{start:{biotype:[cat]}, stop:{}},antisense:{start:{},stop:{}}} | |
| 240 max_value = stop | |
| 241 | |
| 242 # Update the dictionary which represents intervals for every disctinct annotation | |
| 243 else: | |
| 244 # Get intervals on the same strand as the current feature | |
| 245 stranded_intervals = intervals_dict[strand] | |
| 246 start_added = False | |
| 247 for pos in sorted(stranded_intervals.keys()): | |
| 248 #print pos | |
| 249 #print stranded_intervals[pos] | |
| 250 #print | |
| 251 # While the position is below the feature's interval, store the features | |
| 252 if pos < start : | |
| 253 cur_cat_dict = dict(stranded_intervals[pos]) | |
| 254 cur_antisense_dict = dict(intervals_dict[antisense][pos]) | |
| 255 | |
| 256 # If the start position already exists: update it by addind the new feature | |
| 257 elif pos == start : | |
| 258 cur_cat_dict = dict(stranded_intervals[pos]) | |
| 259 cur_antisense_dict = dict(intervals_dict[antisense][pos]) | |
| 260 #print "cur",cur_cat_dict | |
| 261 try : stranded_intervals[pos][biotype] = stranded_intervals[pos][biotype]+[cat] | |
| 262 except KeyError: stranded_intervals[pos][biotype] = [cat] | |
| 263 start_added = True | |
| 264 #print "cur",cur_cat_dict | |
| 265 | |
| 266 elif pos > start : | |
| 267 # Create a new entry for the start position if necessary | |
| 268 if not start_added : | |
| 269 #print "cur",cur_cat_dict | |
| 270 stranded_intervals[start] = dict(cur_cat_dict) | |
| 271 stranded_intervals[start][biotype] = [cat] | |
| 272 #stranded_intervals[pos][biotype].append(cat) | |
| 273 intervals_dict[antisense][start] = cur_antisense_dict | |
| 274 start_added = True | |
| 275 #print "cur",cur_cat_dict | |
| 276 # While the position is below the stop, just add the new feature | |
| 277 if pos < stop : | |
| 278 cur_cat_dict = dict(stranded_intervals[pos]) | |
| 279 cur_antisense_dict = dict(intervals_dict[antisense][pos]) | |
| 280 try: stranded_intervals[pos][biotype] = stranded_intervals[pos][biotype] + [cat] | |
| 281 except KeyError: stranded_intervals[pos][biotype] = [cat] | |
| 282 # Close the created interval : create an entry at the stop position and restore the features | |
| 283 elif pos > stop : | |
| 284 stranded_intervals[stop] = dict(cur_cat_dict) | |
| 285 intervals_dict[antisense][stop] = cur_antisense_dict | |
| 286 break | |
| 287 else : | |
| 288 break | |
| 289 #try: | |
| 290 #cur_cat_dict = list(stranded_intervals[pos][biotype]) | |
| 291 #except KeyError: cur_cat_dict = list() | |
| 292 #print stranded_intervals[pos] | |
| 293 #print | |
| 294 # Extend the dictionary if necessary | |
| 295 if stop > max_value: | |
| 296 max_value = stop | |
| 297 stranded_intervals[stop] = {} | |
| 298 intervals_dict[antisense][stop] = {} | |
| 299 #except KeyError: | |
| 300 #print intervals_dict | |
| 301 #quit() | |
| 302 #intervals_dict[strand] = {strand:{start:{biotype:[cat]}, stop:{}},antisense:{start:{},stop:{}}} | |
| 303 #continue | |
| 304 #for sign in ['-','+']: | |
| 305 #print sign | |
| 306 #for key,val in sorted(intervals_dict[sign].iteritems()): | |
| 307 #print key,'\t',val | |
| 308 #print | |
| 309 #print "-------\n" | |
| 310 | |
| 311 | |
| 312 #Store the categories of the last chromosome | |
| 313 print_chrom(intervals_dict, chrom, stranded_genome_index, unstranded_genome_index, cpt_genome) | |
| 314 print "\rChromosome '" + prev_chrom + "' registered.\nDone!" | |
| 315 | |
| 316 | |
| 317 def create_bedgraph_files(bams,strandness): | |
| 318 samples_files = [] | |
| 319 labels = [] | |
| 320 print "\n### Generating the bedgraph files" | |
| 321 for n in range(0, len(bams), 2): | |
| 322 print "\rProcessing '%s'\n..." %bams[n], | |
| 323 sys.stdout.flush() | |
| 324 #Get the label for this sample | |
| 325 label = bams[n+1] | |
| 326 #Modify it to contain only alphanumeric caracters (avoid files generation with dangerous names) | |
| 327 modified_label = "_".join(re.findall(r"[\w']+", label)) | |
| 328 if strandness in ["reverse","fr-secondstrand"]: | |
| 329 subprocess.call('bedtools genomecov -bg -split -strand - -ibam ' + bams[n] + ' > ' + modified_label + '.plus.bedgraph', shell=True) | |
| 330 subprocess.call('bedtools genomecov -bg -split -strand + -ibam ' + bams[n] + ' > ' + modified_label + '.minus.bedgraph', shell=True) | |
| 331 elif strandness in ["forward","fr-firststrand"]: | |
| 332 subprocess.call('bedtools genomecov -bg -split -strand + -ibam ' + bams[n] + ' > ' + modified_label + '.plus.bedgraph', shell=True) | |
| 333 subprocess.call('bedtools genomecov -bg -split -strand - -ibam ' + bams[n] + ' > ' + modified_label + '.minus.bedgraph', shell=True) | |
| 334 else : | |
| 335 subprocess.call('bedtools genomecov -bg -split -ibam ' + bams[n] + ' > ' + modified_label + '.bedgraph', shell=True) | |
| 336 samples_files.append(modified_label) | |
| 337 labels.append(label) | |
| 338 print "\rDone!" | |
| 339 return samples_files, labels | |
| 340 | |
| 341 def read_gtf(gtf_index_file, sign): | |
| 342 global gtf_line, gtf_chrom, gtf_start, gtf_stop, gtf_features, endGTF | |
| 343 strand = "" | |
| 344 while strand != sign : | |
| 345 gtf_line = gtf_index_file.readline() | |
| 346 # If the GTF file is finished | |
| 347 if not gtf_line: | |
| 348 endGTF = True | |
| 349 return endGTF | |
| 350 splitline = gtf_line.rstrip().split('\t') | |
| 351 try: strand = splitline[3] | |
| 352 # strand information can not be found in the file file header | |
| 353 except IndexError: pass | |
| 354 gtf_chrom = splitline[0] | |
| 355 gtf_features = splitline[4:] | |
| 356 gtf_start, gtf_stop = map(int, splitline[1:3]) | |
| 357 return endGTF | |
| 358 | |
| 359 def read_counts_files(counts_files): | |
| 360 cpt={} | |
| 361 cpt_genome={} | |
| 362 labels=[] | |
| 363 for fcounts in counts_files: | |
| 364 label=os.path.splitext(os.path.basename(fcounts))[0] | |
| 365 labels.append(label) | |
| 366 cpt[label]={} | |
| 367 with open (fcounts,"r") as f: | |
| 368 for line in f: | |
| 369 if line[0]=="#": | |
| 370 continue | |
| 371 line_split=line[:-1].split('\t') | |
| 372 feature=tuple(line_split[0].split(',')) | |
| 373 cpt[label][feature]=float(line_split[1]) | |
| 374 cpt_genome[feature]=float(line_split[2]) | |
| 375 return cpt,cpt_genome,labels | |
| 376 | |
| 377 | |
| 378 def get_chromosome_names_in_index(genome_index): | |
| 379 chrom_list = [] | |
| 380 with open(genome_index, 'r') as findex: | |
| 381 chrom = "" | |
| 382 for line in findex: | |
| 383 cur_chrom = line.split('\t')[0] | |
| 384 if cur_chrom == chrom: | |
| 385 pass | |
| 386 else: | |
| 387 chrom = cur_chrom | |
| 388 if chrom not in chrom_list: | |
| 389 chrom_list.append(chrom) | |
| 390 return set(chrom_list) | |
| 391 | |
| 392 | |
| 393 def intersect_bedgraphs_and_index_to_counts_categories(samples_files,samples_names,prios,genome_index, strandness, biotype_prios = None): | |
| 394 global gtf_line, gtf_chrom, gtf_start, gtf_stop, gtf_cat, endGTF | |
| 395 print "\n### Intersecting files with indexes" | |
| 396 unknown_chrom = [] | |
| 397 cpt = {} # Counter for the nucleotides in the BAM input file(s) | |
| 398 for n in range(len(samples_files)): | |
| 399 sample_file=samples_files[n] | |
| 400 sample_name=samples_names[n] | |
| 401 # Initializing the category counter dict for this sample | |
| 402 init_dict(cpt, sample_name, {}) | |
| 403 if strandness == "unstranded": | |
| 404 strands = [("",".")] | |
| 405 else: | |
| 406 strands = [('.plus','+'), ('.minus','-')] | |
| 407 | |
| 408 # Intersecting the BEDGRAPH and genome index files | |
| 409 print "\rProcessing '%s'\n. . ." %sample_file, | |
| 410 sys.stdout.flush() | |
| 411 | |
| 412 for strand,sign in strands: | |
| 413 prev_chrom = '' | |
| 414 endGTF = False # Reaching the next chr or the end of the GTF index | |
| 415 intergenic_adds = 0.0 | |
| 416 i = 0 | |
| 417 i_chgt = 0 | |
| 418 with open(sample_file + strand + '.bedgraph', 'r') as bam_count_file: | |
| 419 # Running through the BEDGRAPH file | |
| 420 for bam_line in bam_count_file: | |
| 421 i += 1 | |
| 422 if i % 10000 == 0: | |
| 423 print ".", | |
| 424 sys.stdout.flush() | |
| 425 if i % 100000 == 0: | |
| 426 print "\r \r. . .", | |
| 427 sys.stdout.flush() | |
| 428 # Getting the BAM line info | |
| 429 bam_chrom = bam_line.split('\t')[0] | |
| 430 bam_start, bam_stop, bam_cpt = map(float, bam_line.split('\t')[1:4]) | |
| 431 # Skip the line if the chromosome is not in the index | |
| 432 if bam_chrom not in chrom_list: | |
| 433 if bam_chrom not in unknown_chrom: | |
| 434 unknown_chrom.append(bam_chrom) | |
| 435 print "\r \r Chromosome '" + bam_chrom + "' not found in index." | |
| 436 continue | |
| 437 # If this is a new chromosome (or the first one) | |
| 438 if bam_chrom != prev_chrom: | |
| 439 i_chgt = i | |
| 440 intergenic_adds = 0.0 | |
| 441 # (Re)opening the GTF index and looking for the first line of the matching chr | |
| 442 try: gtf_index_file.close() | |
| 443 except UnboundLocalError: pass | |
| 444 gtf_index_file = open(genome_index, 'r') | |
| 445 endGTF = False | |
| 446 read_gtf(gtf_index_file, sign) | |
| 447 while bam_chrom != gtf_chrom: | |
| 448 read_gtf(gtf_index_file, sign) | |
| 449 if endGTF: | |
| 450 break | |
| 451 prev_chrom = bam_chrom | |
| 452 | |
| 453 # Looking for the first matching annotation in the GTF index | |
| 454 while (not endGTF) and (gtf_chrom == bam_chrom) and (bam_start >= gtf_stop): | |
| 455 read_gtf(gtf_index_file, sign) | |
| 456 if gtf_chrom != bam_chrom: | |
| 457 endGTF = True | |
| 458 # Processing BAM lines before the first GTF annotation if there are | |
| 459 if bam_start < gtf_start: | |
| 460 # Increase the 'intergenic' category counter with all or part of the BAM interval | |
| 461 try: | |
| 462 intergenic_adds += min(bam_stop,gtf_start)-bam_start | |
| 463 cpt[sample_name][('intergenic','intergenic')] += (min(bam_stop, gtf_start) - bam_start) * bam_cpt | |
| 464 except KeyError: | |
| 465 cpt[sample_name][('intergenic','intergenic')] = (min(bam_stop, gtf_start) - bam_start) * bam_cpt | |
| 466 # Go to next line if the BAM interval was totally upstream the first GTF annotation, carry on with the remaining part otherwise | |
| 467 if endGTF or (bam_stop <= gtf_start): | |
| 468 continue | |
| 469 else: | |
| 470 bam_start = gtf_start | |
| 471 | |
| 472 # We can start the crossover | |
| 473 while not endGTF: | |
| 474 # Update category counter | |
| 475 add_info(cpt[sample_name], gtf_features, bam_start, min(bam_stop,gtf_stop), coverage = bam_cpt, categ_prios = prios) | |
| 476 # Read the next GTF file line if the BAM line is not entirely covered | |
| 477 if bam_stop > gtf_stop: | |
| 478 # Update the BAM start pointer | |
| 479 bam_start = gtf_stop | |
| 480 endGTF = read_gtf(gtf_index_file, sign) | |
| 481 # If we read a new chromosome in the GTF file then it is considered finished | |
| 482 if bam_chrom != gtf_chrom: | |
| 483 endGTF = True | |
| 484 if endGTF: | |
| 485 break | |
| 486 else: | |
| 487 # Next if the BAM line is entirely covered | |
| 488 bam_start = bam_stop | |
| 489 break | |
| 490 | |
| 491 # Processing the end of the BAM line if necessary | |
| 492 if endGTF and (bam_stop > bam_start): | |
| 493 try: | |
| 494 cpt[sample_name][('intergenic','intergenic')] += (bam_stop - bam_start) * bam_cpt | |
| 495 except KeyError: | |
| 496 cpt[sample_name][('intergenic','intergenic')] = (bam_stop - bam_start) * bam_cpt | |
| 497 gtf_index_file.close() | |
| 498 print "\r \rDone!" | |
| 499 return cpt | |
| 500 | |
| 501 def write_counts_in_files(cpt,genome_counts): | |
| 502 for label,dico in cpt.items(): | |
| 503 label = "_".join(re.findall(r"[\w']+", label)) | |
| 504 with open(label+".categories_counts","w") as fout: | |
| 505 fout.write("#Category,biotype\tCounts_in_bam\tSize_in_genome\n" ) | |
| 506 for feature,counts in dico.items(): | |
| 507 fout.write("%s\t%s\t%s\n" %(','.join(feature),counts,genome_counts[feature])) | |
| 508 | |
| 509 def recategorize_the_counts(cpt,cpt_genome,final): | |
| 510 final_cat_cpt={} | |
| 511 final_genome_cpt={} | |
| 512 for f in cpt: | |
| 513 #print "\nFinal categories for",f,"sample" | |
| 514 final_cat_cpt[f]={} | |
| 515 for final_cat in final: | |
| 516 tot = 0 | |
| 517 tot_genome=0 | |
| 518 for cat in final[final_cat]: | |
| 519 tot += cpt[f][cat] | |
| 520 tot_genome+=cpt_genome[cat] | |
| 521 #output_file.write('\t'.join((final_cat, str(tot))) + '\n') | |
| 522 #print '\t'.join((final_cat, str(tot))) | |
| 523 final_cat_cpt[f][final_cat]=tot | |
| 524 final_genome_cpt[final_cat]=tot_genome | |
| 525 return final_cat_cpt,final_genome_cpt | |
| 526 | |
| 527 def group_counts_by_categ(cpt,cpt_genome,final,selected_biotype): | |
| 528 final_cat_cpt={} | |
| 529 final_genome_cpt={} | |
| 530 filtered_cat_cpt = {} | |
| 531 for f in cpt: | |
| 532 final_cat_cpt[f]={} | |
| 533 filtered_cat_cpt[f] = {} | |
| 534 for final_cat in final: | |
| 535 tot = 0 | |
| 536 tot_filter = 0 | |
| 537 tot_genome=0 | |
| 538 for cat in final[final_cat]: | |
| 539 for key,value in cpt[f].items(): | |
| 540 if key[0] == cat : | |
| 541 tot += value | |
| 542 tot_genome+=cpt_genome[key] | |
| 543 if key[1] == selected_biotype : | |
| 544 tot_filter += value | |
| 545 #output_file.write('\t'.join((final_cat, str(tot))) + '\n') | |
| 546 #print '\t'.join((final_cat, str(tot))) | |
| 547 final_cat_cpt[f][final_cat]=tot | |
| 548 if tot_genome == 0: | |
| 549 final_genome_cpt[final_cat]= 1e-100 | |
| 550 else: | |
| 551 final_genome_cpt[final_cat]=tot_genome | |
| 552 filtered_cat_cpt[f][final_cat]=tot_filter | |
| 553 #if 'antisense' in final_genome_cpt: final_genome_cpt['antisense'] = 0 | |
| 554 return final_cat_cpt,final_genome_cpt,filtered_cat_cpt | |
| 555 | |
| 556 def group_counts_by_biotype(cpt,cpt_genome,biotypes): | |
| 557 final_cpt={} | |
| 558 final_genome_cpt={} | |
| 559 for f in cpt: | |
| 560 final_cpt[f]={} | |
| 561 for biot in biotypes: | |
| 562 tot = 0 | |
| 563 tot_genome=0 | |
| 564 try: | |
| 565 for final_biot in biotypes[biot]: | |
| 566 for key,value in cpt[f].items(): | |
| 567 if key[1] == final_biot : | |
| 568 tot += value | |
| 569 #if key[1] != 'antisense': | |
| 570 tot_genome+=cpt_genome[key] | |
| 571 except: | |
| 572 for key,value in cpt[f].items(): | |
| 573 if key[1] == biot : | |
| 574 tot += value | |
| 575 tot_genome+=cpt_genome[key] | |
| 576 if tot != 0: | |
| 577 final_cpt[f][biot]=tot | |
| 578 final_genome_cpt[biot]=tot_genome | |
| 579 return final_cpt,final_genome_cpt | |
| 580 | |
| 581 #def get_cmap(N): | |
| 582 #'''Returns a function that maps each index in 0, 1, ... N-1 to a distinct | |
| 583 #RGB color.''' | |
| 584 #color_norm = colors.Normalize(vmin=0, vmax=N-1) | |
| 585 #scalar_map = cmx.ScalarMappable(norm=color_norm, cmap='hsv') | |
| 586 #def map_index_to_rgb_color(index): | |
| 587 #return scalar_map.to_rgba(index) | |
| 588 #return map_index_to_rgb_color | |
| 589 | |
| 590 def one_sample_plot(ordered_categs, percentages, enrichment, n_cat, index, index_enrichment, bar_width, counts_type, title) : | |
| 591 ### Initialization | |
| 592 fig = plt.figure(figsize=(13,9)) | |
| 593 ax1 = plt.subplot2grid((2,4),(0,0),colspan=2) | |
| 594 ax2 = plt.subplot2grid((2,4),(1,0),colspan=2) | |
| 595 cmap= plt.get_cmap('Spectral') | |
| 596 cols=[cmap(x) for x in xrange(0,256,256/n_cat)] | |
| 597 if title: | |
| 598 ax1.set_title(title+"in: %s" %samples_names[0]) | |
| 599 else : | |
| 600 ax1.set_title(counts_type+" distribution in mapped reads in: %s" %samples_names[0]) | |
| 601 ax2.set_title('Normalized counts of '+counts_type) | |
| 602 | |
| 603 | |
| 604 ### Barplots | |
| 605 #First barplot: percentage of reads in each categorie | |
| 606 ax1.bar(index, percentages, bar_width, | |
| 607 color=cols) | |
| 608 #Second barplot: enrichment relative to the genome for each categ | |
| 609 # (the reads count in a categ is divided by the categ size in the genome) | |
| 610 ax2.bar(index_enrichment, enrichment, bar_width, | |
| 611 color=cols,) | |
| 612 ### Piecharts | |
| 613 pielabels = [ordered_categs[i] if percentages[i]>0.025 else '' for i in xrange(n_cat)] | |
| 614 sum_enrichment = numpy.sum(enrichment) | |
| 615 pielabels_enrichment = [ordered_categs[i] if enrichment[i]/sum_enrichment>0.025 else '' for i in xrange(n_cat)] | |
| 616 # Categories piechart | |
| 617 ax3 = plt.subplot2grid((2,4),(0,2)) | |
| 618 pie_wedge_collection, texts = ax3.pie(percentages,labels=pielabels, shadow=True, colors=cols) | |
| 619 # Enrichment piechart | |
| 620 ax4 = plt.subplot2grid((2,4),(1,2)) | |
| 621 pie_wedge_collection, texts = ax4.pie(enrichment,labels=pielabels_enrichment, shadow=True, colors=cols) | |
| 622 # Add legends (append percentages to labels) | |
| 623 labels = [" ".join((ordered_categs[i],"({:.1%})".format(percentages[i]))) for i in range(len(ordered_categs))] | |
| 624 ax3.legend(pie_wedge_collection,labels,loc='center',fancybox=True, shadow=True,prop={'size':'medium'}, bbox_to_anchor=(1.7,0.5)) | |
| 625 labels = [" ".join((ordered_categs[i],"({:.1%})".format(enrichment[i]/sum_enrichment))) for i in range(len(ordered_categs))]# if ordered_categs[i] != 'antisense'] | |
| 626 ax4.legend(pie_wedge_collection,labels,loc='center',fancybox=True, shadow=True,prop={'size':'medium'},bbox_to_anchor=(1.7,0.5)) | |
| 627 # Set aspect ratio to be equal so that pie is drawn as a circle | |
| 628 ax3.set_aspect('equal') | |
| 629 ax4.set_aspect('equal') | |
| 630 return fig, ax1, ax2 | |
| 631 | |
| 632 def make_plot(ordered_categs,samples_names,categ_counts,genome_counts,pdf, counts_type, threshold, title = None ,svg = None, png = None): | |
| 633 # From ordered_categs, keep only the features (categs or biotypes) that we can find in at least one sample. | |
| 634 existing_categs = set() | |
| 635 for sample in categ_counts.values(): | |
| 636 existing_categs |= set(sample.keys()) | |
| 637 ordered_categs = filter(existing_categs.__contains__, ordered_categs) | |
| 638 n_cat = len(ordered_categs) | |
| 639 n_exp=len(samples_names) | |
| 640 ##Initialization of the matrix of counts (nrow=nb_experiements, ncol=nb_categorie) | |
| 641 counts=numpy.matrix(numpy.zeros(shape=(n_exp,n_cat))) | |
| 642 for exp in xrange(len(samples_names)): | |
| 643 for cat in xrange(len(ordered_categs)): | |
| 644 try: counts[exp,cat]=categ_counts[samples_names[exp]][ordered_categs[cat]] | |
| 645 except: pass | |
| 646 | |
| 647 ##Normalize the categorie sizes by the total size to get percentages | |
| 648 sizes=[] | |
| 649 sizes_sum=0 | |
| 650 for cat in ordered_categs: | |
| 651 sizes.append(genome_counts[cat]) | |
| 652 sizes_sum+=genome_counts[cat] | |
| 653 if 'antisense' in ordered_categs: | |
| 654 antisense_pos = ordered_categs.index('antisense') | |
| 655 sizes[antisense_pos] = 1e-100 | |
| 656 for cpt in xrange(len(sizes)): | |
| 657 sizes[cpt]/=float(sizes_sum) | |
| 658 | |
| 659 ## Create array which contains the percentage of reads in each categ for every sample | |
| 660 percentages=numpy.array(counts/numpy.sum(counts,axis=1)) | |
| 661 ## Create the enrichment array (counts divided by the categorie sizes in the genome) | |
| 662 enrichment=numpy.array(percentages/sizes) | |
| 663 if 'antisense_pos' in locals(): | |
| 664 for i in xrange(len(samples_names)): | |
| 665 enrichment[i][antisense_pos] = 0 | |
| 666 #enrichment=numpy.log(numpy.array(percentages/sizes)) | |
| 667 for exp in xrange(n_exp): | |
| 668 for i in xrange(n_cat): | |
| 669 val = enrichment[exp][i] | |
| 670 if val > 1: | |
| 671 enrichment[exp][i] = val-1 | |
| 672 elif val == 1 or val == 0: | |
| 673 enrichment[exp][i] = 0 | |
| 674 else: | |
| 675 enrichment[exp][i] = -1/val+1 | |
| 676 | |
| 677 #### Finally, produce the plot | |
| 678 | |
| 679 ##Get the colors from the colormap | |
| 680 ncolor=16 | |
| 681 cmap = ["#e47878", "#68b4e5", "#a3ea9b", "#ea9cf3", "#e5c957", "#a3ecd1", "#e97ca0", "#66d985", "#8e7ae5", "#b3e04b", "#b884e4", "#e4e758", "#738ee3", "#e76688", "#70dddd", "#e49261"] | |
| 682 if n_exp > ncolor: | |
| 683 cmap = plt.get_cmap('Set3',n_exp) | |
| 684 cmap = [cmap(i) for i in xrange(n_exp)] | |
| 685 | |
| 686 ## Parameters for the plot | |
| 687 opacity = 1 | |
| 688 #Create a vector which contains the position of each bar | |
| 689 index = numpy.arange(n_cat) | |
| 690 #Size of the bars (depends on the categs number) | |
| 691 bar_width = 0.9/n_exp | |
| 692 | |
| 693 | |
| 694 ##Initialise the subplot | |
| 695 # if there is only one sample, also plot piecharts | |
| 696 #if n_exp == 1 and counts_type.lower() == 'categories': | |
| 697 #fig, ax1, ax2 = one_sample_plot(ordered_categs, percentages[0], enrichment[0], n_cat, index, bar_width, counts_type, title) | |
| 698 ## If more than one sample | |
| 699 #else: | |
| 700 fig, (ax1,ax2) = plt.subplots(2,figsize=(5+(n_cat+2*n_exp)/3,10)) | |
| 701 # Store the bars objects for enrichment plot | |
| 702 rects = [] | |
| 703 #For each sample/experiment | |
| 704 for i in range(n_exp): | |
| 705 #First barplot: percentage of reads in each categorie | |
| 706 ax1.bar(index+i*bar_width, percentages[i], bar_width, | |
| 707 alpha=opacity, | |
| 708 color=cmap[i], | |
| 709 label=samples_names[i], edgecolor='#FFFFFF', lw=0) | |
| 710 #Second barplot: enrichment relative to the genome for each categ | |
| 711 # (the reads count in a categ is divided by the categ size in the genome) | |
| 712 rects.append( ax2.bar(index+i*bar_width, enrichment[i], bar_width, | |
| 713 alpha=opacity, | |
| 714 color=cmap[i], | |
| 715 label=samples_names[i], edgecolor=cmap[i], lw=0)) | |
| 716 | |
| 717 ## Graphical options for the plot | |
| 718 # Adding of the legend | |
| 719 ax1.legend(loc='best',frameon=False) | |
| 720 #ax2.legend(loc='upper center',bbox_to_anchor=(0.5,-0.1), fancybox=True, shadow=True) | |
| 721 # Main titles | |
| 722 if title: | |
| 723 ax1.set_title(title) | |
| 724 else : | |
| 725 ax1.set_title(counts_type+" distribution in mapped reads") | |
| 726 ax2.set_title('Normalized counts of '+counts_type) | |
| 727 | |
| 728 # Adding enrichment baseline | |
| 729 #ax2.axhline(y=0,color='black',linestyle='dashed',linewidth='1.5') | |
| 730 # Axes limits | |
| 731 ax1.set_xlim(-0.1,len(ordered_categs)+0.1) | |
| 732 if len(sizes) == 1: ax1.set_xlim(-2,3) | |
| 733 ax2.set_xlim(ax1.get_xlim()) | |
| 734 # Set axis limits (max/min values + 5% margin) | |
| 735 ax2_ymin = [] | |
| 736 ax2_ymax = [] | |
| 737 for sample_values in enrichment: | |
| 738 ax2_ymin.append(min(sample_values)) | |
| 739 ax2_ymax.append(max(sample_values)) | |
| 740 ax2_ymax = max(ax2_ymax) | |
| 741 ax2_ymin = min(ax2_ymin) | |
| 742 margin_top, margin_bottom = (abs(0.05*ax2_ymax), abs(0.05*ax2_ymin)) | |
| 743 ax1.set_ylim(0,ax1.get_ylim()[1]*1.05) | |
| 744 if threshold: | |
| 745 threshold_bottom = -abs(float(threshold[0]))+1 | |
| 746 threshold_top = float(threshold[1])-1 | |
| 747 | |
| 748 for i in xrange(n_exp): | |
| 749 for y in xrange(n_cat): | |
| 750 val = enrichment[i][y] | |
| 751 if not numpy.isnan(val) and not (threshold_bottom < val < threshold_top): | |
| 752 rect = rects[i][y] | |
| 753 rect_height = rect.get_height() | |
| 754 if rect.get_y() < 0: | |
| 755 diff = rect_height + threshold_bottom | |
| 756 rect.set_y(threshold_bottom) | |
| 757 ax2_ymin = threshold_bottom | |
| 758 margin_bottom = 0 | |
| 759 else: | |
| 760 diff = rect_height - threshold_top | |
| 761 ax2_ymax = threshold_top | |
| 762 margin_top = 0 | |
| 763 rect.set_height(rect.get_height()-diff) | |
| 764 if margin_top != 0 and margin_bottom != 0: | |
| 765 margin_top, margin_bottom = [max(margin_top, margin_bottom) for i in xrange(2)] | |
| 766 ax2.set_ylim(ax2_ymin-margin_bottom,ax2_ymax+margin_top) | |
| 767 # Y axis title | |
| 768 ax1.set_ylabel('Proportion of reads (%)') | |
| 769 ax2.set_ylabel('Enrichment relative to genome') | |
| 770 # X axis title | |
| 771 ax1.set_xlabel(counts_type) | |
| 772 ax2.set_xlabel(counts_type) | |
| 773 # Add the categories on the x-axis | |
| 774 ax1.set_xticks(index + bar_width*n_exp/2) | |
| 775 ax2.set_xticks(index + bar_width*n_exp/2) | |
| 776 if counts_type.lower() != 'categories': | |
| 777 ax1.set_xticklabels(ordered_categs,rotation='30',ha='right') | |
| 778 ax2.set_xticklabels(ordered_categs,rotation='30',ha='right') | |
| 779 else: | |
| 780 ax1.set_xticklabels(ordered_categs) | |
| 781 ax2.set_xticklabels(ordered_categs) | |
| 782 # Display fractions values in percentages | |
| 783 ax1.set_yticklabels([str(int(i*100)) for i in ax1.get_yticks()]) | |
| 784 # Correct y-axis ticks labels for enrichment subplot | |
| 785 #ax2.set_yticklabels([str(i+1)+"$^{+1}$" if i>0 else 1 if i==0 else str(-(i-1))+"$^{-1}$" for i in ax2.get_yticks()]) | |
| 786 yticks = list(ax2.get_yticks()) | |
| 787 yticks = [ yticks[i]-1 if yticks[i]>9 else yticks[i]+1 if yticks[i]<-9 else yticks[i] for i in xrange(len(yticks))] | |
| 788 ax2.set_yticks(yticks) | |
| 789 ax2.set_yticklabels([str(int(i+1))+"$^{+1}$" if i>0 and i%1==0 else str(i+1)+"$^{+1}$" if i>0 else 1 if i==0 else str(int(-(i-1)))+"$^{-1}$" if i<0 and i%1==0 else str(-(i-1))+"$^{-1}$" for i in ax2.get_yticks()]) | |
| 790 #ax2.set_yticklabels([i+1 if i>0 else 1 if i==0 else "$\\frac{1}{%s}$" %-(i-1) for i in ax2.get_yticks()]) | |
| 791 # Change appearance of 'antisense' bars on enrichment plot since we cannot calculate an enrichment for this artificial category | |
| 792 if 'antisense_pos' in locals(): #ax2.text(antisense_pos+bar_width/2,ax2.get_ylim()[1]/10,'NA') | |
| 793 for i in xrange(n_exp) : | |
| 794 rect = rects[i][antisense_pos] | |
| 795 rect.set_y(ax2.get_ylim()[0]) | |
| 796 rect.set_height(ax2.get_ylim()[1] - ax2.get_ylim()[0]) | |
| 797 rect.set_hatch('/') | |
| 798 rect.set_fill(False) | |
| 799 rect.set_linewidth(0) | |
| 800 #rect.set_color('lightgrey') | |
| 801 #rect.set_edgecolor('#EDEDED') | |
| 802 rect.set_color('#EDEDED') | |
| 803 ax2.text(index[antisense_pos] + bar_width*n_exp/2 - 0.1, (ax2_ymax+ax2_ymin)/2, 'NA') | |
| 804 # Add text for features absent in sample | |
| 805 for i in xrange(n_exp): | |
| 806 for y in xrange(n_cat): | |
| 807 if percentages[i][y] == 0: | |
| 808 txt = ax1.text(y + bar_width*(i+0.5), 0.02, 'Absent in sample', rotation = 'vertical', color = cmap[i], horizontalalignment ='center', verticalalignment = 'bottom') | |
| 809 txt.set_path_effects([PathEffects.Stroke(linewidth=0.5),PathEffects.Normal()]) | |
| 810 elif enrichment[i][y] == 0: | |
| 811 rects[i][y].set_linewidth(1) | |
| 812 | |
| 813 # Remove top/right/bottom axes | |
| 814 for ax in [ax1,ax2]: | |
| 815 ax.spines['top'].set_visible(False) | |
| 816 ax.spines['right'].set_visible(False) | |
| 817 ax.spines['bottom'].set_visible(False) | |
| 818 ax.tick_params(axis='x', which='both', bottom='on', top='off', labelbottom='on') | |
| 819 ax.tick_params(axis='y', which='both', left='on', right='off', labelleft='on') | |
| 820 | |
| 821 | |
| 822 | |
| 823 #ax1.legend(loc='center right', bbox_to_anchor=(1.2, 0),fancybox=True, shadow=True) | |
| 824 ## Showing the plot | |
| 825 plt.tight_layout() | |
| 826 fig.subplots_adjust(wspace=0.2) | |
| 827 if pdf: | |
| 828 pdf.savefig() | |
| 829 plt.close() | |
| 830 elif svg: | |
| 831 if svg == True: | |
| 832 plt.savefig(counts_type+'.svg') | |
| 833 else : | |
| 834 if os.path.splitext(svg)[1] == '.svg': | |
| 835 plt.savefig('.'.join((os.path.splitext(svg)[0],counts_type,'svg'))) | |
| 836 else: | |
| 837 plt.savefig('.'.join((svg,counts_type,'svg'))) | |
| 838 elif png: | |
| 839 if png == True: | |
| 840 plt.savefig(counts_type+'.png') | |
| 841 else : | |
| 842 if os.path.splitext(png)[1] == '.png': | |
| 843 plt.savefig('.'.join((os.path.splitext(png)[0],counts_type,'png'))) | |
| 844 else: | |
| 845 plt.savefig('.'.join((png,counts_type,'png'))) | |
| 846 else: | |
| 847 plt.show() | |
| 848 ## Save on disk it as a png image | |
| 849 #fig.savefig('image_output.png', dpi=300, format='png', bbox_extra_artists=(lgd,), bbox_inches='tight') | |
| 850 | |
| 851 | |
| 852 def filter_categs_on_biotype(selected_biotype,cpt) : | |
| 853 filtered_cpt = {} | |
| 854 for sample in cpt: | |
| 855 filtered_cpt[sample] = {} | |
| 856 for feature,count in cpt[sample].items(): | |
| 857 if feature[1] == selected_biotype: | |
| 858 filtered_cpt[sample][feature[0]] = count | |
| 859 return filtered_cpt | |
| 860 | |
| 861 | |
| 862 ########################################################################## | |
| 863 # MAIN # | |
| 864 ########################################################################## | |
| 865 | |
| 866 def usage_message(name=None): | |
| 867 return ''' | |
| 868 Generate genome indexes: | |
| 869 python ALFA.py -a GTF_FILE [-g GENOME_INDEX] | |
| 870 [--chr_len CHR_LENGTHS_FILE] | |
| 871 Process BAM file(s): | |
| 872 python ALFA.py -g GENOME_INDEX -i BAM1 LABEL1 [BAM2 LABEL2 ...] | |
| 873 [--bedgraph] [-s STRAND] | |
| 874 [-n] [--pdf output.pdf] | |
| 875 [-d {1,2,3,4}] [-t YMIN YMAX] | |
| 876 Index genome + process BAM: | |
| 877 python ALFA.py -a GTF_FILE [-g GENOME_INDEX] | |
| 878 -i BAM1 LABEL1 [BAM2 LABEL2 ...] | |
| 879 [--chr_len CHR_LENGTHS_FILE] | |
| 880 [--bedgraph][-s STRAND] | |
| 881 [-n] [--pdf output.pdf] | |
| 882 [-d {1,2,3,4}] [-t YMIN YMAX] | |
| 883 | |
| 884 Process previously created ALFA counts file(s): | |
| 885 python ALFA.py -c COUNTS1 [COUNTS2 ...] | |
| 886 [-s STRAND] | |
| 887 [-n] [--pdf output.pdf] | |
| 888 [-d {1,2,3,4}] [-t YMIN YMAX] | |
| 889 | |
| 890 ''' | |
| 891 | |
| 892 | |
| 893 #### Parse command line arguments and store them in 'options' | |
| 894 parser = argparse.ArgumentParser(formatter_class=argparse.RawTextHelpFormatter, usage=usage_message()) | |
| 895 parser.add_argument('--version', action='version', version='version 1.0', help="show program's version number and exit\n\n-----------\n\n") | |
| 896 # Options concernant l'index | |
| 897 parser.add_argument('-g','--genome_index', help="Genome index files path and basename for existing index, or path and basename for new index creation\n\n") | |
| 898 parser.add_argument('-a','--annotation', metavar = "GTF_FILE", help='Genomic annotations file (GTF format)\n\n') | |
| 899 parser.add_argument('--chr_len', help='Tabulated file containing chromosome names and lengths\n\n-----------\n\n') | |
| 900 | |
| 901 # Options pour l'étape d'intersection | |
| 902 parser.add_argument('-i','--input','--bam', dest='input', metavar=('BAM_FILE1 LABEL1',""), nargs='+', help='Input BAM file(s) and label(s). The BAM files must be sorted by position.\n\n') | |
| 903 parser.add_argument('--bedgraph', action='store_const',default = False, const = True, help="Use this options if your input file(s) is(are) already in bedgraph format\n\n") | |
| 904 parser.add_argument('-c','--counts',metavar=('COUNTS_FILE',""), nargs='+', help="Use this options instead of '-i/--input' to provide ALFA counts files as input \ninstead of bam/bedgraph files.\n\n") | |
| 905 parser.add_argument('-s','--strandness', dest="strandness", nargs=1, action = 'store', default = ['unstranded'], choices = ['unstranded','forward','reverse','fr-firststrand','fr-secondstrand'], metavar="", help ="Library orientation. Choose within: 'unstranded', 'forward'/'fr-firststrand' \nor 'reverse'/'fr-secondstrand'. (Default: 'unstranded')\n\n-----------\n\n") | |
| 906 | |
| 907 # Options concernant le plot | |
| 908 parser.add_argument('-biotype_filter',nargs=1,help=argparse.SUPPRESS)#"Make an extra plot of categories distribution using only counts of the specified biotype.") | |
| 909 parser.add_argument('-d','--categories_depth', type=int, default='3', choices=range(1,5), help = "Use this option to set the hierarchical level that will be considered in the GTF file (default=3): \n(1) gene,intergenic; \n(2) intron,exon,intergenic; \n(3) 5'UTR,CDS,3'UTR,intron,intergenic; \n(4) start_codon,5'UTR,CDS,3'UTR,stop_codon,intron,intergenic. \n\n") | |
| 910 parser.add_argument('--pdf', nargs='?', default=False, help="Save produced plots in PDF format at specified path ('categories_plots.pdf' if no argument provided)\n\n") | |
| 911 parser.add_argument('--png', nargs='?', default=False, const=True, help="Save produced plots in PNG format with provided argument as basename \nor 'categories.png' and 'biotypes.png' if no argument provided\n\n") | |
| 912 parser.add_argument('--svg', nargs='?', default=False, const=True, help="Save produced plots in SVG format with provided argument as basename \nor 'categories.svg' and 'biotypes.svg' if no argument provided\n\n") | |
| 913 parser.add_argument('-n','--no_plot', dest='quiet', action='store_const', default=False, const=True, help="Do not show plots\n\n") | |
| 914 parser.add_argument('-t','--threshold', dest='threshold', nargs = 2, metavar=("ymin","ymax"), type=float , help="Set axis limits for enrichment plots\n\n") | |
| 915 | |
| 916 if len(sys.argv)==1: | |
| 917 parser.print_usage() | |
| 918 sys.exit(1) | |
| 919 | |
| 920 options = parser.parse_args() | |
| 921 | |
| 922 def required_arg(arg, aliases): | |
| 923 if not arg: | |
| 924 print >> sys.stderr, "\nError: %s argument is missing.\n" %aliases | |
| 925 parser.print_usage() | |
| 926 sys.exit() | |
| 927 | |
| 928 def check_files_enxtension(files): | |
| 929 return | |
| 930 | |
| 931 # Booleans for steps to be executed | |
| 932 make_index = False | |
| 933 intersect_reads = False | |
| 934 process_counts = False | |
| 935 | |
| 936 #### Check arguments conformity and define which steps have to be perfomed | |
| 937 if options.counts : | |
| 938 # Aucun autre argument requis | |
| 939 # Vérifier extension input | |
| 940 | |
| 941 # Action : Faire le plot uniquement | |
| 942 process_counts = True | |
| 943 else: | |
| 944 if options.annotation : | |
| 945 # Vérifier si présence -gi | |
| 946 if options.genome_index : | |
| 947 genome_index_basename = options.genome_index | |
| 948 else: | |
| 949 genome_index_basename = options.annotation.split("/")[-1].split(".gtf")[0] | |
| 950 # Vérifier si un fichier existe déjà: | |
| 951 if os.path.isfile(genome_index_basename+".stranded.index") : | |
| 952 if options.input: | |
| 953 print >> sys.stderr, "\nWarning: a index file named '%s' already exists and will be used. If you want to create a new index, please delete this file or specify an other path." %(genome_index_basename+".stranded.index") | |
| 954 else: | |
| 955 sys.exit("\nError: a index file named %s already exists. If you want to create a new index, please delete this file or specify an other path.\n" %(genome_index_basename+".stranded.index")) | |
| 956 # sinon -> action : index à faire | |
| 957 else : | |
| 958 make_index = True | |
| 959 # si l'index n'est pas à faire : | |
| 960 if options.input: | |
| 961 # Arguments requis: input, genome_index | |
| 962 if 'genome_index_basename' not in locals(): | |
| 963 required_arg(options.genome_index, "-g/--genome_index") | |
| 964 genome_index_basename = options.genome_index | |
| 965 required_arg(options.input, "-i/--input/--bam") | |
| 966 for i in xrange(0, len(options.input), 2): | |
| 967 try : | |
| 968 extension = os.path.splitext(options.input[i+1])[1] | |
| 969 if extension == ".bam" or extension == '.bedgraph' or extension == '.bg': | |
| 970 sys.exit("Error: it seems input files and associated labels are not correctly provided.\n\ | |
| 971 Make sure to follow the expected format : -i Input_file1 Label1 [Input_file2 Label2 ...].") | |
| 972 except: | |
| 973 sys.exit("Error: it seems input files and associated labels are not correctly provided.\nMake sure to follow the expected format : -i Input_file1 Label1 [Input_file2 Label2 ...].") | |
| 974 | |
| 975 intersect_reads = True | |
| 976 # Vérifier input's extension | |
| 977 #TODO | |
| 978 if not (options.counts or options.input or options.annotation): | |
| 979 sys.exit("\nError : some arguments are missing At least '-a', '-c' or '-i' is required. Please refer to help (-h/--help) and usage cases for more details.\n") | |
| 980 if not options.counts: | |
| 981 # Declare genome_index variables | |
| 982 stranded_genome_index = genome_index_basename+".stranded.index" | |
| 983 unstranded_genome_index = genome_index_basename+".unstranded.index" | |
| 984 if options.strandness[0] == "unstranded": | |
| 985 genome_index = unstranded_genome_index | |
| 986 else: | |
| 987 genome_index = stranded_genome_index | |
| 988 | |
| 989 | |
| 990 | |
| 991 #### Initialization of some variables | |
| 992 | |
| 993 # Initializing the category priority order, coding biotypes and the final list | |
| 994 prios = {'start_codon': 7, 'stop_codon': 7, 'five_prime_utr': 6, 'three_prime_utr': 6, 'UTR': 6, 'CDS': 5, 'exon': 4, 'transcript': 3, 'gene': 2, 'antisense': 1,'intergenic': 0} | |
| 995 | |
| 996 biotype_prios = None | |
| 997 #biotype_prios = {"protein_coding":1, "miRNA":2} | |
| 998 | |
| 999 # Possibles groups of categories to plot | |
| 1000 categs_group1={'start': ['start_codon'], '5UTR': ['five_prime_utr','UTR'], 'CDS': ['CDS', 'exon'], '3UTR': ['three_prime_utr'], 'stop': ['stop_codon'], 'introns': ['transcript', 'gene'], 'intergenic': ['intergenic'], 'antisense': ['antisense']} | |
| 1001 categs_group2={'5UTR': ['five_prime_utr', 'UTR'], 'CDS': ['CDS', 'exon','start_codon','stop_codon'], '3UTR': ['three_prime_utr'], 'introns': ['transcript', 'gene'], 'intergenic': ['intergenic'], 'antisense': ['antisense']} | |
| 1002 categs_group3={'exons': ['five_prime_utr', 'three_prime_utr', 'UTR', 'CDS', 'exon','start_codon','stop_codon'], 'introns': ['transcript', 'gene'], 'intergenic': ['intergenic'], 'antisense': ['antisense']} | |
| 1003 categs_group4={'gene': ['five_prime_utr', 'three_prime_utr', 'UTR', 'CDS', 'exon','start_codon','stop_codon', 'transcript', 'gene'], 'intergenic': ['intergenic'], 'antisense': ['antisense']} | |
| 1004 categs_groups = [categs_group4,categs_group3,categs_group2,categs_group1] # Order and merging for the final plot | |
| 1005 cat_list = ['5UTR', 'start', 'CDS', 'stop', '3UTR', 'exons', 'introns', 'gene', 'intergenic', 'antisense'] | |
| 1006 | |
| 1007 | |
| 1008 # biotypes list | |
| 1009 biotypes ={'protein_coding','polymorphic_pseudogene','TR_C_gene','TR_D_gene','TR_J_gene','TR_V_gene','IG_C_gene','IG_D_gene','IG_J_gene','IG_V_gene',"3prime_overlapping_ncrna","lincRNA","macro_lncRNA","miRNA","misc_RNA","Mt_rRNA","Mt_tRNA","processed_transcript","ribozyme","rRNA","scaRNA","sense_intronic","sense_overlapping","snoRNA","snRNA","sRNA","TEC","vaultRNA","antisense","transcribed_processed_pseudogene","transcribed_unitary_pseudogene","transcribed_unprocessed_pseudogene","translated_unprocessed_pseudogene","TR_J_pseudogene","TR_V_pseudogene","unitary_pseudogene","unprocessed_pseudogene","processed_pseudogene","IG_C_pseudogene","IG_J_pseudogene","IG_V_pseudogene","pseudogene","ncRNA","tRNA"} # Type: set (to access quickly) | |
| 1010 | |
| 1011 # Grouping of biotypes: | |
| 1012 biotypes_group1={'protein_coding':['protein_coding'],'pseudogenes':['polymorphic_pseudogene',"transcribed_processed_pseudogene","transcribed_unitary_pseudogene","transcribed_unprocessed_pseudogene","translated_unprocessed_pseudogene","TR_J_pseudogene","TR_V_pseudogene","unitary_pseudogene","unprocessed_pseudogene","processed_pseudogene","IG_C_pseudogene","IG_J_pseudogene","IG_V_pseudogene","pseudogene"],'TR':['TR_C_gene','TR_D_gene','TR_J_gene','TR_V_gene'],'IG':['IG_C_gene','IG_D_gene','IG_J_gene','IG_V_gene'],\ | |
| 1013 'MT_RNA':["Mt_rRNA","Mt_tRNA"],\ | |
| 1014 'ncRNA':["lincRNA","macro_lncRNA","3prime_overlapping_ncrna","ncRNA"],\ | |
| 1015 "others":["misc_RNA","processed_transcript","ribozyme","scaRNA","sense_intronic","sense_overlapping","TEC","vaultRNA"], | |
| 1016 "antisense":["antisense"]} | |
| 1017 for biot in ["miRNA","snoRNA","snRNA","rRNA","sRNA","tRNA"]: | |
| 1018 biotypes_group1[biot]=[biot] | |
| 1019 | |
| 1020 # # Initializing the unkown features lits | |
| 1021 unknown_feature = [] | |
| 1022 | |
| 1023 # Initializing the genome category counter dict | |
| 1024 cpt_genome = {} | |
| 1025 | |
| 1026 | |
| 1027 if process_counts : | |
| 1028 #### If input files are the categories counts, just load them and continue to recategorization step | |
| 1029 cpt,cpt_genome,samples_names = read_counts_files(options.counts) | |
| 1030 else: | |
| 1031 #### Create genome index if needed and get the sizes of categories | |
| 1032 if make_index : | |
| 1033 #### Get the chromosome lengths | |
| 1034 lengths = get_chromosome_lengths(options) | |
| 1035 # Generating the genome index files if the user didn't provide them | |
| 1036 create_genome_index(options.annotation, unstranded_genome_index, stranded_genome_index, cpt_genome, prios, biotypes, lengths) | |
| 1037 | |
| 1038 | |
| 1039 #print '\nChr lengths:', lengths | |
| 1040 | |
| 1041 if intersect_reads: | |
| 1042 # If the indexes already exist, read them to compute the sizes of the categories in the genome and retrieve the chromosome lengths | |
| 1043 if not make_index : | |
| 1044 print "\n### Reading genome indexes\n...\r", | |
| 1045 sys.stdout.flush() | |
| 1046 lengths={} | |
| 1047 with open(genome_index, 'r') as genome_index_file: | |
| 1048 for line in genome_index_file: | |
| 1049 if line[0] == "#": | |
| 1050 lengths[line.split('\t')[0][1:]] = int(line.split('\t')[1]) | |
| 1051 else : | |
| 1052 add_info(cpt_genome, line.rstrip().split('\t')[4:], line.split('\t')[1], line.split('\t')[2], biotype_prios = None, categ_prios = prios) | |
| 1053 | |
| 1054 #### Computing the genome intergenic count: sum of the chr lengths minus sum of the genome annotated intervals | |
| 1055 cpt_genome[('intergenic','intergenic')] = sum(lengths.itervalues()) - sum([v for x,v in cpt_genome.iteritems() if x != ('antisense','antisense')]) | |
| 1056 if not make_index : | |
| 1057 print "Done!" | |
| 1058 #print '\nGenome category counts:' | |
| 1059 #for key,val in cpt_genome.iteritems(): | |
| 1060 #print key,"\t",val | |
| 1061 | |
| 1062 | |
| 1063 #### Create the Bedgraph files if needed and get the files list | |
| 1064 | |
| 1065 if not options.bedgraph: | |
| 1066 # Generating the BEDGRAPH files is the user provided BAM file(s) and get the samples labels (this names will be used in the plot legend) | |
| 1067 samples_files, samples_names = create_bedgraph_files(options.input,options.strandness[0]) | |
| 1068 else: | |
| 1069 # Just initialize the files list with the bedgraph paths | |
| 1070 samples_files = [options.input[i] for i in range(0,len(options.input),2)] | |
| 1071 # and get the labels | |
| 1072 samples_names = [options.input[i] for i in range(1,len(options.input),2)] | |
| 1073 #### Retrieving chromosome names saved in index | |
| 1074 chrom_list = get_chromosome_names_in_index(genome_index) | |
| 1075 #### Processing the BEDGRAPH files: intersecting the bedgraph with the genome index and count the number of aligned positions in each category | |
| 1076 cpt = intersect_bedgraphs_and_index_to_counts_categories(samples_files,samples_names,prios,genome_index, options.strandness[0], biotype_prios = None) | |
| 1077 | |
| 1078 #### Write the counts on disk | |
| 1079 write_counts_in_files(cpt,cpt_genome) | |
| 1080 | |
| 1081 if not (intersect_reads or process_counts) or (options.quiet and options.pdf == False): | |
| 1082 quit("\n### End of program") | |
| 1083 print "\n### Generating plots" | |
| 1084 # Updating the biotypes lists (biotypes and 'biotype_group1'): adding the 'unknow biotypes' found in gtf/index | |
| 1085 if unknown_feature == []: # 'unknown_feature' is define only during the index generation | |
| 1086 # Browse the feature to determine whether some biotypes are 'unknown' | |
| 1087 for sample,counts in cpt.items(): | |
| 1088 for (cat,biot) in counts: | |
| 1089 if biot not in biotypes and cat not in unknown_feature: | |
| 1090 unknown_feature.append(biot) | |
| 1091 for new_biot in unknown_feature: | |
| 1092 biotypes.add(new_biot) | |
| 1093 biotypes_group1["others"].append(new_biot) | |
| 1094 biotypes = sorted(biotypes) | |
| 1095 # move antisense categ to the end of the list | |
| 1096 biotypes.remove('antisense') | |
| 1097 biotypes.append('antisense') | |
| 1098 biotypes_group1 = sorted(biotypes_group1) | |
| 1099 | |
| 1100 | |
| 1101 #print '\nCounts for every category/biotype pair: ',cpt | |
| 1102 | |
| 1103 # Generating plots | |
| 1104 if options.pdf != False: | |
| 1105 if options.pdf == None: | |
| 1106 options.pdf = "categories_plots.pdf" | |
| 1107 pdf = PdfPages(options.pdf) | |
| 1108 else: | |
| 1109 pdf = False | |
| 1110 | |
| 1111 selected_biotype = None | |
| 1112 if options.biotype_filter: | |
| 1113 options.biotype_filter = options.biotype_filter[0] | |
| 1114 for sample in cpt: | |
| 1115 for feature in cpt[sample]: | |
| 1116 biotype = feature[1] | |
| 1117 if options.biotype_filter.lower() == biotype.lower(): | |
| 1118 selected_biotype=biotype | |
| 1119 break | |
| 1120 if selected_biotype == None : | |
| 1121 print "\nError: biotype '"+options.biotype_filter+"' not found. Please check the biotype name and that this biotype exists in your sample(s)." | |
| 1122 sys.exit() | |
| 1123 | |
| 1124 #Print a warning message if the UTRs are not specified as 5' or 3' (they will be ploted as 5'UTR) | |
| 1125 if 'UTR' in [categ[0] for counts in cpt.values() for categ in counts.keys()]: | |
| 1126 print '''\nWARNING: (some) 5'UTR/3'UTR are not precisely defined. Consequently, positions annotated as "UTR" will be counted as "5'UTR"\n''' | |
| 1127 | |
| 1128 #### Make the plot by categories | |
| 1129 #### Recategorizing with the final categories | |
| 1130 final_cats=categs_groups[options.categories_depth-1] | |
| 1131 final_cat_cpt,final_genome_cpt, filtered_cat_cpt = group_counts_by_categ(cpt,cpt_genome,final_cats,selected_biotype) | |
| 1132 #### Display the distribution of specified categories (or biotypes) in samples on a barplot | |
| 1133 # Remove the 'antisense' category if the library type is 'unstranded' | |
| 1134 for dic in cpt.values(): | |
| 1135 if ('antisense','antisense') in dic.keys(): break | |
| 1136 else: | |
| 1137 cat_list.remove('antisense') | |
| 1138 make_plot(cat_list,samples_names,final_cat_cpt,final_genome_cpt,pdf, "categories",options.threshold, svg = options.svg, png = options.png) | |
| 1139 if selected_biotype : | |
| 1140 make_plot(cat_list,samples_names,filtered_cat_cpt,final_genome_cpt,pdf, "categories",options.threshold,title="Categories distribution for '"+selected_biotype+"' biotype", svg = options.svg, png = options.png) | |
| 1141 | |
| 1142 #### Make the plot by biotypes | |
| 1143 #### Recategorizing with the final categories | |
| 1144 final_cat_cpt,final_genome_cpt = group_counts_by_biotype(cpt,cpt_genome,biotypes) | |
| 1145 #### Display the distribution of specified categories (or biotypes) in samples on a barplot | |
| 1146 make_plot(biotypes,samples_names,final_cat_cpt,final_genome_cpt,pdf, "biotypes",options.threshold, svg = options.svg, png = options.png) | |
| 1147 | |
| 1148 | |
| 1149 | |
| 1150 ##### Recategorizing with the final categories | |
| 1151 #final_cat_cpt,final_genome_cpt = group_counts_by_biotype(cpt,cpt_genome,biotypes_group1) | |
| 1152 ##### Display the distribution of specified categories (or biotypes) in samples on a barplot | |
| 1153 #make_plot(biotypes_group1,samples_names,final_cat_cpt,final_genome_cpt,pdf,"Biotype groups", options.threshold, title="Biotypes distribution in mapped reads \n(biotypes are grouped by 'family')", svg = options.svg, png = options.png) | |
| 1154 | |
| 1155 | |
| 1156 if options.pdf: | |
| 1157 pdf.close() | |
| 1158 print "\n### Plots saved in pdf file: %s" %options.pdf | |
| 1159 | |
| 1160 print "\n### End of program" |