Merge branch 'HbRNAi'

Author: Peter Combs

MakeSummaryTable.py

@@ -6,7 +6,7 @@
     2) -c -- also include confidence intervals
 
 """
-from __future__ import division
+from __future__ import division, print_function
 import pandas
 from os import path
 from glob import glob
@@ -132,8 +132,10 @@ def get_stagenum(name, series, dir):
             stage=params.ix[dirname]['Stage'],
             num=get_stagenum(dirname, params.index,
                              params.ix[dirname, 'Direction']))
-        print dirname, '=', new_dirname
+        print(dirname, '=', new_dirname)
         dirname = new_dirname
+    else:
+        print(dirname)
 
     skip = False
     if args.strip_low_reads:
@@ -158,16 +160,17 @@ def get_stagenum(name, series, dir):
             pass
         total_reads += i//4
         skip += (reads / total_reads) < (args.strip_low_map_rate / 100)
-        print(reads, total_reads, reads/total_reads,
-              args.strip_low_map_rate / 100)
+        if skip:
+            print(reads, total_reads, reads/total_reads,
+                  args.strip_low_map_rate / 100)
 
     if skip:
         if args.strip_as_nan:
             from numpy import nan
-            print "NaNing", dirname
+            print("NaNing", dirname, "\t{:,} reads".format(reads))
             table.ix[:] = nan
         else:
-            print "Skipping", dirname
+            print("Skipping", dirname)
             continue
     if df is None:
         df = pandas.DataFrame({dirname+"_FPKM": table.ix[:, args.column]})

Makefile

@@ -13,31 +13,32 @@ MELVERSION = $(word 1, $(subst _FB, ,$(MELRELEASE)))
 VIRVERSION = $(word 1, $(subst _FB, ,$(VIRRELEASE)))
 MELDATE = $(word 2, $(subst _FB, ,$(MELRELEASE)))
 
-MELFASTA = prereqs/dmel-all-chromosome-$(MELVERSION).fasta
-VIRFASTA = prereqs/dvir-all-chromosome-$(VIRVERSION).fasta
+PREREQDIR = prereqs
+MELFASTA = $(PREREQDIR)/dmel-all-chromosome-$(MELVERSION).fasta
+VIRFASTA = $(PREREQDIR)/dvir-all-chromosome-$(VIRVERSION).fasta
 
 REFDIR = Reference
 
 MELFASTA2= $(REFDIR)/dmel_prepend.fasta
 VIRFASTA2= $(REFDIR)/dvir_prepend.fasta
 
-ORTHOLOGS = prereqs/gene_orthologs_fb_$(MELDATE).tsv
+ORTHOLOGS = $(PREREQDIR)/gene_orthologs_fb_$(MELDATE).tsv
 
-CERFASTA = prereqs/S288C_reference_sequence_R64-1-1_20110203.fsa
+CERFASTA = $(PREREQDIR)/S288C_reference_sequence_R64-1-1_20110203.fsa
 CERFASTA2= $(REFDIR)/scer_prepend.fasta
 
-MELGFF   = prereqs/dmel-all-$(MELVERSION).gff
+MELGFF   = $(PREREQDIR)/dmel-all-$(MELVERSION).gff
 MELGTF   = $(REFDIR)/mel_good.gtf
-VIRGFF   = prereqs/dvir-all-$(VIRVERSION).gff
+VIRGFF   = $(PREREQDIR)/dvir-all-$(VIRVERSION).gff
 VIRGTF   = $(REFDIR)/vir_good.gtf
-CERGFF   = prereqs/saccharomyces_cerevisiae_R64-1-1_20110208.gff
+CERGFF   = $(PREREQDIR)/saccharomyces_cerevisiae_R64-1-1_20110208.gff
 MELVIRGTF= $(REFDIR)/melvir.gtf
 MELVIRGTF_FILT= $(REFDIR)/melvir_withgenename.gtf
 MELVIRFASTA=$(REFDIR)/melvir.fa
 MELALLGTF   = $(REFDIR)/mel_all.gtf
 MELBADGTF   = $(REFDIR)/mel_bad.gtf
 
-GENEMAPTABLE = gene_map_table_fb_$(MELDATE).tsv
+GENEMAPTABLE = $(PREREQDIR)/gene_map_table_fb_$(MELDATE).tsv
 
 
 all : $(ANALYSIS_DIR)/summary.tsv $(REFDIR)/$(MELMAJORVERSION) $(REFDIR)/$(MELVERSION)
@@ -124,11 +125,11 @@ $(MELBADGTF): $(MELALLGTF) | $(REFDIR)
 		> $@
 
 
-$(MELFASTA): $(REFDIR)/$(MELMAJORVERSION) | $(REFDIR)
+$(MELFASTA): $(REFDIR)/$(MELMAJORVERSION) | $(REFDIR) $(PREREQDIR)
 	wget -O $@.gz ftp://ftp.flybase.net/genomes/Drosophila_melanogaster/dmel_$(MELRELEASE)/fasta/dmel-all-chromosome-$(MELVERSION).fasta.gz
 	gunzip --force $@.gz
 
-$(MELGFF): $(REFDIR)/$(MELVERSION) | $(REFDIR)
+$(MELGFF): $(REFDIR)/$(MELVERSION) | $(REFDIR) $(PREREQDIR)
 	wget -O $@.gz ftp://ftp.flybase.net/genomes/Drosophila_melanogaster/dmel_$(MELRELEASE)/gff/dmel-all-$(MELVERSION).gff.gz
 	gunzip --force $@.gz
 
@@ -179,7 +180,7 @@ $(REFDIR)/DmelDvir/Genome : $(MELVIRGTF) |  $(REFDIR)/DmelDvir $(MELFASTA2) $(VI
 		--genomeFastaFiles $(MELFASTA2) $(VIRFASTA2) \
 		--sjdbGTFfile $(MELVIRGTF)
 
-$(ORTHOLOGS) :
+$(ORTHOLOGS) : | $(PREREQDIR)
 	wget -O $@.gz -i ftp.flybase.org/releases/FB$(MELDATE)/precomputed_files/genes/gene_orthologs_fb_$(MELDATE).tsv.gz
 	gunzip --force $@.gz
 
@@ -190,6 +191,8 @@ $(REFDIR)/DmelDvir:
 
 $(REFDIR)/DmelScer:
 	mkdir $@
+$(PREREQDIR):
+	mkdir $@
 
 Reference/DmelDper:
 	mkdir $@
@@ -200,8 +203,9 @@ Reference/DmelDvir:
 Reference/DmelDmoj:
 	mkdir $@
 
+
 $(GENEMAPTABLE):
-	wget ftp://ftp.flybase.net/releases/$(MELDATE)/precomputed_files/genes/$(GENEMAPTABLE).gz \
+	wget ftp://ftp.flybase.net/releases/FB$(MELDATE)/precomputed_files/genes/$(notdir $(GENEMAPTABLE)).gz \
 		-O $(GENEMAPTABLE).gz
 	gunzip --force $(GENEMAPTABLE).gz
 

PlotUtils.py

@@ -144,6 +144,7 @@ def svg_heatmap(data, filename, row_labels=None, box_size=4,
                 draw_box=False, draw_name=False, data_names=None,
                 max_width=np.inf,
                 spacers=None,
+                hatch_nan=True, hatch_size=20,
                 first_col='', last_col=''):
     """
     Draw heatmap as an SVG file stored in filename
@@ -200,9 +201,9 @@ def svg_heatmap(data, filename, row_labels=None, box_size=4,
     pat = dwg.pattern(id='hatch', insert=(0, 0), size=(25, 25),
                       patternUnits='userSpaceOnUse')
     g = pat.add(dwg.g(style="fill:none; stroke:#B0B0B0; stroke-width:1"))
-    g.add(dwg.path(('M0,0', 'l20,20')))
-    g.add(dwg.path(('M10,0 l10,10'.split())))
-    g.add(dwg.path(('M0,10 l10,10'.split())))
+    g.add(dwg.path(('M0,0', 'l{hatch},{hatch}'.format(hatch=hatch_size))))
+    g.add(dwg.path(('M{hatch2},0 l{hatch2},{hatch2}'.format(hatch2=hatch_size/2).split())))
+    g.add(dwg.path(('M0,{hatch2} l{hatch2},{hatch2}'.format(hatch2=hatch_size/2).split())))
 
     dwg.add(pat)
 
@@ -302,7 +303,7 @@ def svg_heatmap(data, filename, row_labels=None, box_size=4,
                                .format(*[int(255*x) for x in
                                          c_cmap(norm_data.ix[i, j])])))
                 dwg.add(g)
-                if hatch:
+                if hatch_nan and hatch:
                     g.add(dwg.rect((x_start + box_size*j,
                                     y_start + i*box_height),
                                    (box_size, box_height),

Utils.py

@@ -31,3 +31,21 @@ def to_number(dataval):
 
     except ValueError:
         return dataval
+
+
+def contains(string_or_iterable):
+    if isinstance(string_or_iterable, str):
+        return lambda x: string_or_iterable in x
+    else:
+        return lambda x: any(i in x for i in string_or_iterable)
+
+def startswith(string_or_iterable):
+    if not isinstance(string_or_iterable, str):
+        string_or_iterable = tuple(string_or_iterable)
+    return lambda x: x.startswith(string_or_iterable)
+
+def sel_contains(string_or_iterable):
+    return dict(crit=contains(string_or_iterable), axis=1)
+
+def sel_startswith(string_or_iterable):
+    return dict(crit=startswith(string_or_iterable), axis=1)

configure

@@ -39,12 +39,47 @@ star_str = ('$(ANALYSIS_DIR)/{label}/accepted_hits.bam: '
 mappers = {'star': star_str, 'tophat' : tophat_str}
 
 glob_specs = [
+    '{seqdir}/*/*{index_seq}*_{{read}}_*.fastq*',
     '{seqdir}/*{label}{index}/*_{{read}}_*.fastq*',
     # '{seqdir}/*{index}/*_{{read}}_*.fastq*',
     '{seqdir}/*/*{label}*/*_{{read}}_*.fastq*',
     '{seqdir}/*index{index}/*_{{read}}*.fastq*',
 ]
 
+index_converter = {
+    'S501':'TAGATCGC',
+    'S502':'CTCTCTAT',
+    'S503':'TATCCTCT',
+    'S504':'AGAGTAGA',
+    'S505':'GTAAGGAG',
+    'S506':'ACTGCATA',
+    'S507':'AAGGAGTA',
+    'S508':'CTAAGCCT',
+    #'N701':'TCGCCTTA',
+    #'N702':'CTAGTACG',
+    #'N703':'TTCTGCCT',
+    #'N704':'GCTCAGGA',
+    #'N705':'AGGAGTCC',
+    #'N706':'CATGCCTA',
+    #'N707':'GTAGAGAG',
+    #'N708':'CCTCTCTG',
+    #'N709':'AGCGTAGC',
+    #'N710':'CAGCCTCG',
+    #'N711':'TGCCTCTT',
+    #'N712':'TCCTCTAC',
+    'N701':'TAAGGCGA',
+    'N702':'CGTACTAG',
+    'N703':'AGGCAGAA',
+    'N704':'TCCTGAGC',
+    'N705':'GGACTCCT',
+    'N706':'TAGGCATG',
+    'N707':'CTCTCTAC',
+    'N708':'CAGAGAGG',
+    'N709':'GCTACGCT',
+    'N710':'CGAGGCTG',
+    'N711':'AAGAGGCA',
+    'N712':'GTAGAGGA',
+}
 
 def parse_arguments():
     from argparse import ArgumentParser
@@ -70,6 +105,10 @@ def parse_arguments():
     return args
 
 
+def convert_indices(index_spec):
+    for index in index_converter:
+        index_spec = index_spec.replace(index, index_converter[index])
+    return index_spec
 
 args = parse_arguments()
 config_file = pd.read_table(args.parameters,
@@ -99,7 +138,12 @@ for i, row in config_file.iterrows():
     species = row['CarrierSpecies']
     for glob_spec in glob_specs:
         glob_spec = (glob_spec
-                     .format(seqdir=args.seqdir, label=label, index=index, id=mbepc))
+                     .format(seqdir=args.seqdir,
+                             label=label,
+                             index=index,
+                             index_seq=convert_indices(index),
+                             id=mbepc,
+                            ))
         rf1 = glob(glob_spec.format(read='R1'))
         if rf1 != []:
             out.write("GLOBSPEC = {}\n".format(glob_spec))