I am looking to show how a primer is consistent among some genomic data. I have a primer of about 23bp and looking to compare it to about 5000 genomic sequences of 10kb. I am unsure how to format it the way I need to, since I am just trying to show the nucleotide differences of the primer to the genetic sequences. Here is what I need:

>     1. Cut out the area that my primer is located and about 20bp down each end. 
>     2. Show only the bases that are different from my primer in my analysis.
>      ex: Primer: -----------ATGTGGAAGCAAATATCAAATGA---------
>          Gene:   ATGACCATACG----C--------------T---ATCGTAGGG
>                  ATGAGCATACC-----A----T--------T---TTCGAACGC

The data I am using is all dengue sequences (all serotypes) and the primer with the following code: ATGTGGAAGCAAATATCAAATGA.

I am using Biostrings, msa, and seqinr. I use ncbi to get the genetic sequences and using FASTA files.


Edit: Solution: I used MAFFT to align the data, working very efficiently and put it into a fasta file. Then I cut the nucleotides I did not need by: 1. Finding the part i wanted and cutting 10 nucleotides before it and after it. 2. Putting only the [start:end] into a new list, which I will then manipulate.

However, I am still having difficulty finding a way of showing my data nicely because I can't find a suitable package. Any advice?

Edit: Hi all, Just wanted to show the code I made to solve the problem.

    #So the we will be using a dataframe that has the sequences as a string
    #Then I created a temporary dataframe to store each sequence, with the primer at the top

    tempDF <- data.frame(matrix(primer1, nrow = 1, ncol = 1), 
    stringsAsFactors = FALSE)

    for(i in sequenceAsString){
      #Temporary vector for each sequence after compared to primer. 

      temporary1 <- ""
      temporary1 <- paste(temporary1, substring(i,1,11), sep = "") 

      for(j in 12:(nchar(i)-20)){
        if(substring(i,j,j) == substring(toString(reversePrimer1), j-12,j-12)){
          temporary1 <- paste(temporary1, "-", sep = "")
          temporary1 <- paste(temporary1, substring(i,j,j), sep = "")
      temporary1 <- paste(temporary1, substring(i,nchar(i)-19,nchar(i)), sep = "") 
      tempDF <- rbind(tempDF,temporary1)

So, I essentially take the sequence I chose apart and compare each character to that of the reverseprimer. I then replace the ones that are the same with an '-'

Thanks everyone for their help!

  • 2
    $\begingroup$ “Since that is too much for my computer to do” — No. Performing a short pairwise alignment like this is trivial. Running it against 5000 sequences of length 10knt should take a few seconds at most. Maybe a bit longer in R using pairwiseAlignment but totally feasible. And if performance is really of the essence, you can implement e.g. Myers Bit Vector since the primer is so short, and the whole search could happen in sub-second time. Either way, don’t perform a multiple sequence alignment. $\endgroup$ Sep 11 '18 at 10:10
  • 1
    $\begingroup$ The issue you are having is not entirely clear. When you say "However, it was difficult because to accurately predict if you would need to have it aligned.": Does this mean you do not know how to select the homologous region across all genes that encompass the primer? You may want to be more specific, ie provide detailed example, what you tried and exactly where issue is. $\endgroup$
    – Vince
    Sep 11 '18 at 14:21
  • $\begingroup$ Thanks a lot Konrad and Vince, very helpful. I am clearly just learning the ropes! I will get back with the result. $\endgroup$
    – Colin
    Sep 12 '18 at 16:50
  • $\begingroup$ So, I have figured it out using a combination of MAFFT and R. I still have not found a way to present the data in the format I mentioned, might you have any examples of packages that can be used? $\endgroup$
    – Colin
    Sep 16 '18 at 0:15

If I understand you correctly, after your most recent edit, you have 5,000 sequences that are 63 bp long and in each you simply want to replace the nucleotides from 21 through 43 with - unless that nucleotide does not match the primer. Since - is normally used to represent gaps, I recommend using . instead for matching nucleotides.

I know you're looking for an R solution, but I don't usually use R for string manipulations and don't have time at the moment to figure it out. So, assuming this might be urgent for you, below is a Python solution. This uses your example data which shows 11 nucleotides before the primer and 9 after it (so you'll need to adjust it to work with your actual case in which 20 nucleotides are before and after the primer). In Python, strings are arrays of characters, by the way.

Here are your example sequences are in the file seqs.txt:


Here is the code that gives you the expected output:

#!/bin/env python
import sys

print('-' * 11 + primer + '-' * 9 )
with open('seqs.txt','r') as f:
    for gene in f:
        for gene_index in range(11,34):
            primer_index = gene_index - 11
            if gene[gene_index] == primer[primer_index]:

Running this gives you:

  • 1
    $\begingroup$ Thanks a lot Christopher, I really appreciate that you took the time to help me out. I will definitely use this and potentially try to implement it in R. Thanks! $\endgroup$
    – Colin
    Sep 19 '18 at 22:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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