<p>I assume you have the sequence in seq<sup>1</sup>:</p>

<pre><code># Set the input
seq &lt;- getSequences("GRIN1")
motif &lt;- "ggcc"
# Convert to lower case
seq2 &lt;- tolower(seqinr::s2c(as.character(seq)))
istarts &lt;- seq(from = 1 , to = length(seq2), by = 1)
# Create words of the same size as the motif
oligos &lt;- seq2[istarts]
for (i in 2:nchar(motif)) {
        oligos &lt;- paste(seq2, seq2[istarts + i - 1], sep = "")
}
# Remove exceeding words (words with NA, might need tweaking to motif size
oligos &lt;- oligos[-c(length(oligos):(length(oligos)-nchar(motif)+2))]
# Find motif in the sequence
hits &lt;- grep(motif, oligos, ignore.case = TRUE)
#Prepare the plot
positions &lt;- ifelse(1:length(seq2) %in% hits, 1, 0)
plot(1:length(seq2), positions, col = c("white", "black")[position+1], pch = 19)
</code></pre>

<p>The plot is not the density, but you can choose how to plot it now, as you might be interested in windows sizes of 5 nucleotides or 20. 
<a href="https://i.sstatic.net/6CWmu.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/6CWmu.png" alt="positions where the motif is found in a sequence"></a></p>

<p>From here you can iterate through multiple sequences. If you want to find several motifs you can take advantage of the oligos object if they are of the same size.</p>

<p>There is the <a href="http://bioconductor.org/packages/seqPattern/" rel="nofollow noreferrer">seqPattern</a> package in Bioconductor for more advanced pattern matching. But surprisingly it doesn't have this simple motif search. I think it would be easier to do this with perl and pyhton than in R.</p>

<hr>

<p>1: In the code posted the function getSequence returns an error.</p>

CC BY-SA 3.0
2017-09-25T15:00:22.810Z
I fixed the error in `getSequences()`. You will need to change your first line of code to `seq <- getSequences("Grin1", "cdna")` and then also alter this line to `seq2 <- tolower(seqinr::s2c(as.character(seq$cdna)))`, assuming the seqType of interest is "cdna". I think the final line has a typo: the variable is `positions` not `position`.
syntonicC
CC BY-SA 3.0
2017-09-25T15:04:57.417Z
Thank you for your response! However, though I asked for this distribution in my question title, in the question body I stated that I was interested in determining the intron-exon boundaries as well.
syntonicC
CC BY-SA 3.0
2017-09-25T15:15:58.750Z
Now that I am reading more into this, numbering exons doesn't appear to be that simple. It looks like it is not standardized because there can be internal splice sites within exons as well as exons on different strands. This is probably why I am having so much trouble working this out.
syntonicC
CC BY-SA 3.0
2017-09-25T18:59:55.373Z
So, is your problem how to calculate the distribution of a motif in a sequence or find an intron-exon boundary? Which answer do you want. Please clarify (I spent quite a bit of time trying to answer this question properly but this format is just for one answer for one question)
llrs
CC BY-SA 3.0
2017-09-25T19:14:47.697Z
Well, it's really both. In the question body I said: "In this way, I can tell the distribution of various DNA sequences across the entire transcript visually and see differences among alternatively spliced isoforms." I should have put this in the question title but it's too late now; I should have been more specific. You answered the question in the title exactly as I asked it so I'm going to accept your answer. Thanks!
syntonicC
CC BY-SA 3.0
2017-09-25T20:27:31.407Z
Glad to help a bit! It was a nice problem to solve in R.
llrs
CC BY-SA 3.0
2017-09-28T14:53:53.327Z
@syntonicC In case you haven't seen, the other answer is better as it answers both of yours questions.
llrs
CC BY-SA 3.0
2017-09-25T07:56:54.053Z
2017-09-25T18:57:41.063Z
2017-09-25T18:57:41.063Z
llrs
<p>If I understand the question correctly, you'd like to plot the positions of the matches to you motif along with a gene model that shows the positions of introns and exons for the different transcripts. This can be accomplished fairly easily with ggbio:</p>
<pre><code>library(EnsDb.Mmusculus.v79)
library(ggbio)
library(biomaRt)
library(stringr)
library(dplyr)

# plot gene model with ggbio
gene_model &lt;- autoplot(EnsDb.Mmusculus.v79, ~ symbol == &quot;Grin1&quot;)

ensembl &lt;- useMart(&quot;ensembl&quot;, dataset = &quot;mmusculus_gene_ensembl&quot;)

# You need the full length gene sequence to match the gene model
seq &lt;- biomaRt::getSequence(id = &quot;Grin1&quot;,
                         type = &quot;mgi_symbol&quot;,
                         seqType = &quot;gene_exon_intron&quot;,
                         mart = ensembl)

# retrieve position of gene in genome to match up with gene model
gene_start=biomaRt::getBM(attributes=c(&quot;start_position&quot;), filters = 'external_gene_name', values = 'Grin1', mart=ensembl)

# use stringR to find positions of exact matches
motif &lt;- &quot;ggcc&quot;
motif_positions &lt;- as.data.frame(str_locate_all(seq$gene_exon_intron, toupper(motif))) %&gt;%
  mutate(start=start+gene_start[[1]]-1)

# make plot similar to Llopis's answer 
motif_plot &lt;- ggplot(motif_positions, aes(x=start)) +
    geom_point(y=1) +
    scale_y_continuous(limits=c(0,1))

# combine plots
tracks(motif= motifs, gene_model=gene_model)
</code></pre>
<p><a href="https://i.sstatic.net/JHJ8x.png" rel="noreferrer"><img src="https://i.sstatic.net/JHJ8x.png" alt="Gene Model with Motif matchs" /></a></p>
<p>This gives all motif matches for the gene sequence, not just the transcripts. It would be quite straightforward to filter out positions that overlap introns</p>

CC BY-SA 3.0
2017-09-26T20:28:15.957Z
I couldn't remember the name of the package ggbio (I was beginning to  doubt it existed). This is far more informative than my answer (+1) BTW, It might be better to overlap the plots of the gene model and the motif, with different colors
llrs
CC BY-SA 3.0
2017-09-26T11:49:15.910Z
2017-09-26T11:49:15.910Z
2020-06-18T08:30:41.743Z
heathobrien
<p>I would like to create a density/histogram of the distribution of a particular DNA sequence over the entire transcript using R and/or command line tools. From here, I would like to use the coordinates of the bins to map the intron-exon diagram below the plot. In this way, I can tell the distribution of various DNA sequences across the entire transcript visually and see differences among alternatively spliced isoforms. </p>

<p>I have already found the number of occurrences of the sequences I am interested in using the <code>seqinr</code> and <code>biostrings</code> packages. I have been using <code>biomaRt</code> to get the exon start and stop sequences as well as the full gene sequence. Between all this, I feel I have all the tools I need but I am missing how to put it all together. </p>

<p>The issue I am having is that <code>biomaRt</code> returns multiple exon sequences with identical start positions but different end positions for a total of 37 exons sequences. Meanwhile, <a href="http://genatlas.medecine.univ-paris5.fr/fiche.php?symbol=grin1" rel="nofollow noreferrer">Genatlas</a> and <a href="https://www.ncbi.nlm.nih.gov/IEB/Research/Acembly/av.cgi?db=mouse&amp;c=Gene&amp;l=Grin1" rel="nofollow noreferrer">NCBI</a> are returning anywhere from 19 to 21 exons for Grin1, my gene of interest. How can I put all of these tools together to get the plot I need? Is there a better approach?</p>

<p>Finally, if there is an R package or command line tool that already has this function built in, that would be ideal.</p>

<p>Example code for one transcript. This code is actually from an earlier project where I used multiple transcripts but it works fine with one as well. I masked the sequences I'm searching for as "XXXX" at the request of my project manager.</p>

<pre><code>library(magrittr)
library(biomaRt)
library(org.Mm.eg.db)
library(biostrings)
library(seqinr)

### Select only the longest sequences
getMaxLengthSequences &lt;- function(x) {
  seqMax &lt;- split(x, factor(x$mgi_symbol))
  seqMax &lt;- lapply(seqMax, function(x) max(x$length)) %&gt;% 
    unlist() %&gt;% 
    as.data.frame()
  seqMax$mgi_symbol &lt;- row.names(seqMax)
  row.names(seqMax) &lt;- NULL
  names(seqMax) &lt;- c("length", "mgi_symbol")
  seqUnique &lt;- merge(seqMax, x)
  seqUnique &lt;- seqUnique[order(seqUnique$mgi_symbol), ]

  return(seqUnique)
}

### Get sequence for gene of interest
getSequences &lt;- function(genelist, seqType) {
  ensembl &lt;- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

  sequ &lt;- biomaRt::getSequence(id = genelist,
                               type = "mgi_symbol",
                               seqType = seqType,
                               mart = ensembl)

  ### Remove unavailable sequences
  sequ &lt;- sequ[!grepl("Sequence unavailable", sequ[ ,seqType], fixed = TRUE), ]

  ### Find the length of each sequence then pick the max length sequence if there are multiple sequences returned from BioMart
  sequ &lt;- sequ[order(sequ$mgi_symbol), ]
  sequ$length &lt;- sapply(sequ[ ,seqType], nchar)
  sequ &lt;- getMaxLengthSequences(sequ)

  return(sequ)
}

### Get motif counts of interest (substrings are masked)
getmotifs &lt;- function(x) {

  seq &lt;- x %&gt;% 
    tolower() %&gt;% 
    s2c()

  output &lt;- c("XXXX" = count(seq, 4)[c("xxxx", "xxxx")] %&gt;% sum(),
              "XXXX" = count(seq, 4)[c("xxxx", "xxxx")] %&gt;% sum(),
              "XXX" = count(seq, 3)["xxx"] %&gt;% sum())

  return(output)
}

### Using the above functions: Get sequence and then find motifs
seq &lt;- getSequences("Grin1")
targetMotifs &lt;- matrix(sapply(seq[ ,"coding"], getmotifs), nrow = nrow(targets), ncol = 3, byrow = TRUE)

targetMotifTable &lt;- data.frame(Length = seq$length)
rownames(targetMotifTable) &lt;- make.names(seq$mgi_symbol, unique = TRUE)
targetMotifTable &lt;- cbind(targetMotifTable, targetMotifs) %&gt;% as.data.frame()
names(targetMotifTable) &lt;- c("Length", "XXXX", "XXXX", "XXX")

### Use biomaRt to get exons positions
ensembl &lt;- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

gb &lt;- getBM(attributes = c('ensembl_exon_id', "exon_chrom_start","exon_chrom_end","gene_exon"), 
            filters = "mgi_symbol", 
            values="Grin1", 
            mart = ensembl, 
            bmHeader = TRUE)
</code></pre>

CC BY-SA 3.0
2017-09-15T16:01:21.430Z
It would be much easier to help you if you could provide us with a minimal example that reproduces the issue that we can run at home and play with. Especially the specific transcript and the commands you use to get the different numbers of exons.
terdon
CC BY-SA 3.0
2017-09-15T16:20:24.773Z
@terdon I've added my code. Is this a sufficient MWE?
syntonicC
CC BY-SA 3.0
2017-09-16T11:36:33.390Z
If you have the string you search in X:Y can it be on x+1:y+1? you can look for Sushi, Gviz packages. Also, why don't you simply use grep with your motifs?
llrs
CC BY-SA 3.0
2017-09-16T19:19:41.387Z
@ Llopis Can you clarify what you mean by X and Y? I don't understand. Thank you for the package recs! This might be exactly what I need but it will take me a while to look through it since there are a lot of features available. As for grep... I'm not exactly sure, but I wrote this code last year and I think grep was too slow because I was going over 35000 genes with this code. From what I remember, the `count()` function had some C wizardry going on that made it faster. I think this is why I used it instead of grep but I can't remember exactly.
syntonicC
CC BY-SA 3.0
2017-09-18T06:44:41.583Z
@syntonicC For X and Y I mean positions in the genome. The reason I ask this is to know if you are interested in overlapping motifs or not.
llrs
CC BY-SA 3.0
2017-09-18T13:01:26.090Z
@syntonicC One more thing, could you please provide example motifs/sequence and sequences where they should be to play with it? Or point some ids as examples.
llrs
CC BY-SA 3.0
2017-09-20T14:15:25.147Z
@Llopis Sorry, I'm still new to this and can't figure out what you are asking. Does X and Y indicate start and stop positions along a chromosome or something else? Overlapping motifs: If I search for a substring wouldn't this capture any motifs that overlap? I am interested in overlapping motifs. Some motifs I am interested in: "ggcc", "gac", "agga", "gaga". Is this what you are looking for?
syntonicC
CC BY-SA 3.0
2017-09-20T14:49:53.240Z
@syntonicC Yes, X is the beginning and Y is the end of the region (if this is start or stop of an open reading frame I don't know). Yes, I am interested in this kind of thing, for instance in your data the 3rd and 4th motifs could be overlap in **ag*ga**ga*ta (here in bold and in italics), and some example of gene/transcript (a good model or I will pick one randommly)
llrs
CC BY-SA 3.0
2017-09-20T15:04:36.897Z
@Llopis Ah sorry, I mentioned the gene in my code (and URLs) but did not mention in the question body that it was mouse GRIN1. The coordinates are: 25174720 to 25146695 on the reverse strand. Based on what you said, yes, I do want overlapping motifs as well.
syntonicC
<r><software-recommendation><sequence-annotation>
Software recommendation: find DNA sequence distribution over entire transcript
# Software recommendation: find DNA sequence distribution over entire transcript

(Asked by: username_0 on 2017-09-15)

I would like to create a density/histogram of the distribution of a particular DNA sequence over the entire transcript using R and/or command line tools. From here, I would like to use the coordinates of the bins to map the intron-exon diagram below the plot. In this way, I can tell the distribution of various DNA sequences across the entire transcript visually and see differences among alternatively spliced isoforms. 
I have already found the number of occurrences of the sequences I am interested in using the seqinr and biostrings packages. I have been using biomaRt to get the exon start and stop sequences as well as the full gene sequence. Between all this, I feel I have all the tools I need but I am missing how to put it all together. 
The issue I am having is that biomaRt returns multiple exon sequences with identical start positions but different end positions for a total of 37 exons sequences. Meanwhile, Genatlas and NCBI are returning anywhere from 19 to 21 exons for Grin1, my gene of interest. How can I put all of these tools together to get the plot I need? Is there a better approach?
Finally, if there is an R package or command line tool that already has this function built in, that would be ideal.
Example code for one transcript. This code is actually from an earlier project where I used multiple transcripts but it works fine with one as well. I masked the sequences I'm searching for as "XXXX" at the request of my project manager.
library(magrittr)
library(biomaRt)
library(org.Mm.eg.db)
library(biostrings)
library(seqinr)

### Select only the longest sequences
getMaxLengthSequences <- function(x) {
  seqMax <- split(x, factor(x$mgi_symbol))
  seqMax <- lapply(seqMax, function(x) max(x$length)) %>% 
    unlist() %>% 
    as.data.frame()
  seqMax$mgi_symbol <- row.names(seqMax)
  row.names(seqMax) <- NULL
  names(seqMax) <- c("length", "mgi_symbol")
  seqUnique <- merge(seqMax, x)
  seqUnique <- seqUnique[order(seqUnique$mgi_symbol), ]

  return(seqUnique)
}

### Get sequence for gene of interest
getSequences <- function(genelist, seqType) {
  ensembl <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

  sequ <- biomaRt::getSequence(id = genelist,
                               type = "mgi_symbol",
                               seqType = seqType,
                               mart = ensembl)

  ### Remove unavailable sequences
  sequ <- sequ[!grepl("Sequence unavailable", sequ[ ,seqType], fixed = TRUE), ]

  ### Find the length of each sequence then pick the max length sequence if there are multiple sequences returned from BioMart
  sequ <- sequ[order(sequ$mgi_symbol), ]
  sequ$length <- sapply(sequ[ ,seqType], nchar)
  sequ <- getMaxLengthSequences(sequ)

  return(sequ)
}

### Get motif counts of interest (substrings are masked)
getmotifs <- function(x) {

  seq <- x %>% 
    tolower() %>% 
    s2c()

  output <- c("XXXX" = count(seq, 4)[c("xxxx", "xxxx")] %>% sum(),
              "XXXX" = count(seq, 4)[c("xxxx", "xxxx")] %>% sum(),
              "XXX" = count(seq, 3)["xxx"] %>% sum())

  return(output)
}

### Using the above functions: Get sequence and then find motifs
seq <- getSequences("Grin1")
targetMotifs <- matrix(sapply(seq[ ,"coding"], getmotifs), nrow = nrow(targets), ncol = 3, byrow = TRUE)

targetMotifTable <- data.frame(Length = seq$length)
rownames(targetMotifTable) <- make.names(seq$mgi_symbol, unique = TRUE)
targetMotifTable <- cbind(targetMotifTable, targetMotifs) %>% as.data.frame()
names(targetMotifTable) <- c("Length", "XXXX", "XXXX", "XXX")

### Use biomaRt to get exons positions
ensembl <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

gb <- getBM(attributes = c('ensembl_exon_id', "exon_chrom_start","exon_chrom_end","gene_exon"), 
            filters = "mgi_symbol", 
            values="Grin1", 
            mart = ensembl, 
            bmHeader = TRUE)

**Tags:** r, software-recommendation, sequence-annotation
> **Comments:**
> **username_1**: It would be much easier to help you if you could provide us with a minimal example that reproduces the issue that we can run at home and play with. Especially the specific transcript and the commands you use to get the different numbers of exons.  
> **username_0**: @terdon I've added my code. Is this a sufficient MWE?  
> **username_2**: If you have the string you search in X:Y can it be on x+1:y+1? you can look for Sushi, Gviz packages. Also, why don't you simply use grep with your motifs?  
> **username_0**: @ Llopis Can you clarify what you mean by X and Y? I don't understand. Thank you for the package recs! This might be exactly what I need but it will take me a while to look through it since there are a lot of features available. As for grep... I'm not exactly sure, but I wrote this code last year and I think grep was too slow because I was going over 35000 genes with this code. From what I remember, the `count()` function had some C wizardry going on that made it faster. I think this is why I used it instead of grep but I can't remember exactly.  
> **username_2**: @syntonicC For X and Y I mean positions in the genome. The reason I ask this is to know if you are interested in overlapping motifs or not.  
> **username_2**: @syntonicC One more thing, could you please provide example motifs/sequence and sequences where they should be to play with it? Or point some ids as examples.  
> **username_0**: @Llopis Sorry, I'm still new to this and can't figure out what you are asking. Does X and Y indicate start and stop positions along a chromosome or something else? Overlapping motifs: If I search for a substring wouldn't this capture any motifs that overlap? I am interested in overlapping motifs. Some motifs I am interested in: "ggcc", "gac", "agga", "gaga". Is this what you are looking for?  
> **username_2**: @syntonicC Yes, X is the beginning and Y is the end of the region (if this is start or stop of an open reading frame I don't know). Yes, I am interested in this kind of thing, for instance in your data the 3rd and 4th motifs could be overlap in **ag*ga**ga*ta (here in bold and in italics), and some example of gene/transcript (a good model or I will pick one randommly)  
> **username_0**: @Llopis Ah sorry, I mentioned the gene in my code (and URLs) but did not mention in the question body that it was mouse GRIN1. The coordinates are: 25174720 to 25146695 on the reverse strand. Based on what you said, yes, I do want overlapping motifs as well.  

---

## 2 answers

### Answer by username_3 on 2017-09-26. Score: 5 | Accepted answer

If I understand the question correctly, you'd like to plot the positions of the matches to you motif along with a gene model that shows the positions of introns and exons for the different transcripts. This can be accomplished fairly easily with ggbio:
library(EnsDb.Mmusculus.v79)
library(ggbio)
library(biomaRt)
library(stringr)
library(dplyr)

# plot gene model with ggbio
gene_model <- autoplot(EnsDb.Mmusculus.v79, ~ symbol == "Grin1")

ensembl <- useMart("ensembl", dataset = "mmusculus_gene_ensembl")

# You need the full length gene sequence to match the gene model
seq <- biomaRt::getSequence(id = "Grin1",
                         type = "mgi_symbol",
                         seqType = "gene_exon_intron",
                         mart = ensembl)

# retrieve position of gene in genome to match up with gene model
gene_start=biomaRt::getBM(attributes=c("start_position"), filters = 'external_gene_name', values = 'Grin1', mart=ensembl)

# use stringR to find positions of exact matches
motif <- "ggcc"
motif_positions <- as.data.frame(str_locate_all(seq$gene_exon_intron, toupper(motif))) %>%
  mutate(start=start+gene_start[[1]]-1)

# make plot similar to Llopis's answer 
motif_plot <- ggplot(motif_positions, aes(x=start)) +
    geom_point(y=1) +
    scale_y_continuous(limits=c(0,1))

# combine plots
tracks(motif= motifs, gene_model=gene_model)


This gives all motif matches for the gene sequence, not just the transcripts. It would be quite straightforward to filter out positions that overlap introns

> **Comments:**
> **username_2**: I couldn't remember the name of the package ggbio (I was beginning to  doubt it existed). This is far more informative than my answer (+1) BTW, It might be better to overlap the plots of the gene model and the motif, with different colors  

---

### Answer by username_2 on 2017-09-25. Score: 1

I assume you have the sequence in seq1:
# Set the input
seq <- getSequences("GRIN1")
motif <- "ggcc"
# Convert to lower case
seq2 <- tolower(seqinr::s2c(as.character(seq)))
istarts <- seq(from = 1 , to = length(seq2), by = 1)
# Create words of the same size as the motif
oligos <- seq2[istarts]
for (i in 2:nchar(motif)) {
        oligos <- paste(seq2, seq2[istarts + i - 1], sep = "")
}
# Remove exceeding words (words with NA, might need tweaking to motif size
oligos <- oligos[-c(length(oligos):(length(oligos)-nchar(motif)+2))]
# Find motif in the sequence
hits <- grep(motif, oligos, ignore.case = TRUE)
#Prepare the plot
positions <- ifelse(1:length(seq2) %in% hits, 1, 0)
plot(1:length(seq2), positions, col = c("white", "black")[position+1], pch = 19)

The plot is not the density, but you can choose how to plot it now, as you might be interested in windows sizes of 5 nucleotides or 20. 

From here you can iterate through multiple sequences. If you want to find several motifs you can take advantage of the oligos object if they are of the same size.
There is the seqPattern package in Bioconductor for more advanced pattern matching. But surprisingly it doesn't have this simple motif search. I think it would be easier to do this with perl and pyhton than in R.

1: In the code posted the function getSequence returns an error.

> **Comments:**
> **username_0**: I fixed the error in `getSequences()`. You will need to change your first line of code to `seq <- getSequences("Grin1", "cdna")` and then also alter this line to `seq2 <- tolower(seqinr::s2c(as.character(seq$cdna)))`, assuming the seqType of interest is "cdna". I think the final line has a typo: the variable is `positions` not `position`.  
> **username_0**: Thank you for your response! However, though I asked for this distribution in my question title, in the question body I stated that I was interested in determining the intron-exon boundaries as well.  
> **username_0**: Now that I am reading more into this, numbering exons doesn't appear to be that simple. It looks like it is not standardized because there can be internal splice sites within exons as well as exons on different strands. This is probably why I am having so much trouble working this out.  
> **username_2**: So, is your problem how to calculate the distribution of a motif in a sequence or find an intron-exon boundary? Which answer do you want. Please clarify (I spent quite a bit of time trying to answer this question properly but this format is just for one answer for one question)  
> **username_0**: Well, it's really both. In the question body I said: "In this way, I can tell the distribution of various DNA sequences across the entire transcript visually and see differences among alternatively spliced isoforms." I should have put this in the question title but it's too late now; I should have been more specific. You answered the question in the title exactly as I asked it so I'm going to accept your answer. Thanks!  
> **username_2**: Glad to help a bit! It was a nice problem to solve in R.  
> **username_2**: @syntonicC In case you haven't seen, the other answer is better as it answers both of yours questions.  

---

