Aligning multiple sequences in heterogeneous group

Question

I have a list of ~200 DNA sequences, representing probably 50 different genomic regions but they're all mixed up. For example, if I have seq1, seq2.... seq10, seq1 might align to seq3 and seq8, but be completely unrelated to the others.

There are also some differences in the lengths sampled, so the above example might represent:

Seq1-------------------------------------------------

Seq3---------------------- seq8-----------------

Such that seq3 and seq8 don't align to each other, but both align to seq1

So what I'd like to do is go through and somehow generate a list of groups of sequences that align to each other, as well as the alignments. eg:

Group 1

Seq1-------------------------------------------------

Seq3---------------------- seq8-----------------

Group 2

seq2------------------------------------

. seq6-----------------------------

. seq7--------------------------xxxxxx

Group 3... etc

Trying ClustalW or MUSSLE to align everything doesn't work (or takes an unreasonable amount of time), I'm guessing because there are so many sequences that don't align at all. I tried making a custom BLAST database and then BLASTed each sequence against it, but then I get multiple hits for the same alignment (with the group 2 example above, seq2:seq6, seq2:seq7, seq6:seq2, seq6:seq7, seq7:seq2 and seq7:seq6 all get returned as 6 unique hits, when they should be grouped together.

My current coding knowledge is fairly basic, but I'm willing to read docs and figure stuff out, I just don't want to reinvent the wheel.

Edit2: Really, the grouping is the important part - once I have the groups, I can do the alignment separately with little effort. I'd just like to have groups where each sequence is in a single group.

Answer 1

You can try this:

• BLAST each sequence to every other sequence (pairwise).
• Every alignment (with some defined cutoff) denotes a connection.
• Map all connections.
• If a sequence is connected to some other directly or indirectly, it falls in a group. Put all the sequences that seq1 aligns to, in Group-1, then go to the alignments of these sequences; put all the sequences that these align to again in Group-1; so keep populating the group like this.

Methodology:

• Install standalone blast (if you don't have many sequences then you can run online BLAST as well)
• Make a blast database from your sequences using makeblastdb
• Align these sequences against the database. If you are using online BLAST use BL2seq (align two sequences). It is much better and convenient to use standalone. You can also mention if you want plus-plus or plus-minus or both alignments. You might in some cases want only either of the two.
• In the standalone BLAST you can specify the output format (which fields to include etc- the format you choose just depends on your requirements).

A tabular output format looks something like this:

# BLASTN 2.2.27+
# Query: TCONS_00036712 gene=XLOC_017996
# Database: ../nt_db/nt
# Fields: query id, subject id, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 1014 hits found
TCONS_00036712  gi|191174875|emb|CU655970.6|    95.54   202 9   0   423 624 16680   16479   8e-85    324
TCONS_00036712  gi|51491599|gb|AC144709.2|  95.02   201 10  0   424 624 28443   28243   1e-82    316

Ignore the commented (#) lines; the first field is query ID, the second is the subject ID and there is an alignment between the two; the other fields provide information about the alignment (You can choose these fields).

For parsing I use a fast and easy scripting language called awk which is included in all UNIX based systems. It is also available for windows in GNUWin32 package.

What you have to do is check the first two fields and update the group.

# MakeGroups.awk

BEGIN{FS="\t"}     # Declaring Field separator as Tab

!($1 in grp){          # Check if the seq is a parent group. If not...
    k=1
    for(i in grp){ 
        if($1 in grp[i]){    # Check if the seq is part of any other groups
            parentgrp[$1]=i
            if(!($2 in grp[i]))   # Check if the second field i.e subject already present in the parentgroup
                grp[i][$2]        # if not assign the second field to the parent group
            k=0
            break    # stop checking further
        }
    }
    if(k==1)        # There is no parent group with this label and the seq is not a part of any other group
        grp[$1][$1] # Make a group with the query-id as the label and add this query to this group.
}  

$1 in grp{
    if(!($2 in grp[$1]))
        grp[$1][$2]    # If the second field is not a part of the group with first field as label then assign it to that
}

END{
    for(i in grp){
        x++
        print "Group-"x"\n----------"
        for(j in grp[i])
            print j
    }
    print "\n"
}

Run this script like this in the terminal:
gawk -f MakeGroups.awk blastalignmentfile.txt

Note: This script contains multidimensional arrays. It will not work with all versions of awk. Use gawk version >4.0.

As swarnbes mentions in their answer, there are faster algorithms which do this kind of stuff and are used for sequence assembly. What many of them do is to make a graph (networks called deBruijn graphs), where each connection is an alignment, and compute an Eulerian Path. See this review by Pavel Pevzner for details. Overlapping sequences form contigs and you can easily trace back which sequence came from which contig (which you can call a Group). Each contig/Group is a disjoint deBruijn sub-graph.

Answer 2

Do you really need BLAST? That is, are the sequences different enough from each other that you need an algorithm that looks for large differences between them?

Maybe you could use something like Phrap, which should put together contigs for you, if the sequences which should go together are very close to identical.