2
$\begingroup$

I have an optimization problem where I have a degenerate nucleotide sequence I want to align against subsets of a reference genome (exons, specifically, to make the problem more tractable).

The sequence itself contains subregions of variable length.

The sequence looks like this, containing degenerate bases:

           L           C           R
5' - SSS......SSS|NNN.....NNN|WWW......WWW - 3'

Per IUPAC notation, S denotes C or G, N denotes any base, and W denotes A or T.

I am marking each block L, C, and R to broadly denote "left", "center", and "right" regions, oriented from 5' to 3' ("forward"-strand oriented).

The widths of the L and R blocks can be from 20 to 25 bases long. The width of the C block can be 14 to 16 bases long.

One programmatic way I could think to do this would be to brute-force a solution, generating all combinations of sequences of varied lengths, blastn-aligning each of such combinations to the sequence for the exon, and then looking for the best alignments among such sequences.

The output would be the scores for regions from sequences that align, to be able to say, from point X to point Y of the genome, a sequence here aligns better than other candidate sequences.

To avoid reinventing tools, and avoid brute-force methods, I'd like to know if there is a tool or even a wrapper around tools, which does this, in case it is already written, or another approach that is more elegant/efficient.

Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

I believe mafft can align degenerate sequences. From the documentation for version 7 at: https://mafft.cbrc.jp/alignment/software/anysymbol.html

Acceptable nucleotide symbols in the default mode:
a, c, g, t, u
Ambiguous nucleotides (r, y, w, s, k, m, d, v, h, b; IUPAC-IUB codes) can be used and are scored as:
score(r,a) = ( score(a,a) + score(g,a) ) / 2 
score(r,g) = ( score(a,g) + score(g,g) ) / 2 
score(r,t) = ( score(a,t) + score(g,t) ) / 2 
score(r,c) = ( score(a,c) + score(g,c) ) / 2 
... 
score(r,y) = ( score(r,t) + score(r,c) ) / 2 
... 
score(r,r) = ( score(a,a) + score(g,g) ) / 2 
score(y,y) = ( score(c,c) + score(t,t) ) / 2

MAFFT is designed for multiple sequence alignment, so will attempt to find a local alignment of homologous segments (using FFT), and then attempt to make a global alignment between the input sequences.

Improve this answer
$\endgroup$
4
  • 1
    $\begingroup$ Can you expand your answer to include some support material? $\endgroup$
    – Bioathlete
    Aug 11, 2019 at 14:30
  • $\begingroup$ I thought that a nucleotide blast could do the same. Just wondering if these tools also support querying against variable-width subsequences, or if wrappers around such tools can do this more efficiently in generalized cases. I have worked on parallelized analysis of combinations of substrings, in the meantime. It is brute-force but the scope of the problem seems tractable. $\endgroup$
    – Alex Reynolds
    Aug 12, 2019 at 1:06
  • $\begingroup$ I'm not sure I quite understand your point about variable-width subsequences - as far as I understand, these tools tend to find the longest alignable subsequences they can from the two input sequences, so I guess if you want to test lost of different widths a wrapper would be the way to go. $\endgroup$
    – Jay Moore
    Aug 12, 2019 at 8:52
  • $\begingroup$ @AlexReynolds I don't see why this isn't a 'normal' alignment problem where you could use blast or similar but with relaxed parameters (maybe permissive gap penalties to handle "variable-width" subsequences etc.). Perhaps a little more biological context would help us understand better too $\endgroup$
    – Chris_Rands
    Aug 12, 2019 at 9:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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