diff get_codon_frequency.py @ 17:c87c40e642af draft

Uploaded

--- a/get_codon_frequency.py	Wed May 13 04:34:59 2015 -0400
+++ b/get_codon_frequency.py	Fri May 29 09:17:29 2015 -0400
@@ -23,7 +23,7 @@
 from matplotlib import font_manager
 from matplotlib import colors
 import csv
-from scipy import stats
+from scipy import stats, errstate
 from collections import OrderedDict
 import ribo_functions
 import HTSeq
@@ -60,76 +60,65 @@
                 chrom = feature.iv.chrom
                 start = feature.iv.start
                 stop = feature.iv.end 
-                if start+50 < stop-50 :
+                if start+50 < stop-50 :                
                     region = chrom + ':' + str(start+50) + '-' + str(stop-50)
-                else :
-                    break
-        
-        ## DEPRECATED
-        #for chrom in GFF.iterkeys():
-            #for gene in GFF[chrom] :
-               # codon_dict = init_codon_dict()
-                #start = GFF[chrom][gene]['start']
-                #print start
-                #stop = GFF[chrom][gene]['stop']
-                #print stop
-                #region = chrom + ':' + str(start) + '-' + str(stop)
-        #######
-                # #get all reads in this gene
-                reads = subprocess.check_output(["samtools", "view", bamfile, region])
-                head = subprocess.check_output(["samtools", "view", "-H", bamfile])
-                read_tab = reads.split('\n')        
-                for read in read_tab:
-                    # # search mapper for eliminate multiple alignements
-                    if 'bowtie' in head:
-                        multi_tag = "XS:i:"
-                    elif 'bwa' in  head:
-                        multi_tag = "XT:A:R"
-                    elif 'TopHat' in  head:
-                        tag = "NH:i:1"
-                    else :
-                        stop_err("No PG tag find in "+samfile+". Please use bowtie, bwa or Tophat for mapping")  
-
-                    if len(read) == 0:
-                        continue
-                    len_read = len(read.split('\t')[9])        
-                    # if it's read of good length
-                    if len_read == kmer and (tag in read or multi_tag not in read):
-                        feat = read.split('\t')
-                        seq = feat[9]
-                        # if it's a reverse read
-                        if feat[1] == '16' :
-                            if site == "A" :
-                            # #get A-site
-                                cod = str(Seq(seq[a_pos-5:a_pos-2]).reverse_complement())
-                            elif site == "P" :
-                            # #get P-site
-                                cod = str(Seq(seq[a_pos-2:a_pos+1]).reverse_complement())
-                            else :
-                            # #get site-E
-                                cod = str(Seq(seq[a_pos+1:a_pos+4]).reverse_complement())
-                            # # test if it's a true codon not a CNG codon for example
-                            if codon_dict.has_key(cod) :
-                                codon_dict[cod] += 1
-                        # if it's a forward read
-                        elif feat[1] == '0' :
-                            if site == "A" :
-                            # #get A-site
-                                cod = seq[a_pos:a_pos+3]
-                            elif site == "P" :
-                            # #get P-site
-                                cod = seq[a_pos-3:a_pos]
-                            else :
-                            # #get site-E
-                                cod = seq[a_pos-6:a_pos-3]
-                            if codon_dict.has_key(cod) :
-                                codon_dict[cod] += 1
-                del(read)
+                    # #get all reads in this gene
+                    reads = subprocess.check_output(["samtools", "view", bamfile, region])
+                    head = subprocess.check_output(["samtools", "view", "-H", bamfile])
+                    read_tab = reads.split('\n')        
+                    for read in read_tab:
+                        # # search mapper for eliminate multiple alignements
+                        if 'bowtie' in head:
+                            multi_tag = "XS:i:"
+                        elif 'bwa' in  head:
+                            multi_tag = "XT:A:R"
+                        elif 'TopHat' in  head:
+                            tag = "NH:i:1"
+                        else :
+                            stop_err("No PG tag find in "+samfile+". Please use bowtie, bwa or Tophat for mapping")  
+    
+                        if len(read) == 0:
+                            continue
+                        len_read = len(read.split('\t')[9])        
+                        # if it's read of good length
+                        if len_read == kmer and (tag in read or multi_tag not in read):
+                            feat = read.split('\t')
+                            seq = feat[9]
+                            # if it's a reverse read
+                            if feat[1] == '16' :
+                                if site == "A" :
+                                # #get A-site
+                                    cod = str(Seq(seq[a_pos-5:a_pos-2]).reverse_complement())
+                                elif site == "P" :
+                                # #get P-site
+                                    cod = str(Seq(seq[a_pos-2:a_pos+1]).reverse_complement())
+                                else :
+                                # #get site-E
+                                    cod = str(Seq(seq[a_pos+1:a_pos+4]).reverse_complement())
+                                # # test if it's a true codon not a CNG codon for example
+                                if codon_dict.has_key(cod) :
+                                    codon_dict[cod] += 1
+                            # if it's a forward read
+                            elif feat[1] == '0' :
+                                if site == "A" :
+                                # #get A-site
+                                    cod = seq[a_pos:a_pos+3]
+                                elif site == "P" :
+                                # #get P-site
+                                    cod = seq[a_pos-3:a_pos]
+                                else :
+                                # #get site-E
+                                    cod = seq[a_pos-6:a_pos-3]
+                                if codon_dict.has_key(cod) :
+                                    codon_dict[cod] += 1
+                    del(read)
                 # # add in global dict   
                 for cod, count in codon_dict.iteritems() :
                     codon[cod] += count
-                        
-        return codon
+        if sum(codon.values()) == 0 :
+            stop_err('There are no reads aligning on annotated genes in your GFF file')          
+        else :
+            return codon
         
     except Exception, e:
         stop_err('Error during codon usage calcul: ' + str(e))
@@ -317,12 +306,12 @@
             std_cond2 = []
             max_val = []  # # max value for graph
             for i in codon_sorted:
-                # # cond1 = moyenne of replicats cond1 divided by max 
+                # # cond1 = mean of replicats cond1 divided by max 
                 cond1_val[i] = ((cond1_1[i] / sum11 + cond1_2[i] / sum12) / 2)
                 cond1[i] = ((cond1_1[i] + cond1_2[i]) / 2)
-                # # standard deviation = absolute value of diffence between replicats of cond1
+                # # standard deviation = absolute value of difference between replicats of cond1
                 std_cond1.append(std(array([(cond1_1[i] * 100 / sum11), (cond1_2[i] * 100 / sum12)])))
-                # # cond2 = moyenne of replicats cond1divided by max
+                # # cond2 = mean of replicats cond1divided by max
                 cond2_val[i] = ((cond2_1[i] / sum21 + cond2_2[i] / sum22) / 2)
                 cond2[i] = ((cond2_1[i] + cond2_2[i]) / 2)
                 # # standard deviation = absolute value of difference between replicats of cond2
@@ -388,7 +377,8 @@
                     out.write(i + '\t' + str(cond1[i]) + '\t' + str(cond2[i]) + '\t' + str(cond1_norm[i]) + '\t' + str(cond2_norm[i]) + '\t0.0\n')
                 else: 
                     out.write(i + '\t' + str(cond1[i]) + '\t' + str(cond2[i]) + '\t' + str(cond1_norm[i]) + '\t' + str(cond2_norm[i]) + '\t' + str(cond2_norm[i] / cond1_norm[i]) + '\n')
-            chi = stats.chisquare(observed, expected)
+            with errstate(all='ignore'):
+                chi = stats.chisquare(observed, expected)
             out.write('Khi2 test\n')
             out.write('T : ' + str(chi[0]) + '; p-value : ' + str(chi[1]) + '\n')
 
@@ -435,7 +425,7 @@
             cond1_norm.update ((x, (y / sum1) * 100.0) for x, y in cond1_norm.items())
             cond2_norm.update((x, (y / sum2) * 100.0) for x, y in cond2_norm.items())
         except ZeroDivisionError:
-            stop_err("Not enough reads to compute the codon occupancy")
+            stop_err("Not enough reads to compute the codon occupancy. "+str(sum1)+" and "+str(sum2)+" reads are used for each condition, respectively.\n")
             
         # # compute theorical count in COND2
         cond2_count = []
@@ -497,7 +487,8 @@
                 else:
                     out.write(i + '\t' + str(cond1[i]) + '\t' + str(cond2[i]) + '\t' + str(cond1_norm[i]) + '\t' + str(cond2_norm[i]) + '\t' + str(cond2_norm[i] / cond1_norm[i]) + '\n')
             out.write('Khi2 test\n')
-            chi = stats.chisquare(observed, expected)
+            with errstate(all='ignore'):
+                chi = stats.chisquare(observed, expected)
             out.write('T : ' + str(chi[0]) + '; p-value : ' + str(chi[1]) + '\n')
             
         # # get max value for each codon for histogram   
@@ -603,7 +594,10 @@
     fc = cond1.copy()
 
     for key, value in fc.iteritems():
-        fc[key] = cond2[key]/cond1[key]
+        if cond1[key] == 0:
+            fc[key] = 1
+        else:
+            fc[key] = cond2[key]/cond1[key]
         
     index = arange(len(fc.keys()))
     label = fc.keys()
@@ -695,12 +689,11 @@
         GFF = HTSeq.GFF_Reader(options.gff)
         # # get html file and directory :
         (html, html_dir) = options.dirout.split(',')
-        if os.path.exists(html_dir):
-            raise
-        try:
-            os.mkdir(html_dir)
-        except:
-            raise Exception(html_dir + ' mkdir')
+        if not os.path.exists(html_dir):
+            try:
+                os.mkdir(html_dir)
+            except Exception, e :
+                stop_err('Error running make directory : ' + str(e)) 
         # #RUN analysis
         # #If there are replicats
         if options.rep == "yes" :
@@ -724,6 +717,7 @@
                 check_index_bam (fh)
                 result.append(get_codon_usage(fh, GFF, options.site, options.kmer,options.asite))
             (cond1, cond2, chi_pval) = plot_codon_usage(result, html_dir, options.c1, options.c2, options.outfile,options.color1, options.color2)
+
             # t_pval = compute_FC_plot(cond1,cond2,codons,html_dir)
             plot_fc (cond1, cond2, options.site, html_dir)    
         else :