Tal; I'm trimming the BioC posting. In the R lists it is considered
spamming to cross post. (Please re-read the Posting Guide.)
On Dec 21, 2010, at 4:21 AM, Tal Galili wrote:
Hello everyone,
I am not sure if this should go on the general R mailing list (for
example,
if there is a text mining solution that might work here) or the
bioconductor
mailing list (since I wasn't able to find a solution to my question on
searching their lists) - so this time I tried both, and in the
future I'll
know better (in case it should go to only one of the two).
The task I'm trying to achieve is to align several sequences together.
I don't have a basic pattern to match to. All that I know is that the
"True" pattern should be of length "30" and that the sequences I'm
looking
at, have had missing values introduced to them at random points.
Here is an example of such sequences, were on the left we see what
is the
real location of the missing values, and on the right we see the
sequence
that we will be able to observe. My goal is to reconstruct the left
column
using only the sequences I've got on the right column (based on the
fact
that many of the letters in each position are the same)
Real_sequence The_sequence_we_see
1 CGCAATACTAAC-AGCTGACTTACGCACCG CGCAATACTAACAGCTGACTTACGCACCG
2 CGCAATACTAGC-AGGTGACTTCC-CT-CG CGCAATACTAGCAGGTGACTTCCCTCG
3 CGCAATGATCAC--GGTGGCTCCCGGTGCG CGCAATGATCACGGTGGCTCCCGGTGCG
4 CGCAATACTAACCA-CTAACT--CGCTGCG CGCAATACTAACCACTAACTCGCTGCG
5 CGCACGGGTAAGAACGTGA-TTACGCTCAG CGCACGGGTAAGAACGTGATTACGCTCAG
6 CGCTATACTAACAA-GTG-CTTAGGC-CTG CGCTATACTAACAAGTGCTTAGGCCTG
7 CCCA-C-CTAA-ACGGTGACTTACGCTCCG CCCACCTAAACGGTGACTTACGCTCCG
The agrep function allows one to specify which sort of differences to
consider in calculating a Levenshtein edit distance. Insertions are
one possible distance component. You could take a look at its code (in
C in hte sources) and perhaps rejigger it to spit out the location of
the deletions.
> agrep(seqdat$The_sequence_we_see[1], seqdat$Real_sequence,
max.distance=list(deletions=0, substitutions=0, insertions=0))
integer(0)
> agrep(seqdat$The_sequence_we_see[1], seqdat$Real_sequence,
max.distance=list(deletions=0, substitutions=0, insertions=1))
[1] 1
--
David.
Here is an example code to reproduce the above example:
ATCG <- c("A","T","C","G")
set.seed(40)
original.seq <- sample(ATCG, 30, T)
seqS <- matrix(original.seq,200,30, T)
change.letters <- function(x, number.of.changes = 15,
letters.to.change.with = ATCG)
{
number.of.changes <- sample(seq_len(number.of.changes), 1)
new.letters <- sample(letters.to.change.with , number.of.changes,
T)
where.to.change.the.letters <- sample(seq_along(x) ,
number.of.changes, F)
x[where.to.change.the.letters] <- new.letters
return(x)
}
change.letters(original.seq)
insert.missing.values <- function(x) change.letters(x, 3, "-")
insert.missing.values(original.seq)
seqS2 <- t(apply(seqS, 1, change.letters))
seqS3 <- t(apply(seqS2, 1, insert.missing.values))
seqS4 <- apply(seqS3,1, function(x) {paste(x, collapse = "")})
require(stringr)
# library(help=stringr)
all.seqS <- str_replace(seqS4,"-" , "")
# how do we allign this?
data.frame(Real_sequence = seqS4, The_sequence_we_see = all.seqS)
I understand that if all I had was a string and a pattern I would be
able to
use
library(Biostrings)
pairwiseAlignment(...)
But in the case I present we are dealing with many sequences to
align to one
another (instead of aligning them to one pattern).
Is there a known method for doing this in R?
Thanks,
Tal
----------------Contact
Details:-------------------------------------------------------
Contact me: tal.gal...@gmail.com | 972-52-7275845
Read me: www.talgalili.com (Hebrew) | www.biostatistics.co.il
(Hebrew) |
www.r-statistics.com (English)
----------------------------------------------------------------------------------------------
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
David Winsemius, MD
West Hartford, CT
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
