Htykuut

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.

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





#---------------------------------------

#---------------------------------------

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.

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 '>%sn%sn' $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

                    *    ***      **   **     *       *    * *    

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).