counting motifs in dna sequence

Question

[aguilar.juanpablo]

Hi,

I have these sequences:

GCAGGCATAGTCGGAACTGCTCTAAGCCTATTAATTCGAGCTGAGCTAAGCCAGCCTGGGGCTCTGCTCGGAGATGA AGTGGGCTTGTTGGGACTGGTCTTTCTTTATTAATTCGTTTTGAGTTAGGCACTGTTGGAGTTTTATTAG---ATAA GCAGGAATAGTTGGAACCGCCCTTAGCTTATTAATTCGAGCAGAACTCAGCCAACCTGGTGCCTTATTAGGGGATGA GCTGGCATAGTAGGAACTGCCCTTAGCCTTTTAATTCGAGCAGAGCTCAGTCAACCCGGAGCCCTGCTCGGAGATGA GCAGGAATAGTTGGAACTGCACTAAGCCTTTTAATTCGAGCTGAACTAAGCCAACCCGGAGCATTACTTGGAGACGA

They would be actually longer, but right not it is not important.

I would like to estimate a given value for the sequences given the number of motif/s in. I would like to count a (number of) motif/s like "ATCGCGCGCGCTTTAAA" in a given sequence, and then use that number to estimate a value for that sequence. I know that you can use a logical question to ask whether a given sequence has a motif but I would like to count them.

Thanks

[GouthamAtla]

Do you have actual motif or just a string, ATCGCGCGCGCTTTAAA ?

[Macherki M E]

Biostrings package in R provides a beautiful function for you purpose:gregexpr2

Answer 1

[shenwei356]

The sequences you pasted seemed aligned sequences:

$ cat seqs.txt 
GCAGGCATAGTCGGAACTGCTCTAAGCCTATTAATTCGAGCTGAGCTAAGCCAGCCTGGGGCTCTGCTCGGAGATGA
AGTGGGCTTGTTGGGACTGGTCTTTCTTTATTAATTCGTTTTGAGTTAGGCACTGTTGGAGTTTTATTAG---ATAA
GCAGGAATAGTTGGAACCGCCCTTAGCTTATTAATTCGAGCAGAACTCAGCCAACCTGGTGCCTTATTAGGGGATGA
GCTGGCATAGTAGGAACTGCCCTTAGCCTTTTAATTCGAGCAGAGCTCAGTCAACCCGGAGCCCTGCTCGGAGATGA
GCAGGAATAGTTGGAACTGCACTAAGCCTTTTAATTCGAGCTGAACTAAGCCAACCCGGAGCATTACTTGGAGACGA

Let's do this using seqkit locate (usage) and csvtk freq (usage) in one-line command:

$ cat -n seqs.txt | sed 's/ //g' | seqkit tab2fx | seqkit seq -g | seqkit locate -i -d -p ATC | csvtk -t freq -k -f 1
seqID   frequency
1       1
2       1
3       1
4       1

Both SeqKit and csvtk provide executable binary files for Linux/Windows/Mac OS, so you can just download and run without any compilation and configuration.

Let me explain this step by step:

Step 1: converting your input sequences to FASTA format and remove gaps (-):

$ cat -n seqs.txt | sed 's/ //g' | seqkit tab2fx | seqkit seq -g 
>1
GCAGGCATAGTCGGAACTGCTCTAAGCCTATTAATTCGAGCTGAGCTAAGCCAGCCTGGG
GCTCTGCTCGGAGATGA
>2
AGTGGGCTTGTTGGGACTGGTCTTTCTTTATTAATTCGTTTTGAGTTAGGCACTGTTGGA
GTTTTATTAGATAA

Step 2: searching and locating your sequence motif (e.g., ATC):

$ cat -n seqs.txt | sed 's/ //g' | seqkit tab2fx | seqkit seq -g | seqkit locate -i -d -p ATC | csvtk -t pretty
seqID   patternName   pattern   strand   start   end   matched
1       ATC           ATC       -        73      75    ATC
2       ATC           ATC       -        70      72    ATC
3       ATC           ATC       -        73      75    ATC
4       ATC           ATC       -        73      75    ATC

Notes:

• Flag -d/--degenerate is on, so you can search motif containing degenerate base, e.g., ACTGNNNNCCC.
• Command csvtk -t pretty is used to align table data.

If you have lots of sequence motifs, you can save them to file in FASTA format, and use this:

$ cat motifs.fa
>motif1
ATC
>motif2
actgn


$ cat -n seqs.txt | sed 's/ //g' | seqkit tab2fx | seqkit seq -g | seqkit locate -i -d -f motifs.fa | csvtk -t pretty
seqID   patternName   pattern   strand   start   end   matched
1       motif2        actgn     +        16      20    ACTGC
1       motif1        ATC       -        73      75    ATC
2       motif2        actgn     +        16      20    ACTGG
2       motif2        actgn     +        52      56    ACTGT
2       motif1        ATC       -        70      72    ATC
3       motif1        ATC       -        73      75    ATC
4       motif1        ATC       -        73      75    ATC
4       motif2        actgn     +        16      20    ACTGC
4       motif2        actgn     -        47      51    ACTGA
5       motif2        actgn     +        16      20    ACTGC

Step 3: counting the number of motifs in every sequences with csvtk freq

$  cat -n seqs.txt | sed 's/ //g' | seqkit tab2fx | seqkit seq -g | seqkit locate -p ATC | csvtk -t freq -k -f 1
seqID   frequency
1       1
2       1
3       1
4       1

For case of multiple motifs:

$ cat -n seqs.txt | sed 's/ //g' | seqkit tab2fx | seqkit seq -g | seqkit locate -i -d -f motifs.fa | csvtk -t freq -k -f 1,2 | csvtk -t pretty
seqID   patternName   frequency
1       motif1        1
1       motif2        1
2       motif1        1
2       motif2        2
3       motif1        1
4       motif1        1
4       motif2        2
5       motif2        1

[shenwei356]

Please download the latest version (v0.4.0 +). I've just fixed a bug of seqkit locate. e.g., only find two locations (1-4, 7-10, missing 4-7) of ACGA in ACGACGACGA.

$ echo -ne ">seq\nACGACGACGA\n" | seqkit locate -i -p acga | column -t
seqID  patternName  pattern  strand  start  end  matched
seq    acga         acga     +       1      4    ACGA
seq    acga         acga     +       4      7    ACGA   ## missing in previous version
seq    acga         acga     +       7      10   ACGA

Answer 2

[Istvan Albert]

You may use the dreg emboss tool or seqkit grep and probably many others. Even grep will work once you remove the new lines

• http://emboss.sourceforge.net/apps/cvs/emboss/apps/dreg.html
• http://bioinf.shenwei.me/seqkit/usage/#grep

will put the whole genome on one line, so don't do it for large genomes:

 cat fasta.fa | tr -d '\n' | grep -o "ATGC"

[shenwei356]

Hi, Istvan, it should be seqkit locate not seqkit grep

[Istvan Albert]

Ah ok, thanks for the correction.