Bioinformatics/Sequence mutation: Difference between revisions

=={{header|Nim}}==

<lang Nim>import random

import strformat

import strutils

type

# Enumeration type for bases.

Base {.pure.} = enum A, C, G, T, Other = "other"

# Sequence of bases.

DnaSequence = seq[Base]

# Kind of mutation.

Mutation = enum mutSwap, mutDelete, mutInsert

#---------------------------------------------------------------------------------------------------

proc newDnaSeq(length: Natural): DnaSequence =

## Create a DNA sequence of given length.

result = newSeqOfCap[Base](length)

for _ in 1..length:

result.add(Base(rand(3)))

#---------------------------------------------------------------------------------------------------

proc mutate(dnaSeq: var DnaSequence) =

## Mutate a sequence (it is changed in place).

# Choose randomly the position of mutation.

let idx = rand(dnaSeq.high)

# Choose randomly the kind of mutation.

let mut = Mutation(rand(ord(Mutation.high)))

# Apply the mutation.

case mut

of mutSwap:

let newBase = Base(rand(ord(Base.Other) - 1))

echo fmt"Changing base at position {idx + 1} from {dnaSeq[idx]} to {newBase}"

dnaSeq[idx] = newBase

of mutDelete:

echo fmt"Deleting base {dnaSeq[idx]} at position {idx + 1}"

dnaSeq.delete(idx)

of mutInsert:

let newBase = Base(rand(ord(Base.Other) - 1))

echo fmt"Inserting base {newBase} at position {idx + 1}"

dnaSeq.insert(newBase, idx)

#---------------------------------------------------------------------------------------------------

proc display(dnaSeq: DnaSequence) =

## Display a DNA sequence using EMBL format.

var counts: array[Base, Natural] # Count of bases.

for base in dnaSeq:

inc counts[base]

# Display the SQ line.

var sqline = fmt"SQ {dnaSeq.len} BP; "

for (base, count) in counts.pairs:

sqline &= fmt"{count} {base}; "

echo sqline

# Display the sequence.

var idx = 0

var row = newStringOfCap(80)

var remaining = dnaSeq.len

while remaining > 0:

row.setLen(0)

row.add(" ")

# Add groups of 10 bases.

for group in 1..6:

let nextIdx = idx + min(10, remaining)

for i in idx..<nextIdx:

row.add($dnaSeq[i])

row.add(' ')

dec remaining, nextIdx - idx

idx = nextIdx

if remaining == 0:

break

# Append the number of the last base in the row.

row.add(spaces(72 - row.len))

row.add(fmt"{idx:>8}")

echo row

# Add termination.

echo "//"

#———————————————————————————————————————————————————————————————————————————————————————————————————

randomize()

var dnaSeq = newDnaSeq(200)

echo "Initial sequence"

echo "———————————————\n"

dnaSeq.display()

echo "\nMutations"

echo "—————————\n"

for _ in 1..10:

dnaSeq.mutate()

echo "\nMutated sequence"

echo "————————————————\n"

dnaSeq.display()</lang>

{{out}}

<pre>Initial sequence

———————————————

SQ 200 BP; 53 A; 52 C; 45 G; 50 T; 0 other;

ATGGATAGAA TGGCACCTCA GTGTCACATG TCCTAGGCAC ATTCTACGTA CTAGTTTCTG 60

GAGGTCGATA AATACAAGAT GGAATACTCT TATCAACCGC TAGCAATGAA ACTTAGATAG 120

CCACCCCTCG ATCGCGGTTC GCTATGCGGC ATCGTCAACT GCGTCAAAGC ACTACGTCGT 180

TTCGTGACTG CCAGTCGAGC 200

//

Mutations

—————————

Inserting base A at position 121

Changing base at position 101 from T to G

Deleting base A at position 115

Changing base at position 126 from C to G

Deleting base T at position 155

Deleting base G at position 198

Deleting base T at position 159

Changing base at position 144 from A to T

Inserting base C at position 34

Changing base at position 127 from G to G

Mutated sequence

————————————————

SQ 198 BP; 52 A; 52 C; 46 G; 48 T; 0 other;

ATGGATAGAA TGGCACCTCA GTGTCACATG TCCCTAGGCA CATTCTACGT ACTAGTTTCT 60

GGAGGTCGAT AAATACAAGA TGGAATACTC TTATCAACCG CGAGCAATGA AACTTGATAG 120

ACCACCGCTC GATCGCGGTT CGCTTTGCGG CATCGCAACG CGTCAAAGCA CTACGTCGTT 180

TCGTGACTGC CAGTCGAC 198

//</pre>