Dear Hervé,

Thank you very much for your response.

Your script has a line of code that reads:

 if (gsub("\\[.\\]", "", seq) != ref.sequence)
  stop("'sequences[", i, "]' is not compatible with 'ref.sequence'")

I take it this means that your functions do not allow for polymorphisms and other variations. For instance:

ref.sequence <- "GCAAATT"
sequences <- read.table(textConnection("name sequence
sequence1 GC-AATT
sequence2 GCAAAT[T]T
sequence3 GCAAAGT"), header = TRUE, as.is<http://as.is> = TRUE)
closeAllConnections()

Sequence1 has a deletion (relative to the reference sequence) at position 3, sequence 2 has an insertion (relative to the reference sequence) between positions 5 and 6, and sequence3 has a G instead of a T (relative to the reference sequence) at position 6. (I am numbering the positions from left to right, using the reference sequence to assign the position numbers.)

What I want to get is this:

reference.sequence: "GCAAAT-T
sequence1: GC-AAT-T
sequence2: GCAAATTT
sequence3: GCAAAG-T

When I apply your functions to this new data, I get the following errors:

> stretchRefSequence(ref.sequence, as.character(sequences[,2]))
Error in computeInsSizesInRef(ref.sequence, sequences) :
  'sequences[1]' is not compatible with 'ref.sequence'

> stretchSequences(ref.sequence, as.character(sequences[,2]))
Error in computeInsSizesInRef(ref.sequence, sequences) :
  'sequences[1]' is not compatible with 'ref.sequence'

I tried deleting the following lines of code from your script:

 if (gsub("\\[.\\]", "", seq) != ref.sequence)
  stop("'sequences[", i, "]' is not compatible with 'ref.sequence'")

Reapplying the functions, I get:

> stretchRefSequence(ref.sequence, as.character(sequences[,2]))
[1] "GC-AAA-TT"
> #not correct

> stretchSequences(ref.sequence, as.character(sequences[,2]))
[1] "GC-AAT-T"  "GC-AAATTT" "GC-AAA-GT"
Warning messages:
1: In .Method(..., na.rm = na.rm) :
  an argument will be fractionally recycled
2: In ins_sizes_in_ref - ins_sizes :
  longer object length is not a multiple of shorter object length
3: In seq_along(ins_sizes_in_ref) + cumsum(ins_sizes) :
  longer object length is not a multiple of shorter object length
> #also not correct

Is there a way to modify your codes to make my new toy data give me the answers I'm looking for? Again, I really appreciate your help, because this parsing syntax is way over my head at the moment.

-----------------------------------
Josh Banta, Ph.D
Assistant Professor
Department of Biology
The University of Texas at Tyler
Tyler, TX 75799
Tel: (903) 565-5655
http://plantevolutionaryecology.org


	[[alternative HTML version deleted]]

