How to determine the average distance between known motif in a list of DNA sequences
1
1
Entering edit mode
8.2 years ago
kevinm ▴ 40

Hi everybody !

So there is my problem : I am searching for the average distance between a known motif inside sequence, and extend this to a list of sequences... The first part is done, the second part (extend to a list of sequences) is the problematic one ! So, here the way i am doing the first part :

source("motifOccurrence.R") #https://www.r-bloggers.com/calculate-the-average-distance-between-a-given-dna-motif-within-dna-sequences-in-r/

library("seqinr")

df <- readDNAStringSet("X.fasta")

df2 <- df[[1]]

motif <- c("T", "C", "C", "A")

coord <- coordMotif(df2, motif)

motidist <- computeDistance(coord)

motidist

[1] 152

It's appear that the first sequence of my fasta list have an average distance of 152 nucelotide between two TCCA motif. And, i don't know how automatize this to all my list in df...

Thanks by advance for the help.

Kévin

R sequence Motif Distance • 2.8k views
ADD COMMENT
0
Entering edit mode

Look at apply() and its variants.

ADD REPLY
0
Entering edit mode

Thanks, i will check it !

ADD REPLY
2
Entering edit mode
8.1 years ago
Macherki M E ▴ 120

The Rho function in the seqinr package calculate the indipendence of occurence of a k-words using the expected word count. We can modify it to be run just for our target:

library("seqinr")
data(ec999)

exp_count<-function (sequence, wordsize = 2, alphabet = s2c("acgt")) 
{
    wordcount <- count(sequence, wordsize, freq = FALSE, alphabet = alphabet)
    uni <- count(sequence, 1, freq = TRUE, alphabet = alphabet)
    expected_wordfreq <- function(wordsize, uni) {
        if (wordsize == 1) 
            return(uni)
        else kronecker(uni, expected_wordfreq(wordsize - 1, uni))
    }
    expected_wordcount <- sum(wordcount) * expected_wordfreq(wordsize, 
        uni)

names(expected_wordcount)<-words(wordsize,alphabet = alphabet)
return(expected_wordcount)
}
#exp_count(ec999[[1]],3)
ADD COMMENT

Login before adding your answer.

Traffic: 2683 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6