Bioutils: Difference between revisions

From QiuLab
Jump to navigation Jump to search
imported>Weigang
imported>Weigang
Line 62: Line 62:


==Full documentation==
==Full documentation==
<html xmlns="http://www.w3.org/1999/xhtml">
* bioseq
<head>
* bioaln
<title>bioseq</title>
* biopop
<meta http-equiv="content-type" content="text/html; charset=utf-8" />
* biotree
<link rev="made" href="mailto:root@wallace.hunter.cuny.edu" />
</head>
 
<body style="background-color: white">
 
 
 
<ul id="index">
  <li><a href="#NAME">NAME</a></li>
  <li><a href="#SYNOPSIS">SYNOPSIS</a></li>
  <li><a href="#DESCRIPTION">DESCRIPTION</a></li>
  <li><a href="#OPTIONS">OPTIONS</a></li>
  <li><a href="#EXAMPLES">EXAMPLES</a></li>
  <li><a href="#REQUIRES">REQUIRES</a></li>
  <li><a href="#SEE-ALSO">SEE ALSO</a></li>
  <li><a href="#AUTHORS">AUTHORS</a></li>
  <li><a href="#POD-ERRORS">POD ERRORS</a></li>
</ul>
 
<h1 id="NAME">NAME</h1>
 
<p>bioseq - Fasta sequence editing module based on BioPerl.</p>
 
<h1 id="SYNOPSIS">SYNOPSIS</h1>
 
<p><b>bioseq</b> [options] [sequence file]</p>
 
<h1 id="DESCRIPTION">DESCRIPTION</h1>
 
<p><b>bioseq</b> will read a sequence file and act upon it by doing the following - reformat input (default is fasta) to Genbank or EMBL formats, delete specified sequences, generate overlapping subsequence with a specified window size, generate the reverese complementary sequence, for nucleic acid sequences only, take input list of sequences apart into individual sequence files, extract a specified subset of the sequence, linearize the sequence, remove gaps, find the longest open reading frame (ORF), remove stop codons, give percentage composition of specified amino acids or nuclic acid bases, split the sequences as specified by the user, translate a specific frame of input sequence, or extract a specific gene ID from multiple file sequences. By default, <b>bioseq</b> will assume that both the input and the output are in FASTA format.</p>
 
<h1 id="OPTIONS">OPTIONS</h1>
 
<dl>
 
<dt id="help--h"><b>--help, -h</b></dt>
<dd>
 
<p>Print a brief help message and exit.</p>
 
<pre><code>        Usage: bioseq -h &lt;keyword&gt;</code></pre>
 
</dd>
<dt id="man--m"><b>--man, -m</b></dt>
<dd>
 
<p>Prints the manual page and exit.</p>
 
<pre><code>        Usage: bioseq -m &lt;keyword&gt; </code></pre>
 
</dd>
<dt id="input--i-format"><b>--input, -i</b> &#39;format&#39;</dt>
<dd>
 
<p>Input file format. By default, this is &#39;fasta&#39;. For Genbank format, use &#39;genbank&#39;. For EMBL format, use &#39;embl&#39;.</p>
 
<pre><code>        Usage: bioseq -i &#39;genbank&#39; input_file</code></pre>
 
</dd>
<dt id="output--o-format"><b>--output, -o</b> &#39;format&#39;</dt>
<dd>
 
<p>Output file format. By default, this is &#39;fasta&#39;.For Genbank format, use &#39;genbank&#39;. For EMBL format, use &#39;embl&#39;.</p>
 
<pre><code>        Usage: bioseq -o &#39;EMBL&#39; input_file</code></pre>
 
</dd>
<dt id="noseq--n"><b>--noseq, -n</b></dt>
<dd>
 
<p>Print number of sequences specified by n.</p>
 
<pre><code>        Usage: bioseq -n input_file</code></pre>
 
</dd>
<dt id="sub--s-min-max"><b>--sub, -s</b> &#39;min,max&#39;</dt>
<dd>
 
<p>Select substring (of 1st sequence),</p>
 
<pre><code>        Usage: bioseq -s &#39;&lt;beginning index&gt;, &lt;ending index&gt;&#39; input_file
 
        Example:  bioseq -s&#39;20,80&#39; input_file (or --sub&#39;20,80&#39; or -s=&#39;20,80&#39; or --sub=&#39;20,80&#39;)</code></pre>
 
</dd>
<dt id="lengths--l"><b>--lengths, -l</b></dt>
<dd>
 
<p>Print all sequence lengths.</p>
 
<pre><code>        Usage: bioseq -l</code></pre>
 
</dd>
<dt id="leadgaps--y"><b>--leadgaps, -y</b></dt>
<dd>
 
<p>Count and return the number of leading gaps in each sequence.</p>
 
<pre><code>        Usage: bioseq -y</code></pre>
 
</dd>
<dt id="pick--p-tag:value"><b>--pick, -p</b> &#39;tag:value&#39;</dt>
<dd>
 
<p>Select a single sequence or a comma-separated list of sequences, e.g, --pick &#39;id:foo&#39; by id --pick &#39;order:2&#39; by order --pick &#39;re:REGEX&#39; using a regular expression (only one regex is expected)</p>
 
<p>Specifically for a list of sequences, --pick &#39;id:foo,bar&#39; list by id --pick &#39;order:2,3&#39; list by order --pick &#39;order:2-10&#39; list by range</p>
 
<pre><code>        Usage: bioseq -p &#39;id:&lt;number&gt;&#39; input_file</code></pre>
 
</dd>
<dt id="del--d-tag:value"><b>--del, -d</b> &#39;tag:value&#39;</dt>
<dd>
 
<pre><code> Delete a sequence or a comma-separated list of sequences, eg,
--del &#39;id:foo&#39;        by id
--del &#39;order:2&#39;        by order
--del &#39;length:n&#39;      by min length, where &#39;n&#39; is length
--del &#39;ambig:x&#39;        by min % ambiguous base/aa, where &#39;x&#39; is the %
--del &#39;id:foo,bar&#39;    list by id
--del &#39;re:REGEX&#39;      using a regular expression (only one regex is expected)
 
        Usage: bioseq --del &#39;id:&lt;number&gt;&#39; input_file</code></pre>
 
</dd>
<dt id="revcom--r"><b>--revcom, -r</b></dt>
<dd>
 
<p>Reverse complement</p>
 
<pre><code>        Usage: bioseq -r input_file</code></pre>
 
</dd>
<dt id="slidingwindow--k-window_size"><b>--slidingwindow, -k</b> &#39;window_size&#39;</dt>
<dd>
 
<p>Generate overlapping subsequence with a specified window size (default 1)</p>
 
<pre><code>        Usage: bioseq -k &#39;&lt;index of subsequence&gt;&#39; input_file</code></pre>
 
</dd>
<dt id="translate--t-1-3-6"><b>--translate, -t</b> [1|3|6]</dt>
<dd>
 
<p>Translate in 1, 3, or 6 frames. eg, -t1, -t3, or -t6.</p>
 
<pre><code>        Usage: bioseq -t3 input_file</code></pre>
 
</dd>
<dt id="shred--c"><b>--shred, -c</b></dt>
<dd>
 
<p>Shred into individual sequences</p>
 
<pre><code>        Usage: bioseq -c input_file</code></pre>
 
</dd>
<dt id="linearize--L"><b>--linearize, -L</b></dt>
<dd>
 
<p>Linearize FASTA.</p>
 
<pre><code>        Usage: bioseq -L fasta_file</code></pre>
 
</dd>
<dt id="extract--e"><b>--extract, -e</b></dt>
<dd>
 
<p>Extract in-frame sequences.</p>
 
<pre><code>        Usage: bioseq -e input_file</code></pre>
 
</dd>
<dt id="nogaps--g"><b>--nogaps, -g</b></dt>
<dd>
 
<p>Remove gaps</p>
 
<pre><code>        Usage: bioseq -g input_file</code></pre>
 
</dd>
<dt id="longest-orf--C"><b>--longest-orf, -C</b></dt>
<dd>
 
<p>Output the frame that gives the longest ORF.</p>
 
<pre><code>        Usage: bioseq -C input_file</code></pre>
 
</dd>
<dt id="removestop--x"><b>--removestop, -x</b></dt>
<dd>
 
<p>Remove stop codons (for e.g., PAML input)</p>
 
<pre><code>        Usage: bioseq -x input_file</code></pre>
 
</dd>
<dt id="anonymize--a-n"><b>--anonymize, -a</b> &#39;n&#39;</dt>
<dd>
 
<p>Replace sequence IDs with serial IDs &#39;n&#39; characters long, including a leading &#39;S&#39; (e.g., -a&#39;5&#39; gives S0001). Produces a sed script file with a &#39;.sed&#39; suffix that may be used with sed&#39;s &#39;-f&#39; argument. If the filename is &#39;-&#39;, the sed file is named STDOUT.sed instead. The sed filename is specified on STDERR.</p>
 
<pre><code>        Usage: bioseq -a &lt;number&gt; input_file</code></pre>
 
</dd>
<dt id="prefix-PREFIX"><b>--prefix</b> &#39;PREFIX&#39;</dt>
<dd>
 
<p>Used in conjunction with --anonymize. This lets you specify a custom prefix for the anonymized sequence IDs given when using the --anonymize option. If this is given, then the whole prefix will count toward the total ID length. For example: suppose the prefix chosen is SEQ, and that for --anonymize you supplied 5. Then the maximum id length is 5, so there is room for only two more digits. e.g., SEQ01 atg... SEQ02 atg... SEQ03 atc...</p>
 
<p>If there are enough sequences that the length of the prefix plus the length of the digit portion exceeds the length given to --anonymize, a warning will be given: aln-manipulations.pl -a 4 --prefix=SEQ # output SEQ1 atg... ... SEQ10 atc...</p>
 
<pre><code> # more output
WARNING: Anonymized ID length exceeded requested length: try a different length or prefix.</code></pre>
 
</dd>
<dt id="composition--w"><b>--composition, -w</b></dt>
<dd>
 
<p>Base or AA composition.</p>
 
<pre><code>        Usage: bioseq -w input_file</code></pre>
 
</dd>
<dt id="split--S-split_at"><b>--split, -S</b> &#39;split_at&#39;</dt>
<dd>
 
<p>Split the input sequence file into several smaller files with &#39;split_at&#39; sequences per file. For instance, to split a file into several smaller files with 100000 sequences each, you would run:</p>
 
<pre><code> bioseq -S 100000 seq.fasta</code></pre>
 
<p>The output files will be named with the label &quot;split_N&quot; according to the input file (or with the STDIN prefix if the file is read via standard in), where N denotes the &quot;part&quot; or &quot;split&quot; number.</p>
 
</dd>
<dt id="retrieveseq--z-sequence-retriever-using-GenBank-accession"><b>--retrieveseq, -z</b> &#39;sequence retriever using GenBank accession&#39;</dt>
<dd>
 
<pre><code>    Usage: bioseq -z [Accession]
    Retrieves the sequence from GenBank using the provided GenBank accession. Prints out text in a fasta file.
 
  EXAMPLE:
 
    bioseq -z X83553
 
  OUTPUT:
 
    &gt;X83553 B.garinii (PHei strain) opsC gene.
    ATGAAAAAGAATACATTAAGTGCGATATTAATGACTTTATTTTTATTTATATCTTGTAAT
    AATTCAGGTGGGGATACTGCATCTACTAATCCTGATGAGTCTGCGAAAGGACCTAATCTT
    ATAGAAATAAGCAAAAAAATTACAGATTCTAATGCATTTGTACTGGCTGTGAAAGAAGTT
    GAGGCTTTGATCTCATCTATAGATGAACTTGCTAATAAAGCTATTGGTAAAAAAATAAAT
    CAAAATGGTTTAGATGCTGATGCTAATCACAACGGATCATTGTTAGCAGGAGCCCATGCA
    ATATCAACTCTAATAAAACAAAAAACAGATGGATTGAAAGATCTAGAAGGGTTAAGTAAA
    GAAATTGCAAAGGTGAAGGAATGTTCCGATAAATTTACTAAAAAGCTAACAGATAGTCAT
    GCACAGCTTGGAGCAGTTGGTGGTGCTATTAATGATGATCGTGCAAAAGAAGCTATTTTA
    AAAACACATGGGACTAACGATAAGGGTGCTAAAGAACTTAAAGAGTTATCTGAATCAGTA
    GAAAGCTTGGCAAAAGCAGCTCAAGCAGCATTAGCTAATTCAGTTAAAGAGCTTACAAGT
    CCTGTTGTGGCAGAAAGTCCAAAAAAACCTTAA</code></pre>
 
</dd>
<dt id="dotplot--D-draw_dotplot"><b>--dotplot, -D</b> &#39;draw_dotplot&#39;</dt>
<dd>
 
<p>Extract two sequences from input file and generate a dotplot.</p>
 
<pre><code>    Usage: bioseq -D &#39;id1,id2,window_size,slider&#39; fasta_file
    </code></pre>
 
<p>Id1 and Id2 are extracted with their corresponding sequences. Be sure to use the entire sequence identifer, as this is a whole string match. Window_size corresponds to the number of character you would like to compare (Default window is 10). Slider is the number of windows to compare (Default slider is 10). The sequence corresponding to ID1 will appear on the X axis (row) and ID2 on the Y axis (column). This method will work on both DNA and Amino Acids.</p>
 
<pre><code>  Example:
    Sample Input:
    &gt;id1
    ATACGA
    &gt;id2
    ATACGA
 
    Command: bioseq -D &#39;id1,id2,3&#39;
    Output:
 
        A T A C G A
    A  *
    T    *
    A      *
    C
    G
    A
    A
    C
    A
    A
    T
    G</code></pre>
 
<p>==item <b>--rename, -R</b> &#39;rename_id&#39;</p>
 
<pre><code>    Usage: bioseq -R [file_with_new_names] file_to_be_changed.fasta
 
    Ex: bioseq -D list.txt test-bioseq.nuc
 
  Input:
 
    list.txt:
 
    VS116:7:310:IGS:11      VS116
    B31:1:100:IGS:11        B31
 
    *Left column is the pattern to be replaced by the right column
 
    file_to_be_changed.fasta:
 
    &gt;VS116:7:310:IGS:11
    AATTTCAAAAATATAATATAAAAACAGCTAATCCAATAGAAAAATTTGAAATTTTTCTAT
    TGGATAAATTCTATACAAGAAGGTAAATA
    &gt;B31:1:100:IGS:11
    AATTTTTAAAATATAATATAAAAACAGCTAATCCAATAGAAAAATTTTAAAACTTTTCTA
    TTGGATAGATTTTATACAAAGAAGGTAATA
 
  Output:
    &gt;VS116
    AATTTCAAAAATATAATATAAAAACAGCTAATCCAATAGAAAAATTTGAAATTTTTCTAT
    TGGATAAATTCTATACAAGAAGGTAAATA
    &gt;B31
    AATTTTTAAAATATAATATAAAAACAGCTAATCCAATAGAAAAATTTTAAAACTTTTCTA
    TTGGATAGATTTTATACAAAGAAGGTAATA</code></pre>
 
</dd>
</dl>
 
<h1 id="EXAMPLES">EXAMPLES</h1>
 
<p>Section under construction...</p>
 
<h1 id="REQUIRES">REQUIRES</h1>
 
<p>Perl 5.010, BioPerl</p>
 
<h1 id="SEE-ALSO">SEE ALSO</h1>
 
<pre><code>  perl(1)</code></pre>
 
<h1 id="AUTHORS">AUTHORS</h1>
 
<pre><code> Weigang Qiu at genectr.hunter.cuny.edu
Y&ouml;zen Hern&aacute;ndez yzhernand at gmail dot com
Levy Vargas levy dot vargas at gmail dot com</code></pre>
 
<h1 id="POD-ERRORS">POD ERRORS</h1>
 
<p>Hey! <b>The above document had some coding errors, which are explained below:</b></p>
 
<dl>
 
<dt id="Around-line-382">Around line 382:</dt>
<dd>
 
<p>Non-ASCII character seen before =encoding in &#39;Y&ouml;zen&#39;. Assuming UTF-8</p>
 
</dd>
</dl>
 
 
</body>
 
</html>
----
----


Revision as of 16:41, 29 September 2014

BioPerl-based Sequence Utilities

Figure 1. Design and Methods of bioutils

What is bioutils?

bioutils is a suite of Perl scripts that provide convenient command-line access to popular BioPerl methods. Designed as UNIX utilities, these tools aim to circumvent a constant need (and urge) for composing one-off BioPerl scripts for routine manipulations of sequences, alignments and trees.

The initial release of bioutils consists of four utilities (Figure 1):

  • bioseq: a wrapper of BioPerl class Bio::Seq (with additional methods)
  • bioaln: a wrapper of Bio::SimpleAlign (which inherits Bio::Seq; with additional methods)
  • biopop: a wrapper of Bio::PopGen (which can be converted from Bio::SimpleAlign; with additional methods)
  • biotree: a wrapper of Bio::tree (with additional methods)

These utilities have been in development since 2002 in the lab of Dr Weigang Qiu at Hunter College of the City University of New York. They are the main code base of the Qiu Lab, which specializes in microbial evolutionary genomics. They proved to be convenient, efficient, and popular among students and researchers passing through the lab. By releasing bioutils as an Open Source tool, we hope to (1) share our experience and (2) invite other developers to join the effort of making BioPerl more accessible.

Live Demos

Basic Usage

  • bioseq
bioseq -l foo.fasta # print seq names and lengths from FASTA (default format) file
bioseq -r foo.fasta # reverse complement
bioseq -t1 foo.fasta # translate in the +1 frame
bioseq -t3 foo.fasta # translate in +1, +2, and +3 frames
bioseq -t6 foo.fasta # translate in all 6 frames
bioseq -p'id:seq_1' foo.fasta # pick a sequence by ID
bioseq -p'order:3' # pick the 3rd sequence
bioseq -p're:Human' foo.fasta # pick all sequences labeled as "Human" (by regular expression)
bioseq -g foo.fasta # remove all gaps
bioseq -z'CP003201' -o'genbank' # retrieve a GenBank file with accession
bioseq -z'CP003201' -o'fasta' # same file in FASTA
  • bioaln
bioaln -i'fasta' -o'phylip' foo.fasta # convert a FASTA alignment to PHYLIP
bioaln -l foo.aln # print alignment length of a CLUSTALW (default format) file
bioaln -s'100, 200' foo.aln # obtain an alignment slice
bioaln -m foo.aln # show only variable sites
bioaln -r'seq_2' foo.aln # use "seq_2" as reference (first) sequence
bioaln -g foo.aln # remove gapped sites
bioaln -p'seq_1,seq_3,seq_6' foo.aln # pick a subset of sequences
bioaln -d'seq_1,seq_3,seq_6' foo.aln # delele a subset of sequences
  • biotree
  • biopop

Power usage

# Pipe with the same utility
bioseq -p'order:5' foo.fasta | bioseq -s'100,200' | bioseq -r | bioseq -t1 # pick, subseq, revcom, and translate
# Pipe among utilities
bioaln -o'fasta' foo.aln | bioseq -g # remove gaps within individual sequences

Creative usage

echo -ne ">lookup\nATG\n" | bioseq -t1 # Lookup a codon product
len=$(bioaln -l foo.aln); len_degap=$(bioaln -g foo.aln | bioaln -l); echo "$len-$len_degap" | bc -l # count alignment gaps

Full documentation

  • bioseq
  • bioaln
  • biopop
  • biotree

Release notes

Main contributors

  • Yozen Hernandez
  • Weigang Qiu
  • Pedro Pagan
  • Levy Vargas
Retrieved from "https://wiki.genometracker.org/index.php?title=Bioutils&oldid=2969"