BUG: seq2aln() or seqaln()

[Xinqiu Yao]

It seems example codes for 'seq2aln()' are broken:

example(seq2aln, run.dontrun=TRUE)
seq2ln> aa.1 <- pdbseq( read.pdb("1bg2") )
  Note: Accessing on-line PDB file
  HEADER    MOTOR PROTEIN                           04-JUN-98   1BG2               

seq2ln> aa.2 <- pdbseq( read.pdb("3dc4") )
  Note: Accessing on-line PDB file
  HEADER    MOTOR PROTEIN                           03-JUN-08   3DC4               

seq2ln> aa.3 <- pdbseq( read.pdb("1mkj") )
  Note: Accessing on-line PDB file
  HEADER    TRANSPORT PROTEIN                       29-AUG-02   1MKJ               

seq2ln> aln <- seqaln( seqbind(aa.1,aa.2) )

seq2ln> seq2aln(aa.3, aln)
      1        .         .         .         .         .         .         70 
seq   ADLAECNIKVMCRFRPLNESEVNRGDKYIAKFQGEDTVVIASKPYAFDRVFQSSTSQEQVYNDCAKKIVK
      1        .         .         .         .         .         .         70 

     71        .         .         .         .         .         .         140 
seq   DVLEGYNGTIFAYGQTSSGKTHTMEGKLHDPEGMGIIPRIVQDIFNYIYSMDENLEFHIKVSYFEIYLDK
     71        .         .         .         .         .         .         140 

    141        .         .         .         .         .         .         210 
seq   IRDLLDVSKTNLSVHEDKNRVPYVKGCTERFVCSPDEVMDTIDEGKSNRHVAVTNMNEHSSRSHSIFLIN
    141        .         .         .         .         .         .         210 

    211        .         .         .         .         .         .         280 
seq   VKQENTQTEQKLSGKLYLVDLAGSEKVAKNINKSLSALGNVISALAEGSTYVPYRDSKMTRILQDSLGGN
    211        .         .         .         .         .         .         280 

    281        .         .         .         .         .     336 
seq   CRTTIVICCSPSSYNESETKSTLLFGQRAKTIKNTVCVNVELTAEQWKKKYEKEKE
    281        .         .         .         .         .     336 

Call:
  seq2aln(seq2add = aa.3, aln = aln)

Class:
  fasta

Alignment dimensions:
  1 sequence rows; 336 position columns (336 non-gap, 0 gap) 

+ attr: id, ali, call
Warning message:
In (function (aln, id = NULL, profile = NULL, exefile = "muscle",  :
  nothing to align

However, on releases branch, the example run is okay. We haven't changed 'seq2aln()' for a long time. So, I guess the bug might come from recent changes in 'seqaln()'...

[Lars Skjærven]

Hmm.. I added the following to seqaln:

+  ## nothing to align?
+  if(!nrow(aln$ali)>1) {
+    warning("nothing to align")
+    aln$ali <- aln$ali[ , !is.gap(aln$ali), drop=FALSE]
+    return(aln)
+  }

this is what produces the warning, but I didnt do anything in seq2aln...

[Xinqiu Yao]

This will affect the profile argument in seqaln(), which allows a single sequence aligned to an existing aln and is used by seq2aln(). I guess just add one more check in if for is.null(profile) should solve the problem. Let's see..

[Xinqiu Yao]

See the commit. Let me know if I was wrong..

[Lars Skjærven]

I think this should be fine now. I'm adding a test file for seqaln() as well. Since my "bugfix" leading to this bug was also introduced to releases I've merged your commit to the releases branch.