File: last-pair-probs.html

package info (click to toggle)
last-align 490-1
links: PTS, VCS
area: main
in suites: jessie, jessie-kfreebsd
size: 2,196 kB
ctags: 2,317
sloc: cpp: 20,908; python: 1,636; ansic: 1,466; makefile: 365; sh: 107
file content (513 lines) | stat: -rw-r--r-- 15,158 bytes
<?xml version="1.0" encoding="utf-8" ?>
<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
<meta name="generator" content="Docutils 0.6: http://docutils.sourceforge.net/" />
<title>last-pair-probs</title>
<style type="text/css">

/*
:Author: David Goodger (goodger@python.org)
:Id: $Id: html4css1.css 7056 2011-06-17 10:50:48Z milde $
:Copyright: This stylesheet has been placed in the public domain.

Default cascading style sheet for the HTML output of Docutils.

See http://docutils.sf.net/docs/howto/html-stylesheets.html for how to
customize this style sheet.
*/

/* used to remove borders from tables and images */
.borderless, table.borderless td, table.borderless th {
  border: 0 }

table.borderless td, table.borderless th {
  /* Override padding for "table.docutils td" with "! important".
     The right padding separates the table cells. */
  padding: 0 0.5em 0 0 ! important }

.first {
  /* Override more specific margin styles with "! important". */
  margin-top: 0 ! important }

.last, .with-subtitle {
  margin-bottom: 0 ! important }

.hidden {
  display: none }

a.toc-backref {
  text-decoration: none ;
  color: black }

blockquote.epigraph {
  margin: 2em 5em ; }

dl.docutils dd {
  margin-bottom: 0.5em }

object[type="image/svg+xml"], object[type="application/x-shockwave-flash"] {
  overflow: hidden;
}

/* Uncomment (and remove this text!) to get bold-faced definition list terms
dl.docutils dt {
  font-weight: bold }
*/

div.abstract {
  margin: 2em 5em }

div.abstract p.topic-title {
  font-weight: bold ;
  text-align: center }

div.admonition, div.attention, div.caution, div.danger, div.error,
div.hint, div.important, div.note, div.tip, div.warning {
  margin: 2em ;
  border: medium outset ;
  padding: 1em }

div.admonition p.admonition-title, div.hint p.admonition-title,
div.important p.admonition-title, div.note p.admonition-title,
div.tip p.admonition-title {
  font-weight: bold ;
  font-family: sans-serif }

div.attention p.admonition-title, div.caution p.admonition-title,
div.danger p.admonition-title, div.error p.admonition-title,
div.warning p.admonition-title {
  color: red ;
  font-weight: bold ;
  font-family: sans-serif }

/* Uncomment (and remove this text!) to get reduced vertical space in
   compound paragraphs.
div.compound .compound-first, div.compound .compound-middle {
  margin-bottom: 0.5em }

div.compound .compound-last, div.compound .compound-middle {
  margin-top: 0.5em }
*/

div.dedication {
  margin: 2em 5em ;
  text-align: center ;
  font-style: italic }

div.dedication p.topic-title {
  font-weight: bold ;
  font-style: normal }

div.figure {
  margin-left: 2em ;
  margin-right: 2em }

div.footer, div.header {
  clear: both;
  font-size: smaller }

div.line-block {
  display: block ;
  margin-top: 1em ;
  margin-bottom: 1em }

div.line-block div.line-block {
  margin-top: 0 ;
  margin-bottom: 0 ;
  margin-left: 1.5em }

div.sidebar {
  margin: 0 0 0.5em 1em ;
  border: medium outset ;
  padding: 1em ;
  background-color: #ffffee ;
  width: 40% ;
  float: right ;
  clear: right }

div.sidebar p.rubric {
  font-family: sans-serif ;
  font-size: medium }

div.system-messages {
  margin: 5em }

div.system-messages h1 {
  color: red }

div.system-message {
  border: medium outset ;
  padding: 1em }

div.system-message p.system-message-title {
  color: red ;
  font-weight: bold }

div.topic {
  margin: 2em }

h1.section-subtitle, h2.section-subtitle, h3.section-subtitle,
h4.section-subtitle, h5.section-subtitle, h6.section-subtitle {
  margin-top: 0.4em }

h1.title {
  text-align: center }

h2.subtitle {
  text-align: center }

hr.docutils {
  width: 75% }

img.align-left, .figure.align-left, object.align-left {
  clear: left ;
  float: left ;
  margin-right: 1em }

img.align-right, .figure.align-right, object.align-right {
  clear: right ;
  float: right ;
  margin-left: 1em }

img.align-center, .figure.align-center, object.align-center {
  display: block;
  margin-left: auto;
  margin-right: auto;
}

.align-left {
  text-align: left }

.align-center {
  clear: both ;
  text-align: center }

.align-right {
  text-align: right }

/* reset inner alignment in figures */
div.align-right {
  text-align: inherit }

/* div.align-center * { */
/*   text-align: left } */

ol.simple, ul.simple {
  margin-bottom: 1em }

ol.arabic {
  list-style: decimal }

ol.loweralpha {
  list-style: lower-alpha }

ol.upperalpha {
  list-style: upper-alpha }

ol.lowerroman {
  list-style: lower-roman }

ol.upperroman {
  list-style: upper-roman }

p.attribution {
  text-align: right ;
  margin-left: 50% }

p.caption {
  font-style: italic }

p.credits {
  font-style: italic ;
  font-size: smaller }

p.label {
  white-space: nowrap }

p.rubric {
  font-weight: bold ;
  font-size: larger ;
  color: maroon ;
  text-align: center }

p.sidebar-title {
  font-family: sans-serif ;
  font-weight: bold ;
  font-size: larger }

p.sidebar-subtitle {
  font-family: sans-serif ;
  font-weight: bold }

p.topic-title {
  font-weight: bold }

pre.address {
  margin-bottom: 0 ;
  margin-top: 0 ;
  font: inherit }

pre.literal-block, pre.doctest-block, pre.math {
  margin-left: 2em ;
  margin-right: 2em }

span.classifier {
  font-family: sans-serif ;
  font-style: oblique }

span.classifier-delimiter {
  font-family: sans-serif ;
  font-weight: bold }

span.interpreted {
  font-family: sans-serif }

span.option {
  white-space: nowrap }

span.pre {
  white-space: pre }

span.problematic {
  color: red }

span.section-subtitle {
  /* font-size relative to parent (h1..h6 element) */
  font-size: 80% }

table.citation {
  border-left: solid 1px gray;
  margin-left: 1px }

table.docinfo {
  margin: 2em 4em }

table.docutils {
  margin-top: 0.5em ;
  margin-bottom: 0.5em }

table.footnote {
  border-left: solid 1px black;
  margin-left: 1px }

table.docutils td, table.docutils th,
table.docinfo td, table.docinfo th {
  padding-left: 0.5em ;
  padding-right: 0.5em ;
  vertical-align: top }

table.docutils th.field-name, table.docinfo th.docinfo-name {
  font-weight: bold ;
  text-align: left ;
  white-space: nowrap ;
  padding-left: 0 }

h1 tt.docutils, h2 tt.docutils, h3 tt.docutils,
h4 tt.docutils, h5 tt.docutils, h6 tt.docutils {
  font-size: 100% }

ul.auto-toc {
  list-style-type: none }
/* Style sheet for LAST HTML documents */
h1 { color: navy }
h2 { color: teal }
div.document { margin-left: auto; margin-right: auto; max-width: 40em }
li { margin-bottom: 1ex }

</style>
</head>
<body>
<div class="document" id="last-pair-probs">
<h1 class="title">last-pair-probs</h1>

<p>This script reads alignments of paired DNA reads to a genome, and:</p>
<ol class="arabic simple">
<li>estimates the distribution of distances between paired reads,</li>
<li>estimates the probability that each alignment represents the
genomic source of the read.</li>
</ol>
<p>The script takes as input two files of alignments, where one read from
each pair is in one file, and the other read is in the other file.</p>
<p>Typical usage:</p>
<pre class="literal-block">
lastdb -m1111110 humandb human/chr*.fa
lastal -Q1 -e120 -i1 humandb reads1.fastq &gt; temp1.maf
lastal -Q1 -e120 -i1 humandb reads2.fastq &gt; temp2.maf
last-pair-probs temp1.maf temp2.maf &gt; results.maf
</pre>
<p>If your reads come from potentially-spliced RNA molecules, use the -r
option:</p>
<pre class="literal-block">
last-pair-probs -r temp1.maf temp2.maf &gt; results.maf
</pre>
<p>Without -r, it assumes the distances between paired reads follow a
normal distribution.  With -r, it assumes the distances follow a
skewed (log-normal) distribution, which is much more appropriate for
spliced RNA.</p>
<div class="section" id="details">
<h2>Details</h2>
<ul class="simple">
<li>The script writes the alignments with &quot;mismap&quot; probabilities,
i.e. the probability that the alignment does not represent the
genomic source of the read.  By default, it discards alignments with
mismap probability &gt; 0.01.</li>
<li>It assumes that each pair of reads comes from opposite strands of a
DNA molecule.  The &quot;distance&quot; between them means the distance
between their 5' ends.  Positive distance indicates tail-to-tail
orientation (with the 5' end being the head and the 3' end being the
tail).  Negative distance indicates head-to-head orientation.
Negative distances are not considered when -r is used, nor for
circular chromosomes.</li>
<li>The script reads one batch of alignments at a time from each file
(by looking for lines starting with &quot;# batch&quot;).  It assumes that
each pair of batches has the alignments for one pair of reads.  To
achieve that, the read pairs should be in the same order in the
fastq files, and the lastal -i1 option ensures there is one query
per batch.</li>
<li>The input files may be in either format produced by lastal (maf or
tabular).  They must include header lines (of the kind produced by
lastal) describing the alignment parameters.</li>
<li>By default, the script makes two passes over each file.  So they
must be real files (not pipes).  However, if you use option -e, or
both -f and -s, it makes just one pass over each file.</li>
<li>If a read name ends in neither &quot;/1&quot; nor &quot;/2&quot;, the script appends
&quot;/1&quot; if it comes from the first file or &quot;/2&quot; if it comes from the
second.</li>
</ul>
</div>
<div class="section" id="options">
<h2>Options</h2>
<blockquote>
<table class="docutils option-list" frame="void" rules="none">
<col class="option" />
<col class="description" />
<tbody valign="top">
<tr><td class="option-group">
<kbd><span class="option">-h</span>, <span class="option">--help</span></kbd></td>
<td>Print a help message and exit.</td></tr>
<tr><td class="option-group">
<kbd><span class="option">-r</span>, <span class="option">--rna</span></kbd></td>
<td>Specifies that the fragments are from potentially-spliced RNA.</td></tr>
<tr><td class="option-group">
<kbd><span class="option">-e</span>, <span class="option">--estdist</span></kbd></td>
<td>Just estimate the distribution of distances.</td></tr>
<tr><td class="option-group" colspan="2">
<kbd><span class="option">-m <var>M</var></span>, <span class="option">--mismap=<var>M</var></span></kbd></td>
</tr>
<tr><td>&nbsp;</td><td>Don't write alignments with mismap probability &gt; M.</td></tr>
<tr><td class="option-group" colspan="2">
<kbd><span class="option">-f <var>BP</var></span>, <span class="option">--fraglen=<var>BP</var></span></kbd></td>
</tr>
<tr><td>&nbsp;</td><td>The mean distance in bp.  (With -r, the mean of
ln[distance].)  If this is not specified, the script will
estimate it from the alignments.</td></tr>
<tr><td class="option-group" colspan="2">
<kbd><span class="option">-s <var>BP</var></span>, <span class="option">--sdev=<var>BP</var></span></kbd></td>
</tr>
<tr><td>&nbsp;</td><td>The standard deviation of distance in bp.  (With -r, the
standard deviation of ln[distance].)  If this is not
specified, the script will estimate it from the alignments.</td></tr>
<tr><td class="option-group" colspan="2">
<kbd><span class="option">-d <var>PROB</var></span>, <span class="option">--disjoint=<var>PROB</var></span></kbd></td>
</tr>
<tr><td>&nbsp;</td><td>The prior probability that a pair of reads comes from
disjoint locations (e.g., different chromosomes).  This may
arise from real differences between the genome and the source
of the reads, or from errors in obtaining the reads or the
genome sequence.</td></tr>
<tr><td class="option-group" colspan="2">
<kbd><span class="option">-c <var>CHROM</var></span>, <span class="option">--circular=<var>CHROM</var></span></kbd></td>
</tr>
<tr><td>&nbsp;</td><td>Specifies that the chromosome named CHROM is circular.  You
can use this option more than once (e.g., -c chrM -c chrP).
As a special case, &quot;.&quot; means all chromosomes are circular.
If this option is not used, &quot;chrM&quot; is assumed to be circular
(but if it is used, only the specified CHROMs are assumed to
be circular.)</td></tr>
</tbody>
</table>
</blockquote>
</div>
<div class="section" id="tips">
<h2>Tips</h2>
<ul>
<li><p class="first">For greater speed and convenience, first estimate the distance
distribution from a sample of your data (using -e), then analyze all
your data with this estimate (using -f and -s).</p>
</li>
<li><p class="first">You can avoid temp files by using named pipes:</p>
<pre class="literal-block">
mkfifo pipe1 pipe2
lastal -Q1 -e120 -i1 humandb reads1.fastq &gt; pipe1 &amp;
lastal -Q1 -e120 -i1 humandb reads2.fastq &gt; pipe2 &amp;
last-pair-probs -f250 -s30 pipe1 pipe2 &gt; results.maf
</pre>
<p>This streams the alignments from lastal to last-pair-probs.</p>
</li>
<li><p class="first">To go faster, try a fast version of Python such as PyPy.</p>
</li>
<li><p class="first">To go really fast, try gapless alignment (add -j1 to the lastal
options).  Often, this is only minusculely less accurate than gapped
alignment.</p>
</li>
<li><p class="first">Tabular output (lastal option -f0) is smaller and faster.</p>
</li>
</ul>
</div>
<div class="section" id="using-multiple-cpus">
<h2>Using multiple CPUs</h2>
<p>With large datasets, it's important to go faster by using multiple
CPUs.  One way to do that is by using GNU parallel
(<a class="reference external" href="http://www.gnu.org/software/parallel/">http://www.gnu.org/software/parallel/</a>), as follows.</p>
<ul>
<li><p class="first">Split reads1.fastq into multiple files, with (say) 100000 reads per
file:</p>
<pre class="literal-block">
split -l400000 -a5 reads1.fastq 1
</pre>
<p>This assumes that a new fastq record begins every 4th line, i.e. no
line wrapping.  The created files will be called 1aaaaa, 1aaaab,
etc.</p>
</li>
<li><p class="first">Split reads2.fastq:</p>
<pre class="literal-block">
split -l400000 -a5 reads2.fastq 2
</pre>
</li>
<li><p class="first">Make a file called (say) last-pair.sh, with the following content
(or similar - you might want to use -r, -j1, etc):</p>
<pre class="literal-block">
#! /bin/sh
lastal -Q1 -e120 -i1 &quot;$3&quot; &quot;$4&quot; &gt; /tmp/$$.1
lastal -Q1 -e120 -i1 &quot;$3&quot; &quot;$5&quot; &gt; /tmp/$$.2
last-pair-probs -f &quot;$1&quot; -s &quot;$2&quot; /tmp/$$.1 /tmp/$$.2
rm /tmp/$$.*
</pre>
</li>
<li><p class="first">Set execute permission:</p>
<pre class="literal-block">
chmod +x last-pair.sh
</pre>
</li>
<li><p class="first">Run it in parallel on all your CPU cores:</p>
<pre class="literal-block">
parallel --xapply ./last-pair.sh 250 30 humandb ::: 1* ::: 2* &gt; results.maf
</pre>
<p>Here we have specified a mean distance of 250 and a standard
deviation of 30.</p>
</li>
</ul>
</div>
<div class="section" id="reference">
<h2>Reference</h2>
<p>For more information, please see this article:</p>
<blockquote>
An approximate Bayesian approach for mapping paired-end DNA reads to
a reference genome.  Shrestha AM, Frith MC.  Bioinformatics 2013
29(8):965-972.</blockquote>
</div>
</div>
</body>
</html>