PEP8-ify many files

Lots of whitespace errors, line length issues, etc.

Author: Peter Combs

AssignReads2.py

@@ -5,25 +5,31 @@
 from progressbar import ProgressBar, ETA, Bar, Percentage
 from argparse import Namespace
 
+
 class my_defaultdict(dict):
     def __init__(self, default_factory, basename, other_args):
         self.default_factory = default_factory
         self.basename = basename
         self.other_args = other_args
+
     def __missing__(self, key):
         self[key] = value = self.default_factory(self.basename % key,
                                                  **self.other_args)
         return value
 
+
 def get(read, tag):
     return {tname.upper(): val for tname, val in read.tags}[tag.upper()]
 
+
 def get_nh(read):
-    return {tag.upper(): val for tag,val in read.tags}['NH']
+    return {tag.upper(): val for tag, val in read.tags}['NH']
+
 
 def get_species(read):
     return references[read.rname].split('_')[0]
 
+
 def process_read(read):
     if not read.tags:
         print read
@@ -32,13 +38,14 @@ def process_read(read):
     nh = get_nh(read)
     species = get_species(read)
     if nh == 1:
-        #assigned.write(read)
+        # assigned.write(read)
         specific_files[species].write(read)
         species_counts[species] += 1
         return
     else:
         resolve_multiread(read, nh, species)
 
+
 def resolve_multiread(read, nh, species):
     nm = get(read, 'NM')
     has_multi_frags = bool(0x1 & read.flag)
@@ -70,12 +77,13 @@ def resolve_multiread(read, nh, species):
             on_last_multiread(dbs, read)
     else:
         pass
-        #print "didi we not have multiple frags?"
-        #print has_multi_frags
-        #print read.is_read1, read.is_read2
-        #print read.qname
-        #assert False
-        # WTF are we doign here?
+        # print "didi we not have multiple frags?"
+        # print has_multi_frags
+        # print read.is_read1, read.is_read2
+        # print read.qname
+        # assert False
+        #  WTF are we doign here?
+
 
 def on_last_multiread(dbs, read):
     # Sort out the reads
@@ -84,19 +92,19 @@ def on_last_multiread(dbs, read):
         # Report the best, or if equal quality, the first (which
         # tophat would've given anyways)
         species = get_species(read)
-        #assigned.write(read)
+        # assigned.write(read)
         specific_files[species].write(read)
         species_counts[species] += 1
     else:
         # Hits from multiple species
         vals = sorted([(val, spec) for spec, val in
-                        dbs.to_be_resolved_vals[read.qname].iteritems()])
+                       dbs.to_be_resolved_vals[read.qname].iteritems()])
         diff_val = vals[1][0] - vals[0][0]
         ambig_counts[diff_val] += 1
         if diff_val > ambig_threshold:
             species = vals[0][1]
             best_read = dbs.to_be_resolved_reads[read.qname][species]
-            #assigned.write(best_read)
+            # assigned.write(best_read)
             specific_files[species].write(best_read)
             species_counts[species] += 1
         else:
@@ -110,7 +118,7 @@ def on_last_multiread(dbs, read):
                     ambig_names.append(spec)
 
             ambig_names = tuple(ambig_names)
-            ambig_types[ambig_names]+=1
+            ambig_types[ambig_names] += 1
             for amb_read in dbs.to_be_resolved_reads[read.qname].itervalues():
                 ambig.write(amb_read)
 
@@ -126,31 +134,31 @@ def on_last_multiread(dbs, read):
     samfile = pysam.Samfile(fname, 'rb')
     references = samfile.references
     dir = path.dirname(fname)
-    #assigned = pysam.Samfile(path.join(dir, 'assigned.bam'), 'wb',
-                             #template=samfile)
+    # assigned = pysam.Samfile(path.join(dir, 'assigned.bam'), 'wb',
+    # template=samfile)
     ambig = pysam.Samfile(path.join(dir, 'ambiguous.bam'), 'wb',
-                             template=samfile)
+                          template=samfile)
     specific_files = my_defaultdict(pysam.Samfile,
                                     path.join(dir, 'assigned_%s.bam'),
                                     {'template': samfile,
                                      'mode': 'wb'})
 
     to_be_resolved_reads = defaultdict(dict)
-    to_be_resolved_vals = defaultdict(lambda : defaultdict(lambda : 1000))
+    to_be_resolved_vals = defaultdict(lambda: defaultdict(lambda: 1000))
     to_be_resolved_counts = Counter()
     to_be_resolved_reads2 = defaultdict(dict)
-    to_be_resolved_vals2 = defaultdict(lambda : defaultdict(lambda : 1000))
+    to_be_resolved_vals2 = defaultdict(lambda: defaultdict(lambda: 1000))
     to_be_resolved_counts2 = Counter()
     species_counts = Counter()
     ambig_counts = Counter()
     ambig_types = Counter()
 
     print "Measuring file size"
     start = samfile.tell()
-    maxval = path.getsize(fname) * 2**16 # I don't know why it's off by 2^16
+    maxval = path.getsize(fname) * 2**16  # I don't know why it's off by 2^16
     pbar = ProgressBar(maxval=maxval - start + 2**16,
-                       widgets = [fname, ': ', Percentage(), ' ', Bar(), ' ',
-                                  ETA(), ' '])
+                       widgets=[fname, ': ', Percentage(), ' ', Bar(), ' ',
+                                ETA(), ' '])
     pbar.start()
 
     for read in samfile:
@@ -161,5 +169,3 @@ def on_last_multiread(dbs, read):
     print "Species assignments in %s: %s" % (fname, species_counts)
     print "Ambiguity distribution: ", ambig_counts
     print "Ambiguity types: ", ambig_types.most_common(50)
-
-

CalculateAbundances.py

@@ -19,36 +19,37 @@
 cuffcmp = ('cuffcompare -o {cuffname} -s {fasta} -CG -r {gtf} {gtf}')
 cuffname = path.join(path.dirname(gtf_ref), 'cuffcmp')
 
-runstr = cuffcmp.format(cuffname = cuffname,
-                        fasta = fasta_ref,
-                        gtf = gtf_ref)
+runstr = cuffcmp.format(cuffname=cuffname,
+                        fasta=fasta_ref,
+                        gtf=gtf_ref)
 print runstr
 stdout.flush()
 call(runstr.split())
 
-gtf_ref =  cuffname + '.combined.gtf'
+gtf_ref = cuffname + '.combined.gtf'
 
 design_file = pd.read_table(design_fname)
 files = defaultdict(list)
 for ix, row in design_file.iterrows():
     files[row['condition']].append(path.join(analysis_dir,
                                              row['Sample'],
-                                             bamfile_base)
-                                  )
+                                             bamfile_base))
 conditions = sorted(files.keys())
 
-cd = cd_base.format(conditions = ','.join(conditions),
-                    outdir = analysis_dir,
-                    fasta = fasta_ref,
-                    gtf = gtf_ref,
-                    bams = ' '.join([','.join(files[key]) for key in conditions]))
+cd = cd_base.format(conditions=','.join(conditions),
+                    outdir=analysis_dir,
+                    fasta=fasta_ref,
+                    gtf=gtf_ref,
+                    bams=' '.join([','.join(files[key])
+                                   for key in conditions]))
 
 conds_nobcd = [c for c in conditions if 'Bcd' not in c]
-cd_nobcd = cd_base.format(conditions = ','.join(conds_nobcd),
-                    outdir = analysis_dir+'-nobcd',
-                    fasta = fasta_ref,
-                    gtf = gtf_ref,
-                    bams = ' '.join([','.join(files[key]) for key in conds_nobcd]))
+cd_nobcd = cd_base.format(conditions=','.join(conds_nobcd),
+                          outdir=analysis_dir+'-nobcd',
+                          fasta=fasta_ref,
+                          gtf=gtf_ref,
+                          bams=' '.join([','.join(files[key])
+                                         for key in conds_nobcd]))
 print cd_nobcd
 stdout.flush()
 call(cd_nobcd.split())
@@ -57,10 +58,10 @@
 call(cd.split())
 for condition in files:
     for file in files[condition]:
-        cl = cl_base.format(outdir = path.dirname(file),
-                           fasta = fasta_ref,
-                           gtf = gtf_ref,
-                           bamfile = file)
+        cl = cl_base.format(outdir=path.dirname(file),
+                            fasta=fasta_ref,
+                            gtf=gtf_ref,
+                            bamfile=file)
         print '-'*30
         print cl
         call(cl.split())

CalculateConcentrations.py

@@ -1,16 +1,17 @@
 from __future__ import division
 import pandas
 
+
 def main():
-    df = pandas.read_table('CountConfig.tab', index_col = 5)
+    df = pandas.read_table('CountConfig.tab', index_col=5)
     df.dropna(how='any')
 
     for rowname, row in df.iterrows():
         mel_reads = row['mel_reads']
         carrier_reads = row['carrier_reads']
         carrier_conc = row['carrier_conc']
 
-        mel_conc = carrier_conc * mel_reads / carrier_reads 
+        mel_conc = carrier_conc * mel_reads / carrier_reads
         print rowname,
         print mel_conc, 'ng total RNA'
     return df

CheckCoverage.py

@@ -14,7 +14,7 @@
     print """Usage: python CheckCoverage.py <GTF-File> BAMfile [BAMfile ...]"""
     sys.exit(1)
 
-gtf_fname = sys.argv[1] #'Reference/dmel-all-r5.42.gtf'
+gtf_fname = sys.argv[1]  # 'Reference/dmel-all-r5.42.gtf'
 analysis_dir = 'analysis'
 
 starts = set()
@@ -23,11 +23,12 @@
 
 cutoff = 0
 
+
 def analyze_bamfile(bam_fname):
     print bam_fname,
     bam_file = pysam.Samfile(bam_fname, 'rb')
 
-    coverages = defaultdict(lambda : [0,0, set()])
+    coverages = defaultdict(lambda: [0, 0, set()])
     parent = ''
     curr_len = 0
     coverage = 0
@@ -36,28 +37,30 @@ def analyze_bamfile(bam_fname):
     f = open(gtf_fname)
     for line in f:
         pass
-    pb = pbar.ProgressBar(widgets=[bam_fname, pbar.Bar(), pbar.ETA()],maxval=f.tell()).start()
+    pb = pbar.ProgressBar(widgets=[bam_fname, pbar.Bar(), pbar.ETA()],
+                          maxval=f.tell()).start()
     f.seek(0)
     for line in f:
         pb.update(f.tell())
-        if line.startswith('#'): continue
-        if line.startswith('>'): break
+        if line.startswith('#'):
+            continue
+        if line.startswith('>'):
+            break
         data = line.split()
         chrom = data[0]
         kind = data[2]
         start = int(data[3]) - 1
         stop = int(data[4])
         fbtr_finder = re.compile('FBtr[0-9]*')
-        #parent = fbtr_finder.findall(line)[0]
-
-
+        # parent = fbtr_finder.findall(line)[0]
 
         if kind == 'exon':
             fbtrs = fbtr_finder.findall(line)
-            if not fbtrs: continue
+            if not fbtrs:
+                continue
             fbtr = fbtrs[0]
-            coverages[fbtr][1] += (stop - start )
-            
+            coverages[fbtr][1] += (stop - start)
+
             starts = set()
             coverage = 0
             for read in bam_file.fetch(chrom, start, stop):
@@ -66,34 +69,34 @@ def analyze_bamfile(bam_fname):
             coverages[fbtr][0] += coverage
             coverages[fbtr][2].update(starts)
 
-
     pb.finish()
     curr_lens, rpks, uniques = zip(*coverages.itervalues())
     dir, fname = path.split(bam_fname)
 
-
     return dir, rpks, uniques, curr_lens
 
 if __name__ == "__main__":
     import multiprocessing as mp
 
     POOL = mp.Pool(20)
-    res = POOL.map(analyze_bamfile, [f for f in sys.argv[2:] if f.endswith('.bam')])
+    res = POOL.map(analyze_bamfile,
+                   [f for f in sys.argv[2:] if f.endswith('.bam')])
     all_dirs, all_rpks, all_pct_uniques, all_lens = zip(*res)
 
     import cPickle as pickle
     out_fh = open('checkcoverage.pkl', 'w')
-    pickle.dump({'dirs':all_dirs, 'rpks':all_rpks, 'pct_uniques':all_pct_uniques,
-                 'lens': all_lens}, out_fh)
+    pickle.dump({'dirs': all_dirs, 'rpks': all_rpks,
+                 'pct_uniques': all_pct_uniques, 'lens': all_lens},
+                out_fh)
     for fname, rpks, uniques, curr_lens in res:
         print fname
         try:
             xs = array(rpks)
-            ys = array([len(u)/(curr_len + 1) 
+            ys = array([len(u)/(curr_len + 1)
                        for u, curr_len in zip(uniques, curr_lens)])
-            cutoff = max(xs[ys<.1])
-            reg = stats.linregress(log(xs[(xs < cutoff) * (xs > 0) * (ys > 0)]),
-                                   log(ys[(xs < cutoff) * (xs > 0) * (ys > 0)]))
+            cutoff = max(xs[ys < .1])
+            reg = stats.linregress(log(xs[(xs < cutoff)*(xs > 0)*(ys > 0)]),
+                                   log(ys[(xs < cutoff)*(xs > 0)*(ys > 0)]))
             print "exp(%f) * x ** %f" % (reg[1], reg[0])
             print "Duplicate badness score: ", exp(-reg[1]-.38)
         except Exception as exc:

CountAllReads.py

@@ -18,4 +18,3 @@
         print dirname, '\t{:15,}'.format(int(n)), ('*' if n % 1 else '')*10
     except:
         print dirname, "ERR!"
-

GetSpeciesFromBlast.py

@@ -20,8 +20,10 @@
 
 blast_recs5 = [r for r in NCBIXML.parse(open('5.blastout.xml'))]
 blast_recs6 = [r for r in NCBIXML.parse(open('6.blastout.xml'))]
-c5 = cs.Counter([tuple(r.alignments[0].hit_def.split()[:2]) for r in blast_recs5])
-c6 = cs.Counter([tuple(r.alignments[0].hit_def.split()[:2]) for r in blast_recs6])
+c5 = cs.Counter([tuple(r.alignments[0].hit_def.split()[:2])
+                 for r in blast_recs5])
+c6 = cs.Counter([tuple(r.alignments[0].hit_def.split()[:2])
+                 for r in blast_recs6])
 
 
 print(c5)