6
$\begingroup$

I have two DNA strings:

GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC

and

AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG

I want a tool that allows me to do something like this on the command line:

$ aligner  GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG

and receive an ASCII visualization of a pairwise alignment. Something like this would work:

GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC---------
---------AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG

I prefer a tool that is packaged on bioconda. A tool that also does protein and RNA sequences is even more preferable.

$\endgroup$
  • $\begingroup$ Is this a toy example? If so, maybe update it to another one. These sequences are low complexity and the actual correct alignment output is unclear to me $\endgroup$ – Chris_Rands Nov 26 '18 at 16:07
  • $\begingroup$ The sequences are not important. Edits are welcome. $\endgroup$ – winni2k Nov 26 '18 at 16:11
  • 1
    $\begingroup$ To get the alignment in your example, you have to use zero (or very small) gap cost at the ends of sequences. That is not the standard NW. $\endgroup$ – user172818 Nov 27 '18 at 3:43
  • $\begingroup$ I'm not sure that's true? To quote from the EMBOSS package manual: "In a Needleman-Wunsch global alignment, the entire length of each sequence is aligned. The sequences might be partially overlapping or one sequence might be aligned entirely internally to the other. There is no penalty for the hanging ends of the overlap. In bioinformatics, it is usually reasonable to assume that the sequences are incomplete and there should be no penalty for failing to align the missing bases." $\endgroup$ – winni2k Apr 15 at 10:33
11
$\begingroup$

You are looking for the needle program from the EMBOSS suite. Available in bioconda.

http://emboss.sourceforge.net/apps/release/6.6/emboss/apps/needle.html

To read sequences from the commandline, you need so specify the format as asis. To get the output on the screen, you'll need -stdout and to use the default alignment parameters (gap penalty 10, extend_penalty 0.5) you'll need -auto. Thus your query above would be:

$ needle -asequence asis:GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC -bsequence asis:AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG -auto -stdout
########################################
# Program: needle
# Rundate: Tue 27 Nov 2018 10:48:50
# Commandline: needle
#    -asequence asis:GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC
#    -bsequence asis:AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG
#    -auto
#    -stdout
# Align_format: srspair
# Report_file: stdout
########################################

#=======================================
#
# Aligned_sequences: 2
# 1: asis
# 2: asis
# Matrix: EDNAFULL
# Gap_penalty: 10.0
# Extend_penalty: 0.5
#
# Length: 56
# Identity:      37/56 (66.1%)
# Similarity:    37/56 (66.1%)
# Gaps:          18/56 (32.1%)
# Score: 181.0
# 
#
#=======================================

asis               1 GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC---     47
                              |||||||||||||||||||||||||||||||||||||.   
asis               1 ---------AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAG     41

asis              47 ------     47

asis              42 AGGAGG     47


#---------------------------------------
#---------------------------------------
$\endgroup$
9
$\begingroup$

Install biopython with conda and use the following script:

#!/usr/bin/env python
import sys
from Bio import Align
aligner = Align.PairwiseAligner()
aligner.mode = "local"
alignments = aligner.align(sys.argv[1], sys.argv[2])
print(alignments[0])

The output is then:

.GGA-GGAGGGAG--AAGGAGGGAGGGAAGA-GGAGGGAG--AAGGAGGGAGGC
.|-|-||||||||--|||-|||-|||||-||-||||||||--|||-|||-|||.
AG-AAGGAGGGAGGGAAG-AGG-AGGGA-GAAGGAGGGAGGGAAG-AGG-AGG.

There are a LOT of alignments with that same score and you can tweak the details of the mode by changing aligner.mode.

$\endgroup$
  • $\begingroup$ yes, this should be pretty fast too since the alignment work is done in C (rather than python) $\endgroup$ – Chris_Rands Nov 26 '18 at 15:53
4
$\begingroup$

I'd just use something like t_coffee and write a little bash function to run it on the strings you give:

aligner(){
  tmp=$(mktemp);
  args=("$@")
  for ((i=0; i<$#; i++)); do
    printf '>%s\n%s\n' $i ${args[i]} >> "$tmp"
  done
  t_coffee "$tmp"  2>/dev/null | /bin/grep -Pv '^CLUSTAL|^\s*$'
  rm "$tmp"

}

If you add that to your shell's initialization file (e.g. ~/.bashrc if using bash on Linux, ~/.profile if using bash on macOS), you can run align2seqs with as many sequences as you want to align:

$ aligner 'GGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC' 'AGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG'
0               GG--AGGAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC
1               nAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGG-GAAGAGGAGG
                    *********      *******  *********  * ** *  * 
$ aligner 'GAGGGAGAAGGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGC' 'GGAGGGAGGGAAGAGGAGGGAGAAGGAGGGAGGGAAGAGGAGG' GAGGGAGGGAAGAGGAGGGAGAAGGATCCGAGGAAGAGGAGG AGAAGGAGGGAGGGAAGAGGCCCGAGATATATAAGGATCCGAGGAA
0               GA-GG----GAGAAGGAGGGAGGGAAGAGGAGGGAGAAGG--AGGGAGGC
1               GGAGG----GAGGGAAGAGGAGGGAGAAGGAGGGAG-GGA--AGAGGAGG
2               GA-GG----GAGGGAAGAGGAGGGAGAAGGATCCGA-GGA--AGAGGAGG
3               AGAAGGAGGGAGGGAAGAGGCCCGAGATATAT--AA-GGATCCGAG-GAA
                    *    ***      **   **     *       *    * *    
$\endgroup$
3
$\begingroup$

Here's another option: https://github.com/noporpoise/seq-align; I've used it in the past as it is small and has no dependencies (just make and you're set).

$\endgroup$
  • $\begingroup$ This tool actually prompted my question! Do you know if anyone has packaged it? $\endgroup$ – winni2k Nov 29 '18 at 8:23
  • $\begingroup$ Ah, no... I didn't notice the last part of your question. $\endgroup$ – Kirill G Nov 29 '18 at 21:27
  • $\begingroup$ Could you edit and include an example of how the OP could use this tool to do what they ask for? That will help future readers of this question. $\endgroup$ – terdon Dec 4 '18 at 13:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.