--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dexseq-hts_1.0/README	Tue Oct 08 08:22:45 2013 -0400
@@ -0,0 +1,52 @@
+-----------------------------------------
+A Galaxy wrapper for DEXSeq version 1.6.0
+-----------------------------------------
+
+DEXSeq is focused on finding differential exon usage using RNA-seq exon counts between samples with different experimental designs. 
+DEXSeq: http://www.bioconductor.org/packages/release/bioc/html/DEXSeq.html
+
+Requirements:
+-------------
+    Python, HTSeq package :- Preprocessing of sequencing reads and GFF files
+    R, Bio-conductor package :- Required for DEXSeq 
+    SAMTOOLS :- Sequencing read processing 
+
+Contents:
+---------
+
+./bin 
+    Contains a config file for running DEXSeq in different settings, setup_dexseq-hts.sh
+    will help to create the config file according to your local work station. 
+
+./setup_dexseq-hts.sh
+    Setup script, which helps to set the right path to the config file present in the bin
+    folder. 
+    
+./README 
+    Basic information about the wrapper 
+
+./galaxy 
+    Galaxy tool configuration file can be found here. 
+
+./src
+    All relevant scripts for DEXSeq execution, Please follow the shell script to 
+    understand the details.
+
+Licenses:
+---------
+
+    If DEXSeq is used to obtain results for scientific publications it should be cited as [1].
+
+    This wrapper program is free software; you can redistribute it and/or modify it under the 
+    terms of the GNU General Public License as published by the Free Software Foundation; either 
+    version 3 of the License, or (at your option) any later version.
+
+    Copyright(C) 2013 cBio Memorial Sloan Kettering Cancer Center, New York City, USA.
+
+References
+----------
+    [1] Simon Anders and Alejandro Reyes and Wolfgang Huber (2012): Detecting differential usage of exons from RNA-seq data. 
+
+Contact:
+--------
+    support [at] oqtans.org

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dexseq-hts_1.0/bin/dexseq_config.sh	Tue Oct 08 08:22:45 2013 -0400
@@ -0,0 +1,6 @@
+#!/bin/bash 
+export DEXSEQ_VERSION=1.6.0
+export SAMTOOLS_DIR=
+export PYTHON_PATH=/usr/bin/python
+export PYTHONPATH=
+export R_PATH=

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dexseq-hts_1.0/bin/dexseq_config.sh.sample	Tue Oct 08 08:22:45 2013 -0400
@@ -0,0 +1,6 @@
+#!/bin/bash 
+export DEXSEQ_VERSION=1.6.0
+export SAMTOOLS_DIR=/home/vipin/samtools-0.18-beta
+export PYTHON_PATH=/usr/bin/python
+export PYTHONPATH=/home/vipin/lib/python2.7/site-packages/:$PYTHONPATH
+export R_PATH=/home/vipin/bin/R

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dexseq-hts_1.0/galaxy/DEXSeq.xml	Tue Oct 08 08:22:45 2013 -0400
@@ -0,0 +1,117 @@
+<tool id="dexseq-hts" name="DEXSeq" version="1.6.0">
+    <description> Inference of differential exon usage in RNA-Seq</description>
+    <command interpreter="bash"> 
+    ./../src/DEXseq-hts.sh 
+    $anno_input
+    $sPaired
+    $library_type
+    $minQUAL
+    $dexseq_out $dexseq_out.extra_files_path
+#for $i in $replicate_groups
+#for $j in $i.replicates
+$j.bam_alignment:#slurp
+#end for
+
+#end for
+    >> $Log_File </command>
+    <inputs>
+    
+    <!-- Genome annotation in GFF/GTF/GFF3 file --> 
+	<param format="gff,gtf,gff3" name="anno_input" type="data" label="Reference genome annotation file" help="A tab delimited format for storing sequence features and annotations (GFF/GTF/GFF3)."/>
+
+    <!-- RNA-Seq read alignment files-->  
+    <repeat name="replicate_groups" title="Replicate group" min="2">
+        <repeat name="replicates" title="Replicate">
+        <param format="bam" name="bam_alignment" type="data" label="BAM file of aligned RNA-Seq reads" help="Can be generated from SAM file using the SAMTools."/> 
+        </repeat>
+    </repeat>
+
+    <!--Paired -->
+    <param name="sPaired" type="select" display="radio" label="Is this library mate-paired?" help="Indicates whether the data is paired-end (default: single-end)">
+        <option value="no" selected="true">Single-end</option>
+        <option value="yes">Paired-end</option>
+    </param>
+
+    <!-- RNA-Seq library type-->
+    <param name="library_type" type="select" display="radio" label="Strand-specific library type" help="Indicates RNA-Seq library preparation protocol (default: strand-specific assay)">
+        <option value="yes">Stranded</option>
+        <option value="no">Non stranded</option>
+        <option value="reverse">Reverse</option>
+    </param>
+
+    <!-- Read alignment quality-->
+    <param name="minQUAL" type="integer" min="0" value="10" label="Minimum alignment quality" help="Skip all reads with alignment quality lower than the given minimum value (default: 10)"/>
+  </inputs>
+
+  <outputs>
+    <data format="txt" name="dexseq_out" label="${tool.name} on ${on_string}: Differential Expression"/>
+    <data format="txt" name="Log_File" label="${tool.name} on ${on_string}: Log"/>
+  </outputs>
+
+  <tests>
+    <test> 
+    </test>
+  </tests> 
+
+  <help>
+
+.. class:: infomark
+
+**What it does** 
+
+DEXSeq_ is focused on finding differential exon usage using RNA-seq exon counts between samples with different experimental designs. 
+
+.. _DEXSeq: http://www.bioconductor.org/packages/release/bioc/html/DEXSeq.html
+
+`DEXSeq` requires:
+
+Genome annotation in GFF file type, containing the necessary information about the transcripts that are to be quantified.
+
+The BAM alignment files grouped into replicate groups, each containing several replicates. BAM files store the read alignments, The program will also work with only two groups containing only a single replicate each. However, this analysis has less statistical power and is therefore not recommended!
+
+------
+
+**Licenses**
+
+If **DEXSeq** is used to obtain results for scientific publications it
+should be cited as [1]_.
+
+**References** 
+
+.. [1] Simon Anders and Alejandro Reyes and Wolfgang Huber (2012): `Detecting differential usage of exons from RNA-seq data`_. 
+
+.. _Detecting differential usage of exons from RNA-seq data: http://genome.cshlp.org/content/early/2012/06/21/gr.133744.111.full.pdf+html
+
+------
+
+.. class:: infomark
+
+**About formats**
+
+**GFF/GTF format** General Feature Format/Gene Transfer Format is a format for describing genes and other features associated with DNA, RNA and protein sequences. GFF3 lines have nine tab-separated fields:
+
+1. seqid - The name of a chromosome or scaffold.
+2. source - The program that generated this feature.
+3. type - The name of this type of feature. Some examples of standard feature types are "gene", "CDS", "protein", "mRNA", and "exon". 
+4. start - The starting position of the feature in the sequence. The first base is numbered 1.
+5. stop - The ending position of the feature (inclusive).
+6. score - A score between 0 and 1000. If there is no score value, enter ".".
+7. strand - Valid entries include '+', '-', or '.' (for don't know/care).
+8. phase - If the feature is a coding exon, frame should be a number between 0-2 that represents the reading frame of the first base. If the feature is not a coding exon, the value should be '.'.
+9. attributes - All lines with the same group are linked together into a single item.
+
+For more information see http://www.sequenceontology.org/gff3.shtml
+
+**BAM format** The Sequence Alignment/Map (SAM) format is a
+tab-limited text format that stores large nucleotide sequence
+alignments. BAM is the binary version of a SAM file that allows for
+fast and intensive data processing. The format specification and the
+description of SAMtools can be found on
+http://samtools.sourceforge.net/.
+
+------
+
+DEXSeq-hts Wrapper Version 0.1 (Sept. 2013)
+
+</help>
+</tool>

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dexseq-hts_1.0/setup_dexseq-hts.sh	Tue Oct 08 08:22:45 2013 -0400
@@ -0,0 +1,73 @@
+#!/bin/bash
+set -e 
+
+DIR=`dirname $0`
+. ${DIR}/./bin/dexseq_config.sh
+
+echo ==========================================
+echo  DEXSeq-hts setup script \(DEXSeq version $DEXSEQ_VERSION\) 
+echo ==========================================
+echo
+echo SAMTools directory \(currently set to \"$SAMTOOLS_DIR\", system version used if left empty\)
+read SAMTOOLS_DIR
+if [ "$SAMTOOLS_DIR" == "" ];
+then
+	if [ "$(which samtools)" != "" ] ;
+	then
+		SAMTOOLS_DIR=$(dirname $(which samtools)) 
+	else
+		echo samtools not found
+		exit -1 ;
+	fi
+fi
+echo '=>' Setting SAMTools directory to \"$SAMTOOLS_DIR\"
+echo
+
+echo Path to the python binary \(currently set to \"$PYTHON_PATH\", system version used, if left empty\)
+read PYTHON_PATH
+if [ "$PYTHON_PATH" == "" ];
+then
+    PYTHON_PATH=`which python`
+	if [ "$PYTHON_PATH" == "" ];
+	then
+		echo python not found
+		exit -1 
+	fi
+fi
+echo '=>' Setting Python path to \"$PYTHON_PATH\"
+echo
+
+echo Path to HTSeq library files \(currently set to \"$PYTHONPATH\", system version is used if left empty\)
+read PYTHONPATH
+echo '=>' Setting HTSeq path to \"$PYTHONPATH\"
+echo
+
+echo Path to the R binary \(currently set to \"$R_PATH\", system version used, if left empty\)
+read R_PATH
+if [ "$R_PATH" == "" ];
+then
+    R_PATH=`which R`
+	if [ "$R_PATH" == "" ];
+	then
+		echo R not found
+		exit -1 
+	fi
+fi
+echo '=>' Setting R path to \"$R_PATH\"
+echo
+
+cp -p bin/dexseq_config.sh bin/dexseq_config.sh.bk
+grep -v -e SAMTOOLS_DIR -e PYTHON_PATH -e PYTHONPATH -e R_PATH -e $DEXSEQ_VERSION bin/dexseq_config.sh.bk > bin/dexseq_config.sh
+echo
+echo
+echo generating config file
+
+echo export DEXSEQ_VERSION=$DEXSEQ_VERSION >> bin/dexseq_config.sh
+echo export SAMTOOLS_DIR=$SAMTOOLS_DIR >> bin/dexseq_config.sh
+echo export PYTHON_PATH=$PYTHON_PATH >> bin/dexseq_config.sh
+echo export PYTHONPATH=$PYTHONPATH >> bin/dexseq_config.sh
+echo export R_PATH=$R_PATH >> bin/dexseq_config.sh
+
+echo
+echo Done.
+echo 

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dexseq-hts_1.0/src/DEXseq-hts.sh	Tue Oct 08 08:22:45 2013 -0400
@@ -0,0 +1,98 @@
+#/bin/bash
+##
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 3 of the License, or
+# (at your option) any later version.
+#
+# Copyright (C) 2013 Memorial Sloan-Kettering Cancer Center 
+##
+
+set -e 
+
+PROG=`basename $0`
+DIR=`dirname $0`
+
+. ${DIR}/../bin/dexseq_config.sh
+
+echo
+echo ${PROG}: Oqtans http://galaxy.cbio.mskcc.org Galaxy wrapper for the DEXSeq version ${DEXSEQ_VERSION}.
+echo
+echo DEXSeq: Detecting differential usage of exons from RNA-seq data.
+echo 
+
+## input arguments from the interface 
+GFF_IN=${1}
+shift
+MATE_PAIR=${1}
+shift
+LIBTP=${1}
+shift
+minQL=${1}
+shift
+RES_FILE=${1}
+shift
+RES_WD=${1}
+shift
+
+## associated array with sequencing type.
+declare -A SEQ_TYPE=( [no]=SE [yes]=PE )
+
+echo %%%%%%%%%%%%%%%%%%%%%%%
+echo % 1. Data preparation %
+echo %%%%%%%%%%%%%%%%%%%%%%%
+echo
+
+mkdir -p ${RES_WD}
+echo extra file path $RES_WD
+tmpGTF=`mktemp --tmpdir=/tmp`
+
+echo load the genome annotation in GFF file 
+
+${PYTHON_PATH} ${DIR}/dexseq_prepare_annotation.py ${GFF_IN} ${tmpGTF}
+echo genome annotation stored in ${tmpGTF} 
+echo
+
+echo %%%%%%%%%%%%%%%%%%%%
+echo % 2. Read counting %
+echo %%%%%%%%%%%%%%%%%%%%
+echo
+
+tmpFILE=`mktemp --tmpdir=/tmp`
+echo $tmpFILE
+echo -e '\t'condition'\t'libType > ${tmpFILE}_CONDITIONS.tab
+
+COND=0
+for REPLICATE_GROUP in $@
+do
+    IFS=':'
+    COND=$((COND+1))
+    for BAM_FILE in ${REPLICATE_GROUP}
+    do
+        ## different group information 
+        REPNAME=$(basename ${BAM_FILE%.dat})
+        echo -e ${REPNAME}"\t"$COND"\t"${SEQ_TYPE[$MATE_PAIR]} >> ${tmpFILE}_CONDITIONS.tab
+        
+        ## counting the reads 
+        ${SAMTOOLS_DIR}/samtools view -h $BAM_FILE | ${PYTHON_PATH} ${DIR}/dexseq_count.py -p ${MATE_PAIR} -s ${LIBTP} -a ${minQL} ${tmpGTF} - ${RES_WD}/${REPNAME}
+
+        echo 
+    done
+    echo conuted condition ${COND} 
+done
+echo counted reads map to each exon.
+echo 
+
+echo %%%%%%%%%%%%%%%%%%%%%%%%%%%
+echo % 3. Differential testing %
+echo %%%%%%%%%%%%%%%%%%%%%%%%%%%
+echo
+
+echo "cat ${DIR}/run_DEXseq.R | $R_PATH --slave --args $tmpFILE $RES_WD $tmpGTF ${RES_FILE} $#" 
+cat ${DIR}/run_DEXseq.R | $R_PATH --slave --args $tmpFILE $RES_WD $tmpGTF ${RES_FILE} 
+
+## clean up
+rm -fr ${RES_WD} ${tmpGTF} ${tmpFILE} 
+echo %%%%%%%%
+echo % Done %
+echo %%%%%%%%

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dexseq-hts_1.0/src/dexseq_count.py	Tue Oct 08 08:22:45 2013 -0400
@@ -0,0 +1,173 @@
+import sys, itertools, optparse
+
+optParser = optparse.OptionParser( 
+   
+   usage = "python %prog [options] <flattened_gff_file> <sam_file> <output_file>",
+   
+   description=
+      "This script counts how many reads in <sam_file> fall onto each exonic " +
+      "part given in <flattened_gff_file> and outputs a list of counts in " +
+      "<output_file>, for further analysis with the DEXSeq Bioconductor package. " +
+      "(Notes: The <flattened_gff_file> should be produced with the script " +
+      "dexseq_prepare_annotation.py). <sam_file> may be '-' to indicate standard input.",
+      
+   epilog = 
+      "Written by Simon Anders (sanders@fs.tum.de), European Molecular Biology " +
+      "Laboratory (EMBL). (c) 2010. Released under the terms of the GNU General " +
+      "Public License v3. Part of the 'DEXSeq' package." )
+      
+optParser.add_option( "-p", "--paired", type="choice", dest="paired",
+   choices = ( "no", "yes" ), default = "no",
+   help = "'yes' or 'no'. Indicates whether the data is paired-end (default: no)" )
+
+optParser.add_option( "-s", "--stranded", type="choice", dest="stranded",
+   choices = ( "yes", "no", "reverse" ), default = "yes",
+   help = "'yes', 'no', or 'reverse'. Indicates whether the data is " +
+      "from a strand-specific assay (default: yes ). " +
+      "Be sure to switch to 'no' if you use a non strand-specific RNA-Seq library " +
+      "preparation protocol. 'reverse' inverts strands and is neede for certain " +
+      "protocols, e.g. paired-end with circularization."  )
+   
+optParser.add_option( "-a", "--minaqual", type="int", dest="minaqual",
+   default = 10,
+   help = "skip all reads with alignment quality lower than the given " +
+      "minimum value (default: 10)" )
+   
+if len( sys.argv ) == 1:
+   optParser.print_help()
+   sys.exit(1)
+
+(opts, args) = optParser.parse_args()
+
+if len( args ) != 3:
+   sys.stderr.write( sys.argv[0] + ": Error: Please provide three arguments.\n" )
+   sys.stderr.write( "  Call with '-h' to get usage information.\n" )
+   sys.exit( 1 )
+
+try:
+   import HTSeq
+except ImportError:
+   sys.stderr.write( "Could not import HTSeq. Please install the HTSeq Python framework\n" )   
+   sys.stderr.write( "available from http://www-huber.embl.de/users/anders/HTSeq\n" )   
+   sys.exit(1)
+
+gff_file = args[0]
+sam_file = args[1]
+out_file = args[2]
+stranded = opts.stranded == "yes" or opts.stranded == "reverse"
+reverse = opts.stranded == "reverse"
+is_PE = opts.paired == "yes"
+minaqual = opts.minaqual
+
+if sam_file == "-":
+   sam_file = sys.stdin
+
+# Step 1: Read in the GFF file as generated by aggregate_genes.py
+# and put everything into a GenomicArrayOfSets
+
+features = HTSeq.GenomicArrayOfSets( "auto", stranded=stranded )     
+for f in  HTSeq.GFF_Reader( gff_file ):
+   if f.type == "exonic_part":
+      f.name = f.attr['gene_id'] + ":" + f.attr['exonic_part_number']
+      features[f.iv] += f
+
+# initialise counters
+num_reads = 0
+counts = {}
+counts[ '_empty' ] = 0
+counts[ '_ambiguous' ] = 0
+counts[ '_lowaqual' ] = 0
+counts[ '_notaligned' ] = 0
+
+# put a zero for each feature ID
+for iv, s in features.steps():
+   for f in s:
+      counts[ f.name ] = 0
+
+#We need this little helper below:
+def reverse_strand( s ):
+   if s == "+":
+      return "-"
+   elif s == "-":
+      return "+"
+   else:
+      raise SystemError, "illegal strand"
+
+# Now go through the aligned reads
+
+if not is_PE:
+
+   num_reads = 0
+   for a in HTSeq.SAM_Reader( sam_file ):
+      if not a.aligned:
+         counts[ '_notaligned' ] += 1
+         continue
+      if a.aQual < minaqual:
+         counts[ '_lowaqual' ] += 1
+         continue
+      rs = set()
+      for cigop in a.cigar:
+         if cigop.type != "M":
+            continue
+         if reverse:
+            cigop.ref_iv.strand = reverse_strand( cigop.ref_iv.strand )
+         for iv, s in features[cigop.ref_iv].steps( ):
+            rs = rs.union( s )
+      set_of_gene_names = set( [ f.name.split(":")[0] for f in rs ] )
+      if len( set_of_gene_names ) == 0:
+         counts[ '_empty' ] += 1
+      elif len( set_of_gene_names ) > 1:
+         counts[ '_ambiguous' ] +=1
+      else:
+         for f in rs:
+            counts[ f.name ] += 1
+      num_reads += 1
+      if num_reads % 100000 == 0:
+         sys.stdout.write( "%d reads processed.\n" % num_reads )
+
+else: # paired-end
+
+   num_reads = 0
+   for af, ar in HTSeq.pair_SAM_alignments( HTSeq.SAM_Reader( sam_file ) ):
+      rs = set()
+      if af and ar and not af.aligned and not ar.aligned:
+         counts[ '_notaligned' ] += 1
+         continue
+      if af and ar and not af.aQual < minaqual and ar.aQual < minaqual:
+         counts[ '_lowaqual' ] += 1
+         continue
+      if af and af.aligned and af.aQual >= minaqual and af.iv.chrom in features.chrom_vectors.keys():
+         for cigop in af.cigar:
+            if cigop.type != "M":
+               continue
+            if reverse:
+               cigop.ref_iv.strand = reverse_strand( cigop.ref_iv.strand )
+            for iv, s in features[cigop.ref_iv].steps():
+               rs = rs.union( s )
+      if ar and ar.aligned and ar.aQual >= minaqual and ar.iv.chrom in features.chrom_vectors.keys():
+         for cigop in ar.cigar:
+            if cigop.type != "M":
+               continue
+            if not reverse:
+               cigop.ref_iv.strand = reverse_strand( cigop.ref_iv.strand )
+            for iv, s in features[cigop.ref_iv].steps():
+                  rs = rs.union( s )
+      set_of_gene_names = set( [ f.name.split(":")[0] for f in rs ] )
+      if len( set_of_gene_names ) == 0:
+         counts[ '_empty' ] += 1
+      elif len( set_of_gene_names ) > 1:
+         counts[ '_ambiguous' ] = 0
+      else:
+         for f in rs:
+            counts[ f.name ] += 1
+      num_reads += 1
+      if num_reads % 100000 == 0:
+         sys.stderr.write( "%d reads processed.\n" % num_reads )
+
+ 
+# Step 3: Write out the results
+
+fout = open( out_file, "w" )
+for fn in sorted( counts.keys() ):
+   fout.write( "%s\t%d\n" % ( fn, counts[fn] ) )
+fout.close()

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dexseq-hts_1.0/src/dexseq_prepare_annotation.py	Tue Oct 08 08:22:45 2013 -0400
@@ -0,0 +1,131 @@
+import sys, collections, itertools, os.path
+
+if len( sys.argv ) != 3:
+   sys.stderr.write( "Script to prepare annotation for DEXSeq.\n\n" )
+   sys.stderr.write( "Usage: python %s <in.gtf> <out.gff>\n\n" % os.path.basename(sys.argv[0]) )
+   sys.stderr.write( "This script takes an annotation file in Ensembl GTF format\n" )
+   sys.stderr.write( "and outputs a 'flattened' annotation file suitable for use\n" )
+   sys.stderr.write( "with the count_in_exons.py script.\n" )
+   sys.exit(1)
+
+try:
+   import HTSeq
+except ImportError:
+   sys.stderr.write( "Could not import HTSeq. Please install the HTSeq Python framework\n" )   
+   sys.stderr.write( "available from http://www-huber.embl.de/users/anders/HTSeq\n" )   
+   sys.exit(1)
+
+gtf_file = sys.argv[1]
+out_file = sys.argv[2]
+
+## make sure that it can handle GFF files.
+parent_child_map = dict()
+for feature in HTSeq.GFF_Reader( gtf_file ):
+   if feature.type in ['mRNA', 
+      'transcript', 
+      'ncRNA', 
+      'miRNA', 
+      'pseudogenic_transcript', 
+      'rRNA', 
+      'snoRNA', 
+      'snRNA', 
+      'tRNA', 
+      'scRNA']:
+      parent_child_map[feature.attr['ID']] = feature.attr['Parent']
+
+# Step 1: Store all exons with their gene and transcript ID 
+# in a GenomicArrayOfSets
+
+exons = HTSeq.GenomicArrayOfSets( "auto", stranded=True )
+for f in HTSeq.GFF_Reader( gtf_file ):
+   if not f.type in ["exon", "pseudogenic_exon"]:
+      continue
+   if not f.attr.get('gene_id'):
+      f.attr['gene_id'] = parent_child_map[f.attr['Parent']]
+      f.attr['transcript_id'] = f.attr['Parent']
+   f.attr['gene_id'] = f.attr['gene_id'].replace( ":", "_" )
+   exons[f.iv] += ( f.attr['gene_id'], f.attr['transcript_id'] )
+
+# Step 2: Form sets of overlapping genes
+
+# We produce the dict 'gene_sets', whose values are sets of gene IDs. Each set
+# contains IDs of genes that overlap, i.e., share bases (on the same strand).
+# The keys of 'gene_sets' are the IDs of all genes, and each key refers to
+# the set that contains the gene.
+# Each gene set forms an 'aggregate gene'.
+
+gene_sets = collections.defaultdict( lambda: set() )
+for iv, s in exons.steps():
+   # For each step, make a set, 'full_set' of all the gene IDs occuring
+   # in the present step, and also add all those gene IDs, whch have been
+   # seen earlier to co-occur with each of the currently present gene IDs.
+   full_set = set()
+   for gene_id, transcript_id in s:
+      full_set.add( gene_id )
+      full_set |= gene_sets[ gene_id ]
+   # Make sure that all genes that are now in full_set get associated
+   # with full_set, i.e., get to know about their new partners
+   for gene_id in full_set:
+      assert gene_sets[ gene_id ] <= full_set
+      gene_sets[ gene_id ] = full_set
+
+
+# Step 3: Go through the steps again to get the exonic sections. Each step
+# becomes an 'exonic part'. The exonic part is associated with an
+# aggregate gene, i.e., a gene set as determined in the previous step, 
+# and a transcript set, containing all transcripts that occur in the step.
+# The results are stored in the dict 'aggregates', which contains, for each
+# aggregate ID, a list of all its exonic_part features.
+
+aggregates = collections.defaultdict( lambda: list() )
+for iv, s in exons.steps( ):
+   # Skip empty steps
+   if len(s) == 0:
+      continue
+   # Take one of the gene IDs, find the others via gene sets, and
+   # form the aggregate ID from all of them   
+   gene_id = list(s)[0][0]
+   assert set( gene_id for gene_id, transcript_id in s ) <= gene_sets[ gene_id ] 
+   aggregate_id = '+'.join( gene_sets[ gene_id ] )
+   # Make the feature and store it in 'aggregates'
+   f = HTSeq.GenomicFeature( aggregate_id, "exonic_part", iv )   
+   f.source = os.path.basename( sys.argv[1] )
+   f.attr = {}
+   f.attr[ 'gene_id' ] = aggregate_id
+   transcript_set = set( ( transcript_id for gene_id, transcript_id in s ) )
+   f.attr[ 'transcripts' ] = '+'.join( transcript_set )
+   aggregates[ aggregate_id ].append( f )
+
+
+# Step 4: For each aggregate, number the exonic parts
+
+aggregate_features = []
+for l in aggregates.values():
+   for i in xrange( len(l)-1 ):
+      assert l[i].name == l[i+1].name, str(l[i+1]) + " has wrong name"
+      assert l[i].iv.end <= l[i+1].iv.start, str(l[i+1]) + " starts too early"
+      if l[i].iv.chrom != l[i+1].iv.chrom:
+         raise ValueError, "Same name found on two chromosomes: %s, %s" % ( str(l[i]), str(l[i+1]) )
+      if l[i].iv.strand != l[i+1].iv.strand:
+         raise ValueError, "Same name found on two strands: %s, %s" % ( str(l[i]), str(l[i+1]) )
+   aggr_feat = HTSeq.GenomicFeature( l[0].name, "aggregate_gene", 
+      HTSeq.GenomicInterval( l[0].iv.chrom, l[0].iv.start, 
+         l[-1].iv.end, l[0].iv.strand ) )
+   aggr_feat.source = os.path.basename( sys.argv[1] )
+   aggr_feat.attr = { 'gene_id': aggr_feat.name }
+   for i in xrange( len(l) ):
+      l[i].attr['exonic_part_number'] = "%03d" % ( i+1 )
+   aggregate_features.append( aggr_feat )
+      
+      
+# Step 5: Sort the aggregates, then write everything out
+
+aggregate_features.sort( key = lambda f: ( f.iv.chrom, f.iv.start ) )
+
+fout = open( out_file, "w" ) 
+for aggr_feat in aggregate_features:
+   fout.write( aggr_feat.get_gff_line() )
+   for f in aggregates[ aggr_feat.name ]:
+      fout.write( f.get_gff_line() )
+
+fout.close()      

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dexseq-hts_1.0/src/run_DEXseq.R	Tue Oct 08 08:22:45 2013 -0400
@@ -0,0 +1,73 @@
+### load DEXSeq package
+suppressMessages(require("DEXSeq"))
+
+### get arguments 1: INFILE, 2: OUTFILE 3:SIZE
+args <- commandArgs()
+INFILE<-args[4]
+EXTRAPATH<-args[5]
+annodb<-args[6]
+OUTFILE<-args[7]
+
+INFILE_CONDS=c(paste(INFILE, "_CONDITIONS.tab", sep=""))
+
+### read count data from file
+condsTable <- read.delim( INFILE_CONDS, header=TRUE, stringsAsFactors=TRUE, row.names=1 )
+#head(condsTable)
+
+conditions<-factor( condsTable[ , 1] )
+#print(conditions)
+
+## unique condition to define the pair of tests 
+uniq_conds <- unique(conditions)
+#print(uniq_conds)
+
+## all possible pairs of conditions
+pw_tests <- list()
+for(i in 1:(length(uniq_conds)-1)) {
+    for(j in (i+1):length(uniq_conds)) {
+        pw_tests[[length(pw_tests)+1]] <- c(uniq_conds[i], uniq_conds[j])
+    }
+}
+#print(pw_tests)
+
+tab <- NULL
+## testing all possible pairs of conditions
+for(i in 1:length(pw_tests)) {
+    ## header name 
+    test_pair_name <- c(paste(pw_tests[[i]][1], "__vs__", pw_tests[[i]][2], sep=""))
+    print(test_pair_name)
+
+    sub.data <- subset(condsTable, (conditions %in% c(pw_tests[[i]][1],pw_tests[[i]][2])))
+    sub.data[[1]]<-as.factor(sub.data[[1]])
+
+    ecs = read.HTSeqCounts(countfiles=file.path(EXTRAPATH, row.names(sub.data)), design=sub.data, flattenedfile=annodb)
+    
+    ## Normalisation
+    ecs <- estimateSizeFactors(ecs)
+    print(sizeFactors(ecs))
+
+    suppressMessages(require("parallel"))
+    ## Dispersion estimation 
+    ecs <- estimateDispersions(ecs, nCores=2)
+    ecs <- fitDispersionFunction(ecs)
+
+    ## Testing for differential exon usage
+    ecs <- testForDEU(ecs, nCores=2)
+
+    ## estimate the fold change 
+    ecs <- estimatelog2FoldChanges(ecs, nCores=2)
+
+    ## Result table
+    resultTable <- DEUresultTable(ecs)
+    print(head(resultTable))
+
+    colnames(resultTable) <- paste(test_pair_name, colnames(resultTable), sep=":")
+    #print(colnames(resultTable))
+
+    if(is.null(tab)) {
+        tab<- resultTable
+    }
+    else tab<- cbind(tab, resultTable)
+}
+## printing the result to out file 
+write.table(tab, OUTFILE, quote=F, sep="\t", row.names=F)