Merge branch 'BcdRNAi'

Author: Peter Combs

MakeSummaryTable.py

@@ -142,16 +142,17 @@ def get_stagenum(name, series, dir):
             reads = 0
             for read in sf:
                 reads += not read.is_secondary
-                if reads > args.strip_low_reads and not args.strip_low_map_rate:
+                if ((reads > args.strip_low_reads)
+                    and not args.strip_low_map_rate):
                     break
             skip = reads < args.strip_low_reads
         else:
             skip = sf.mapped < args.strip_low_reads
     if args.strip_low_map_rate and args.has_params and not skip:
-        rfs = sorted(glob(path.join('sequence',
-                                    '*{}*'.format(params.ix[old_dirname]['Index']),
-                                    '*_R1_*.fastq.gz'))
-                    )
+        rfs = [entry for entry in
+               sf.header['PG'][0]['CL'].split()
+              if entry.endswith('.gz') or entry.endswith('.fastq')][0]
+        rfs = sorted(rfs.split(','))
         total_reads = 4e6 * (len(rfs) - 1)
         for i, line in enumerate(gzip.open(rfs[-1])):
             pass

Makefile

@@ -6,7 +6,7 @@ STARCONFIG = Parameters/STAR_params.in
 ANALYSIS_DIR = analysis
 
 # Reference FASTA and GFF files from FlyBase and SGD
-MELRELEASE = r6.01_FB2014_04
+MELRELEASE = r6.02_FB2014_05
 VIRRELEASE = r1.2_FB2012_01
 MELMAJORVERSION = $(word 1, $(subst ., , $(MELRELEASE)))
 MELVERSION = $(word 1, $(subst _FB, ,$(MELRELEASE)))
@@ -34,6 +34,8 @@ CERGFF   = prereqs/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
 
@@ -71,8 +73,14 @@ $(ANALYSIS_DIR)/summary.tsv : MakeSummaryTable.py $(FPKMS) $(RUNCONFIG) Makefile
 	@echo 'Calculating Abundances'
 	@echo '============================='
 	touch $@
-	cufflinks --num-threads 8 --output-dir $(@D) -u \
-		--frag-bias-correct $(MELFASTA2) -G $(MELGTF) $<
+	cufflinks \
+		--num-threads 8 \
+		--output-dir $(@D) \
+		--multi-read-correct \
+		--frag-bias-correct $(MELFASTA2) \
+		--GTF $(MELGTF) \
+		--mask-file $(MELBADGTF) \
+		$<
 
 %/accepted_hits_sorted.bam: %/accepted_hits.bam
 	touch $@
@@ -92,11 +100,9 @@ $(ANALYSIS_DIR)/summary.tsv : MakeSummaryTable.py $(FPKMS) $(RUNCONFIG) Makefile
 	rm $(@D)/assigned_dmel_sorted.bam
 	samtools index $@
 
-
-$(MELGTF): $(MELGFF) | $(REFDIR)
+$(MELALLGTF): $(MELGFF) | $(REFDIR)
 	gffread $< -E -T -o- | \
-		awk '{print "dmel_"$$0}' | \
-		grep -vP '(snoRNA|CR[0-9]{4}|Rp[ILS]|mir-|tRNA|unsRNA|snRNA|snmRNA|scaRNA|rRNA|RNA:|mt:)' > \
+		awk '{print "dmel_"$$0}' > \
 		$@
 
 $(VIRGTF): $(VIRGFF) $(ORTHOLOGS) | $(REFDIR)
@@ -107,11 +113,22 @@ $(VIRGTF): $(VIRGFF) $(ORTHOLOGS) | $(REFDIR)
 		| python FilterOrthologs.py $(ORTHOLOGS) \
 		> $@
 
+$(MELGTF): $(MELALLGTF) | $(REFDIR)
+	cat $< \
+		| grep -vP '(snoRNA|CR[0-9]{4}|Rp[ILS]|mir-|tRNA|unsRNA|snRNA|snmRNA|scaRNA|rRNA|RNA:|mt:)' \
+		> $@
+
+$(MELBADGTF): $(MELALLGTF) | $(REFDIR)
+	cat $< \
+		| grep -P '(snoRNA|CR[0-9]{4}|Rp[ILS]|mir-|tRNA|unsRNA|snRNA|snmRNA|scaRNA|rRNA|RNA:|mt:)' \
+		> $@
+
+
 $(MELFASTA): $(REFDIR)/$(MELMAJORVERSION) | $(REFDIR)
 	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)
 	wget -O $@.gz ftp://ftp.flybase.net/genomes/Drosophila_melanogaster/dmel_$(MELRELEASE)/gff/dmel-all-$(MELVERSION).gff.gz
 	gunzip --force $@.gz
 

PlotUtils.py

@@ -200,8 +200,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', 'l25,25')))
-    g.add(dwg.path(('M25,0 l-25,25'.split())))
+    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())))
 
     dwg.add(pat)
 
@@ -246,6 +247,11 @@ def svg_heatmap(data, filename, row_labels=None, box_size=4,
             norm_data = frame.divide(frame.dropna(axis=1).mean(axis=1), axis=0)
         elif normer is 'max':
             norm_data = frame.divide(frame.dropna(axis=1).max(axis=1), axis=0)
+        elif normer is 'center0':
+            norm_data = (0.5 +
+                         0.5 * frame.divide(frame.dropna(axis=1).abs().max(axis=1),
+                                      axis=0)
+                        )
         elif index is not None and hasattr(normer, "ix"):
             norm_data = frame.divide(normer.ix[index], axis=0)
         elif hasattr(normer, "__len__") and len(normer) == rows:
@@ -345,5 +351,5 @@ def svg_heatmap(data, filename, row_labels=None, box_size=4,
     if draw_row_labels:
         for i in range(rows):
             dwg.add(dwg.text(row_labels[i],
-                             (x_start, y_start + i*box_size+box_height),))
+                             (x_start, y_start + i*box_height+box_height),))
     dwg.saveas(filename)

configure

@@ -38,9 +38,12 @@ star_str = ('$(ANALYSIS_DIR)/{label}/accepted_hits.bam: '
 
 mappers = {'star': star_str, 'tophat' : tophat_str}
 
-glob_specs = ['{seqdir}/*/*{label}*/*_{{read}}_*.fastq*',
-              '{seqdir}/*index{index}/*_{{read}}*.fastq*',
-             ]
+glob_specs = [
+    '{seqdir}/*{label}{index}/*_{{read}}_*.fastq*',
+    # '{seqdir}/*{index}/*_{{read}}_*.fastq*',
+    '{seqdir}/*/*{label}*/*_{{read}}_*.fastq*',
+    '{seqdir}/*index{index}/*_{{read}}*.fastq*',
+]
 
 
 def parse_arguments():