5
$\begingroup$

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)
$\endgroup$
9
  • 1
    $\begingroup$ 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. $\endgroup$
    – terdon
    Sep 15, 2017 at 16:01
  • $\begingroup$ @terdon I've added my code. Is this a sufficient MWE? $\endgroup$
    – syntonicC
    Sep 15, 2017 at 16:20
  • $\begingroup$ 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? $\endgroup$
    – llrs
    Sep 16, 2017 at 11:36
  • $\begingroup$ @ 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. $\endgroup$
    – syntonicC
    Sep 16, 2017 at 19:19
  • 1
    $\begingroup$ @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. $\endgroup$
    – syntonicC
    Sep 20, 2017 at 15:04

2 Answers 2

5
$\begingroup$

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

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

$\endgroup$
1
  • 1
    $\begingroup$ 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 $\endgroup$
    – llrs
    Sep 26, 2017 at 20:28
1
$\begingroup$

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

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.

$\endgroup$
7
  • $\begingroup$ 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. $\endgroup$
    – syntonicC
    Sep 25, 2017 at 15:00
  • $\begingroup$ 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. $\endgroup$
    – syntonicC
    Sep 25, 2017 at 15:04
  • $\begingroup$ 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. $\endgroup$
    – syntonicC
    Sep 25, 2017 at 15:15
  • $\begingroup$ 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) $\endgroup$
    – llrs
    Sep 25, 2017 at 18:59
  • $\begingroup$ 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! $\endgroup$
    – syntonicC
    Sep 25, 2017 at 19:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.