This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Generate random DNA sequences with variable %identity based on a reference sequence

Hello there!

I have a query sequence (~500 bases long) and I want to generate multiple sequences which have a specific % identity to this query sequence. For example, I would like to have 3 sequences which are 80%, 60%, and 20% similar to this query sequence. Are there any tools that could generate them?

Thank you!

alignment sequence

1 answer

This short python script will read a fasta file and randomly mutate the requested fraction of nucleotides. Save it as mutate_dna.py and run:

python mutate_dna.py dna.fa 0.2 > mutated_dna.fa

You will need to run it multiple times with different mutation rates to get what you want.

Starting with this sequence:

>dna
AGGGATCAGCGGATGTTGCTTATAGGACTCCGCTGGCACCTTATGAGAAA
TCAAAGTTTTTGGGTTCCGGgGGGaGtATGGTCGCaAGGcTGAAACTTAA
AGGAATTGACGGAAGGGCACCACCAGGAGTGGAGCCTGCGGCTTAATTTG
ACGCAACACGGGGAAACTTACCAGGTCCAGACATAGTAAGGATTGACAGA
CTGAGAGCTCTTTCTTGATTCTATGGGTGGTGGTGCATGGCCGTTCTTAG
TTGGTGGAGCGATTTGTCTGGTTAATTCCGTTAACGAACGAGACCTCAGC
CTGCTAACTAGCTACGTGGAGGCATCCCTTCACGGCCGGCTTCTTAGAGG
GACTATGGCCGTTTAGGCCAAGGAAGTTTGAGGCAATAACAGGTCTGTGA
TGCCCTTAGATGTTCTGGGCCGCACGCGCGCTACACTGATGTATTCAACG
AGTTCACACCTTGGCCGACAGGCCCGGGTAATCTTTGAAATTTCATCGTG
ATGGGGATAGATCATTGCAATTGTTGGTCTTCAACGAGGAATTCCTAGTA
AGCGCGAGTCATCAGCTCGCGTTGACTACGTCCCTGCCCTTTGTACACAC
CGCCCGTCGCTCCTACCGATTGAATGATCCGGTGAAGTGTTCGGATCGCG
GCGACGTGGGTGGTTCGCCGCCCGCGACGTCGCGAGAAGTCCACTAAACC
TTATCATTTAGAGGAAGGAGAAGTCGTAACAAGGTTTCCGTAGGTGAACC
TGCGG

It makes this sequence:

>mutated_0.2_dna
TGGGATCAGCGGATGTTGCTTCTAGGTCCCCACTGGCACCTTATGAGAATGGAAAGTTTTTGGATTCCGG
gGGGaGtATGGTGGCaTGGcTGATACTTAAATGAATTGTCGGGAGGGCTCCACCCGGAGTGGATCTTGCG
TCTTAATTCGACTTAAGAGGGGGAACCTTACCCGGTCCAGACAGAGTAAGTATTGACAAACTGAGGGCTC
TTTCTAGAAACCATGGGTGGTGCAGCATGCTCGTTCTTTGTTCGTGGAGCGATTTGTCTGGTTAATTTCC
TTGCCGAACGAGACCTTAGCCAGCTTACTAGCGACGTAGAGCCGTCCCTGGACGTCGGGCATCTTAGAAG
GACTAAGGCCGTTTGGACCAAGGAAGTTGAAGGCATTAACAGGTCTGTGATCACCATCGATGTTCGGGGC
CGCACCCGCGCTACTCTGATGTATTTAACGAGTTCCCGTCCTGGCACTCTGGCCCAGGTAATCCTAGATT
TTTCATCATGCTGGGGTTAAAGCATTCCAGTTGTAGGTCTTATGCGAGGGATTTCTAGTAAGCACGAGTC
TTCAGTTCACGTAGACTATGTCCTTGCCCGTTCTACACACCACCCGGCGCTTCTACCGTATGAATGTTCG
ACTGAAGTGTTCGCATCTCGGCGACGTGGGCGGTTCGGCGCCCCCGACGTCGCTAGAAGTCTACGAAACC
TTATCATCTAAAAGAAGGACAAGTCGTATAAAGGCTTCCGTAGGTGAGCCTGCGG

import os
import sys
from random import randint
import textwrap
import numpy as np
from Bio import SeqIO
from Bio.Seq import MutableSeq

def mutate(letter):
    pos = randint(0, 2)
    if (letter == 'A') or (letter == 'a'):
        bases = ['C', 'G', 'T']
    if (letter == 'C') or (letter == 'c'):
        bases = ['A', 'G', 'T']
    if (letter == 'G') or (letter == 'g'):
        bases = ['A', 'C', 'T']
    if (letter == 'T') or (letter == 't'):
        bases = ['A', 'C', 'G']
    return str(bases[pos])

if os.access(sys.argv[1], os.R_OK):
    FastaFile = open(sys.argv[1], 'rU')
else:
    print('\n !!! Input file "%s" does not exist in this directory !!!\n' %
          sys.argv[1])
    sys.exit(1)

mut_rate = float(sys.argv[2])
if (mut_rate > 1) or (mut_rate < 0):
    print('\n Sequence mutation rate must be in 0-1 range !\n')
    sys.exit(1)

for rec in SeqIO.parse(FastaFile, 'fasta'):
    name = rec.id
    seq = MutableSeq(list(rec.seq))
    seq_len = len(rec)
    mut_len = int(seq_len * mut_rate)
    mut_pos = np.random.randint(0, seq_len, mut_len)
    for x in range(mut_len):
        seq[mut_pos[x]] = mutate(seq[mut_pos[x]])
    rec.seq = seq
    print('>mutated_%s_%s' % (str(mut_rate), name))
    for lines in textwrap.wrap(str(rec.seq)):
        print(lines)

FastaFile.close()

Log in to answer this question.