[RC] Bioinformatics/Sequence mutation

[tsh73]

rosettacode.org/wiki/Bioinformatics/Sequence_mutation

following post by tenochtitlanuk
libertybasiccom.proboards.com/post/12332/thread

run example

  0: CCAGAGCTCCTTAAAGCCATACGCGCTTTGGGATGCAAGAATGCGGCAGG
 50: GTTTAGGACGGAGTATGTAATTGCAGCGTTCGAGCCGTTGGACTTTACCT
100: TTCAATTTCGTCGGGTTCACAGCCGATGCCCTTTAGGAAGGCGAAGCCAC
150: ACCTGGCTATAGTGAGAGTATTACCCGATAACCCGATTTCCCAGCCGTTC

BASECOUNT:
  A: 46
  C: 51
  G: 53
  T: 50
TOT= 200

Mutations
1 At pos  93 D   - C
2 At pos  14 S   A -> A
3 At pos 108 D   - C
4 At pos 137 D   - A
5 At pos 158 D   - A
6 At pos 127 D   - C
7 At pos  29 I   + A
8 At pos 137 S   G -> G
9 At pos 195 S   T -> A
10 At pos 161 D   - A

  0: CACAGAGCTCCTTAAAGCCATACGCGCTTTGGGATGCAAGAATGCGGCAG
 50: GGTTTAGGACGGAGTATGTAATTGCAGCGTTCGAGCCGTTGGATTTACCT
100: TTCAATTTGTCGGGTTCACAGCCGATGCCTTTAGGAGGCGAAGCCACACC
150: TGGCTATGTGGAGTATTACCCGATAACCCGATTTCCCAGCCGTAC

BASECOUNT:
  A: 45
  C: 48
  G: 53
  T: 49
TOT= 195


code
(and I think single select case is enough. But I use character bases as index point for counter, probably that could use another select case instead)
'http://rosettacode.org/wiki/Bioinformatics/Sequence_mutation

global bases$
bases$ = "ACGT"

'generate, once
seq$=""
for i = 1 to 200
   seq$=seq$+randBase$()
next
'print seq$

gosub [prettyPrt]
gosub [baseCount]

print
print "Mutations"
N=10    'mutations
mt$="SDI"
for i = 1 to N
   rPos = int(rnd(0)*len(seq$))+1
   print using("##", i);" At pos ";using("###", rPos);" " ;
   r = int(rnd(0)*3)+1
   r$=mid$(mt$, r, 1)
   print r$;"   ";
   select case r$
   case "S"
       r$ = randBase$()
       print mid$(seq$, rPos, 1);" -> ";r$
       seq$=left$(seq$, rPos-1) + r$ + mid$(seq$, rPos+1)
   case "D"
       print "- "; mid$(seq$, rPos, 1)
       seq$=left$(seq$, rPos-1) + mid$(seq$, rPos+1)
   case "I"
       r$ = randBase$()
       print "+ ";r$
       seq$=left$(seq$, rPos-1) + r$ + mid$(seq$, rPos)
   end select
next

gosub [prettyPrt]
gosub [baseCount]


end

[prettyPrt]
   for i = 1 to len(seq$)
       if (i-1) mod 50 =0 then print: print using("####", i-1);": ";
       print mid$(seq$, i, 1);
   next
   print
   return

[baseCount]
   dim count(len(bases$))    'clear array
   for i = 1 to len(seq$)
       c$=mid$(seq$, i, 1)
       pos = instr(bases$, c$)
       count(pos)=count(pos)+1
   next
   print
   print " BASECOUNT:"
   for i = 1 to len(bases$)
       print "   ";mid$(bases$, i, 1);": "; count(i)
   next
   print " TOT= ";len(seq$)
return

function randBase$()
   r = int(rnd(0)*len(bases$))+1
   randBase$= mid$(bases$, r, 1)
end function