+ added some text and references to the documentation
authorFabio <fabio@congo.fml.local>
Tue, 16 Sep 2008 20:28:26 +0000 (22:28 +0200)
committerFabio <fabio@congo.fml.local>
Tue, 16 Sep 2008 20:28:26 +0000 (22:28 +0200)
+ added some logic to the last step in the pipeline

doc/qpalma.tex
qpalma/gridtools.py
scripts/createAlignmentFileFromPrediction.py
scripts/qpalma_main.py

index 3e0abbb..35f87df 100644 (file)
 \section{Intro}
 
 \QP is an alignment tool targeted to align spliced reads produced by ``Next
-Generation'' sequencing platforms such as $Illumina Genome Analyzer$, $454$ or
-$SOLid$. We refer to the whole pipeline as the \QP pipeline and \QP
-respectively.
+Generation'' sequencing platforms such as $Illumina Genome Analyzer$ or $454$.
+The basic idea is to use an extended Smith-Waterman algorithm for local
+alignments that uses the base quality information of the reads directly during
+the alignment step. Optimal alignment parameters i.e. scoring matrices are
+inferred using a machine learning technique similar to \emph{Support Vector
+Machines}. For further details on \QP itself consult the paper \cite{DeBona08}.
+For details about the learning method \cite{Tsochantaridis04}.
+%$SOLid$. 
+%We refer to the whole pipeline as the \QP pipeline and \QP respectively.
+
+\section{Quicktour}
 
 Basically \QP assumes that you have the following data:
 \begin{itemize}
@@ -30,9 +38,9 @@ Basically \QP assumes that you have the following data:
 \item splice site scores predicted for the genomic sequences.
 \end{itemize}
 
-The project has two central configuration files: PipelineConf.py and
-QPalmaConf.py where the former stores the pipeline settings and the latter the
-alignment specific options.
+\QP has one central configuration file: 
+
+Suppose you have a
 
 The project results directory (\emph{result\_dir}) contains then the subdirectories
 \begin{itemize}
@@ -152,36 +160,28 @@ The sequence information tuple itself consists of
 \item Stop of the genomic region
 \end{enumerate}
 
-
-
-            currentSeqInfo 
-
 %
 %
 %
-\section{Data Standards}
+\section{Format Specifications}
 
+This section introduces all formats and conventions that are assumed to be met
+by the users in order to make \QP work.
 
-\subsection{Read format and internal representation}
+\subsection{Format of the configuration file}
 
-\begin{itemize}
-\item Which nucleotide bears the score of the splice site ?
-\item What exactly are the exon/intron boundaries pointing to (also are they 0/1-based) ?
-\end{itemize}
+The configuration file includes are settings \QP needs to perform an analysis.
+This includes paths to file where the raw data exists as well as settings which
+sequencing platform is being used,the number of cluster nodes to employ etc. .
 
-We address the first question as follows:
-\texttt{
-----gt-------ag------
-    *         *
-}
-the score sits at the g's of the splice sites.
-    
+Its values are in the form 
+\begin{center}
+key = value
+\end{center}
+and ``#'' for lines containing comments.
 
-%
-%
-%
-\section{Format Specifications}
 
+\subsection{File Formats}
 The format of the file containing the mapped short reads is as follows.  Each
 line corresponds to one short read. Each line has six tab separated entries,
 namely:
@@ -194,9 +194,40 @@ namely:
 \item read quality
 \end{enumerate}
 
+\subsection{Read format and internal representation}
+
+\begin{itemize}
+\item Which nucleotide bears the score of the splice site ?
+\item What exactly are the exon/intron boundaries pointing to (also are they 0/1-based) ?
+\end{itemize}
+
+We address the first question as follows:
+\texttt{
+----gt-------ag------
+    *         *
+}
+the score is assumed to be saved at the position of the g's of the splice sites.
 
 \subsection{Splice Scores}
 
+The splice site scores where generated by ... . However if you train by
+yourself then you can use splice site predictions of any kind... as long as the
+prediction of one site is a real value.
+
+
+\begin{thebibliography}{1}
+
+\bibitem[1]{DeBona08} 
+   De~Bona~F.~and~Ossowski~S.~and~Schneeberger~K.~and~G.~R{\"a}tsch
+   \newblock Optimal Spliced Alignment of Short Sequence Reads
+   \newblock {\em ECCB 2008}
+
+\bibitem[2]{Tsochantaridis04} 
+   Ioannis~Tsochantaridis~and~Thomas~Hofmann~and~Thorsten~Joachims~and~Yasemin~Altun
+   \newblock Support Vector Machine Learning for Interdependent and Sturcutured Output Spaces
+   \newblock {\em Proceedings of the 16th International Conference on Machine Learning}, 2004
+
+\end{thebibliography}
 %
 %
 %
index e3dddcd..912e93f 100644 (file)
@@ -256,28 +256,25 @@ class PostprocessingTask(ClusterTask):
       ClusterTask.__init__(self)
 
 
-   def TaskStarter(chunk_fn,result_fn):
-      create_alignment_file(chunk_fn,result_fn)
-
-
    def CreateJobs(self):
-      run_dir     = '/fml/ag-raetsch/home/fabio/tmp/vmatch_evaluation/spliced_1/prediction'
-      result_dir  = '/fml/ag-raetsch/home/fabio/tmp/vmatch_evaluation/spliced_1/alignment'
+      run_dir     = self.settings['prediction_dir']
+      self.result_dir  = self.settings['alignment_dir']
 
       chunks_fn = []
       for fn in os.listdir(run_dir):
          if fn.startswith('chunk'):
             chunks_fn.append(fn)
 
-      print chunks_fn
-
       functionJobs=[]
 
+      self.result_files = []
       for chunk_fn in chunks_fn:
          chunk_name  = chunk_fn[:chunk_fn.find('.')]
-         result_fn   = jp(result_dir,'%s.align_remap'%chunk_name)
+         result_fn   = jp(self.result_dir,'%s.align'%chunk_name)
          chunk_fn = jp(run_dir,chunk_fn) 
 
+         self.result_files.append(result_fn)
+
          current_job = KybJob(grid_alignment.TaskStarter,[chunk_fn,result_fn])
          current_job.h_vmem = '15.0G'
          current_job.express = 'True'
@@ -285,6 +282,16 @@ class PostprocessingTask(ClusterTask):
          print "job #1: ", current_job.nativeSpecification
 
 
+   def collectResults(self):
+      combined_fn = jp(self.result_dir,'all_alignments.align')
+      combine_files(self.result_files,combined_fn)
+
+
+def TaskStarter(chunk_fn,result_fn):
+   create_alignment_file(chunk_fn,result_fn)
+
+
+
 class DummyTask(ClusterTask):
    """
    This class represents a dummy task to make debugging easier.
index c772b7e..8fdc822 100644 (file)
@@ -24,12 +24,16 @@ from palma.output_formating import print_results
 
 from qpalma.sequence_utils import SeqSpliceInfo,DataAccessWrapper
 
+# a little auxiliary function
+pp = lambda x : str(x)[1:-1].replace(' ','')
+
 
 def alignment_reconstruct(current_prediction,num_exons):
    """
    We reconstruct the exact alignment to infer the positions of indels and the
    sizes of the respective exons.
    """
+
    translation = '-acgtn_' # how aligned est and dna sequence is displayed (gap_char, 5 nucleotides, intron_char) 
 
    SpliceAlign = current_prediction['spliceAlign']
@@ -82,17 +86,8 @@ def create_alignment_file(current_loc,out_fname):
    print 'Loaded %d examples' % len(allPredictions)
 
    # fetch the data needed
-   g_dir    = '/fml/ag-raetsch/share/projects/genomes/A_thaliana_best/genome/'
-   acc_dir  = '/fml/ag-raetsch/home/fabio/tmp/interval_query_files/acc'
-   don_dir  = '/fml/ag-raetsch/home/fabio/tmp/interval_query_files/don'
-
-   g_fmt = 'chr%d.dna.flat'
-   s_fmt = 'contig_%d%s'
-
-   num_chromo = 6
-
-   accessWrapper = DataAccessWrapper(g_dir,acc_dir,don_dir,g_fmt,s_fmt)
-   seqInfo = SeqSpliceInfo(accessWrapper,range(1,num_chromo))
+   accessWrapper = DataAccessWrapper(settings)
+   seqInfo = SeqSpliceInfo(accessWrapper,settings['allowed_fragments'])
 
    # Uniqueness filtering using alignment score for reads with several
    # alignments
@@ -162,10 +157,6 @@ def create_alignment_file(current_loc,out_fname):
             Tgapbs += predExons[1,0] - predExons[0,1]
             Tgapcnt = 1
 
-         # &strand, Qname, &Qsize, 
-         # &Qstart, &Qend, Tname, &Tsize,
-         # &Tstart, &Tend, &blockcnt,exonsizes_str,Qstarts_str,Tstarts_str
-         #
          # ATTENTION
          #  
          # we enforce equal exons sizes for q and t because we are missing indel
@@ -197,36 +188,8 @@ def create_alignment_file(current_loc,out_fname):
          blockcnt,exonsizes_str,Qstarts_str,Tstarts_str) #str(predExons))
 
       else:
-         #try:
          exonIdentities,exonGaps = alignment_reconstruct(current_prediction,num_exons)
-         #except:
-         #   print 'Bug in example %d (idx: %d)' %\
-         #   (current_prediction['id'],pred_idx)
-         #   continue
-
-         #t_size = tEnd - tStart
-         #read_size = 38
-         #if strand == '-':
-         #   total_size = seqInfo.chromo_sizes[chromo]
 
-         #   start_pos = total_size - start_pos
-
-         #   qExonSizes.reverse()
-         #   qStarts  = [read_size-e for e in qEnds]
-         #   qStarts.reverse()
-         #   qEnds  = [read_size-e for e in qStarts]
-         #   qEnds.reverse()
-
-         #   tExonSizes.reverse()
-         #   tStarts  = [t_size-e for e in tEnds]
-         #   tStarts.reverse()
-         #   tEnds    = [t_size-e for e in tStarts]
-         #   tEnds.reverse()
-
-         #   exonIdentities.reverse()
-         #   exonGaps.reverse()
-
-         pp = lambda x : str(x)[1:-1].replace(' ','')
 
          new_line = '%d\t%d\t%s\t%s\t%s\t%d\t%d\t%d\t%d\t%d\t%d\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n' %\
          (id,chromo,strand,seq,str(q1)[1:-1],start_pos,qStart,qEnd,tStart,tEnd,num_exons,\
@@ -241,7 +204,48 @@ def create_alignment_file(current_loc,out_fname):
    print 'Wrote out %d lines' % written_lines
 
 
-if __name__ == '__main__':
-   current_dir = sys.argv[1]
-   out_fh = sys.argv[2]
-   create_alignment_file(current_dir,out_fh)
+def run(chunk_dir,outfile):
+   """   
+   The run...
+   """
+
+   out_fh = open(outfile,'w+')
+
+   # fetch the data needed
+   accessWrapper = DataAccessWrapper(settings)
+   seqInfo = SeqSpliceInfo(accessWrapper,settings['allowed_fragments'])
+
+
+   for line in open(result_fn):
+      sl = line.split()
+
+      chromo = int(sl[1])
+      strand = sl[2]
+      start_pos = int(sl[5])-2
+      starts   = sl[15].split(',')
+      ends     = sl[16].split(',')
+
+      ids = sl[-2].split(',')
+      ids   = [int(elem) for elem in ids]
+      gaps = sl[-1].split(',')
+      gaps   = [int(elem) for elem in gaps]
+
+      if strand == '+':
+         pass
+      elif strand == '-':
+         total_size = seqInfo.chromo_sizes[chromo]
+         start_pos = total_size - start_pos
+      else:
+         assert False
+
+      starts   = [int(elem) + start_pos for elem in starts]
+      ends     = [int(elem) + start_pos for elem in ends]
+
+      if strand == '+':
+         starts   = [e+1 for e in starts]
+         ends     = [e+1 for e in ends]
+      else:
+         starts   = [e-3001 for e in starts]
+         ends     = [e-3001 for e in ends]
+
+      out_fh.write("%d\t%s\t%s\t%s\t%s\t%s\n"%(chromo,strand,pp(starts),pp(ends),pp(ids),pp(gaps)))
index cf07c55..e70e6d2 100644 (file)
@@ -634,7 +634,8 @@ class QPalma:
       _newSpliceAlign   = array.array('B',newSpliceAlign)
       _newEstAlign      = array.array('B',newEstAlign)
        
-      alignment = get_alignment(_newSpliceAlign,_newEstAlign, dna_array, read_array) #(qStart, qEnd, tStart, tEnd, num_exons, qExonSizes, qStarts, qEnds, tExonSizes, tStarts, tEnds)
+      #(qStart, qEnd, tStart, tEnd, num_exons, qExonSizes, qStarts, qEnds, tExonSizes, tStarts, tEnds)
+      alignment = get_alignment(_newSpliceAlign,_newEstAlign, dna_array, read_array) 
 
       dna_array   = array.array('B',dna_array)
       read_array  = array.array('B',read_array)