Aligning sequence and comparing it against primer

Question

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.

Thanks!

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 = "")
        }else{
          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!

Answer 1

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:

ATGACCATACGATGTCGAAGCAAATATCAATTGAATCGTAGGG
ATGAGCATACCATGTGAAAGCTAATATCAATTGATTCGAACGC

Here is the code that gives you the expected output:

#!/bin/env python
import sys
primer = 'ATGTGGAAGCAAATATCAAATGA'

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

Running this gives you:

-----------ATGTGGAAGCAAATATCAAATGA---------
ATGACCATACG....C..............T...ATCGTAGGG
ATGAGCATACC.....A....T........T...TTCGAACGC