I've got a new issue...
I would like to compute the average distance between two distinct motifs in the same sequences list than before. Have you some clue on how manage this ???
I just come to do it for one motif like that :
source("motifOccurrence.R")
motif <- c("T", "C", "A", "A")
motidist <- sapply(df, FUN=function(df, motif) {
computeDistance(coordMotif(df, motif))
}, motif = motif)
This R code give me the average distance between a define motif inside all the sequences of my list. And I would like to do the same but with two motif... Can someone help me ?
To extend informations about what i wanted to do :
i worked with a fasta file in input :
'> 1
GACTCTACTATAAACGGGAGATAGCAATCTAACGCAGTGCTTCAACTCCTCCTCCATCTGAACACCCTTCAACCTTTGATACTCAGACGTTTTAGGTCGG
'> 2
ACCACCCCTTTGTCCAGAAATAGGACTCTTGGGCCTGTTGCCTGAATAAAGTCCAACCACCACAACCACTACACTACCATATGTAAGCTTCACTGATGGT
'> 3
CACCACAAGTGCGCGCCACGACGTGCATAGCCTCTAGATCGGCAACTCAGGCGAGAAGTGTTTTATTTCGGTGTGGCCGGTCCTGGGCATTTTACGGAAA
'> 4
GTTAGTGTACAAGTCCGAATAGAGTCACGAAAGACCCACACAACCACGTAATGACCTCGCTGTAATGAGATCAGTTGGCTCATGAAGGAAGAACGTAATG
'> 5
TGAGCGTTCGCCAATAACCATCCCTCTCGTTCCTTGTAACTGTACTATGATAGCGGGCGCCCCCCTAATTAAATAGCGGACGCCCTGACCTATTGTATGA
'> 6
TGATATATCTACTCGATAAGGATATAGAGGTCTAATTGTTGAGAAGTGTACCACCTTAGAGCACGAGTTTAGGATACTTAGTAGGTTCTTGCGAAGGATA
etcetc
Then the motidist object look like this :
1 2
152 94
3 4
36 138
5 6
92 113
And the distances given by the function stand for one motif, and now, i would like to do the same but for the disatnce between two motifs like this :
atcgacatagacgactgatcgtcag MOTIF1 acggtagacagt MOTIF2 agcagatgacta # And this for all sequences in the file !
Thanks by advance
2 answers
Assuming I understood your problem this is my attempt to a solution. I do not use the functions in the script you pointed out, but the Biostrings Bioconductor package. Also, I use the Bioconductor package BSgenome.Scerevisiae.UCSC.sacCer3 for the yeast genome as example sequence data.
# load required packages.
library(ggplot2) # for qplot().
library(Biostrings)
library(BSgenome.Scerevisiae.UCSC.sacCer3)
# get sequence data (DNAStringSet)
seq <- getSeq(Scerevisiae)
# create motif dictionary.
dict0 <- DNAStringSet(
x = c("TCAA",
"GGAT",
"ATAT",
"GGCC")
)
pdict0 <- PDict(dict0)
# search with the dictionary every sequence in seq.
res <- lapply(seq, function(s) matchPDict(pdict0, s))
res[[1]]
MIndex object of length 4
[[1]]
IRanges object with 1533 ranges and 0 metadata columns:
start end width
<integer> <integer> <integer>
[1] 114 117 4
[2] 161 164 4
[3] 456 459 4
[4] 719 722 4
[5] 776 779 4
... ... ... ...
[1529] 228963 228966 4
[1530] 229173 229176 4
[1531] 229260 229263 4
[1532] 229817 229820 4
[1533] 229825 229828 4
...
<3 more elements>
# calculate combinations of motifs to compare
# assuming you want to compute different motifs to each other.
motifcomp <- t(combn(seq_len(length(pdict0)), 2))
colnames(motifcomp) <- c("motif_i", "motid_j")
motifcomp
motif_i motid_j
[1,] 1 2
[2,] 1 3
[3,] 1 4
[4,] 2 3
[5,] 2 4
[6,] 3 4
# iterate over the comparison matrix. for each pair of motifs and
# compute the distance to the nearest one over all sequences.
# this gives a list with each elements being the distances.
foo <- lapply(seq_len(nrow(motifcomp)), function(i) {
m <- motifcomp[i,]
d <- sapply(res, function(r) {
mcols(distanceToNearest(r[[m[1]]], r[[m[2]]]))$distance
})
unlist(d)
})
# we can check these distances are not normally distributed.
qplot(foo[[1]])
# now compute mean (and maybe other summaries, here standard deviation).
hoo <- lapply(foo, function(x) {
data.frame(mean = mean(x), sd = sd(x))
})
hoo <- do.call(rbind, hoo)
data.frame(motifcomp, hoo)
motif_i motid_j mean sd
1 1 2 150.15481 154.67573
2 1 3 71.63235 85.10428
3 1 4 358.76577 390.83777
4 2 3 69.26831 87.50602
5 2 4 348.63816 389.88737
6 3 4 378.64977 409.81555
I think that I already give answer here. Moreover, the distance between two motifs A and B is usually formulated as geometric distribution Where:
P(A,B,distance =k)=p(A)P(B)*(1-P(B))^k
You can simulate k (exemple 1:10000), and the Expected value will be the sum of pi*ni
Log in to answer this question.
Can you detailled your answer with a running example
I've edited my post with more details.
Do you need it in R or any other language would work?
R would be the easier way for me, but i can handle with python or perl if you have an idea in those languages !
What does
source("motifOccurrence.R")does? What doesdflooks like? A minimum working example would help.For source("motifOccurrence.R") ==> https://www.r-bloggers.com/calculate-the-average-distance-between-a-given-dna-motif-within-dna-sequences-in-r/
And df is just DNAStringSet instance looking like this :
Can't you just loop over the sequence, record the index of the motif if found and subtract that when you find the next one? Additional question: what should be done if you see:
1) atcgacatagacgactgatcgtcag MOTIF1 acggtagacagt MOTIF2 agcagatgacta MOTIF2 acgtcgtagctgatgctcggct (twice motif 2 after each other)
2) atcgacatagacgactgatcgtcag MOTIF1 acggtagacagt MOTIF2 agcagatgactaacgtgtgtgtg MOTIF1 acgtcgtagctgatgctcggct (motif2 sandwhiched by motif1 with different lengths)
3) atcgacatagacgactgatcgtcag MOTIF1 acggtagacagtagcagatgactaacgtcgtagctgatgctcggct (just a motif1 without motif2)
What would you mean by loop over the sequence ?
And for the additional questions :
1) I would like to know only the distance between the closer motifs motif 1 and first motif 2.
2) In this case the two informations are interesting and in a first time the average distance will be a sufficient info.
3) If there is no motif 2, return me 0.
You should have a look at the
repython module, regular expression, and just get the index of the position(s) in which a motif is found. This should be pretty easy.