Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-09-25T15:00:22.810Z
Answers.Text: 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`.
Answers.UserDisplayName: syntonicC
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-09-25T15:04:57.417Z
Answers.Text: 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.
Answers.UserDisplayName: syntonicC
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-09-25T15:15:58.750Z
Answers.Text: 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.
Answers.UserDisplayName: syntonicC
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-09-25T18:59:55.373Z
Answers.Text: 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)
Answers.UserDisplayName: llrs
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-09-25T19:14:47.697Z
Answers.Text: 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!
Answers.UserDisplayName: syntonicC
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-09-25T20:27:31.407Z
Answers.Text: Glad to help a bit! It was a nice problem to solve in R.
Answers.UserDisplayName: llrs
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-09-28T14:53:53.327Z
Answers.Text: @syntonicC In case you haven't seen, the other answer is better as it answers both of yours questions.
Answers.UserDisplayName: llrs
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-09-25T07:56:54.053Z
Answers.LastActivityDate: 2017-09-25T18:57:41.063Z
Answers.LastEditDate: 2017-09-25T18:57:41.063Z
Answers.UserDisplayName: llrs
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-09-26T20:28:15.957Z
Answers.Text: 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
Answers.UserDisplayName: llrs
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-09-26T11:49:15.910Z
Answers.LastActivityDate: 2017-09-26T11:49:15.910Z
Answers.LastEditDate: 2020-06-18T08:30:41.743Z
Answers.UserDisplayName: heathobrien
Body: <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>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-09-15T16:01:21.430Z
Comments.Text: 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.
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-09-15T16:20:24.773Z
Comments.Text: @terdon I've added my code. Is this a sufficient MWE?
Comments.UserDisplayName: syntonicC
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-09-16T11:36:33.390Z
Comments.Text: 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?
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-09-16T19:19:41.387Z
Comments.Text: @ 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.
Comments.UserDisplayName: syntonicC
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-09-18T06:44:41.583Z
Comments.Text: @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.
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-09-18T13:01:26.090Z
Comments.Text: @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.
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-09-20T14:15:25.147Z
Comments.Text: @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?
Comments.UserDisplayName: syntonicC
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-09-20T14:49:53.240Z
Comments.Text: @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)
Comments.UserDisplayName: llrs
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2017-09-20T15:04:36.897Z
Comments.Text: @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.
Comments.UserDisplayName: syntonicC
Tags: <r><software-recommendation><sequence-annotation>
Title: Software recommendation: find DNA sequence distribution over entire transcript
ThreadText: # 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](http://genatlas.medecine.univ-paris5.fr/fiche.php?symbol=grin1) and [NCBI](https://www.ncbi.nlm.nih.gov/IEB/Research/Acembly/av.cgi?db=mouse&c=Gene&l=Grin1) 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**: @username_1 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**: @username_0 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**: @username_0 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**: @username_0 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)
  
  ```
  
  [![Gene Model with Motif matchs](https://i.sstatic.net/JHJ8x.png)](https://i.sstatic.net/JHJ8x.png)
  
  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.
  [![positions where the motif is found in a sequence](https://i.sstatic.net/6CWmu.png)](https://i.sstatic.net/6CWmu.png)
  
  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](http://bioconductor.org/packages/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**: @username_0 In case you haven't seen, the other answer is better as it answers both of yours questions.  
  
  ---
  
  