Sequence operations¶

Modeled on the BioPython sequence tutorial.

A sequence is a series of letters representing an organism's DNA, RNA, or protein. In BioKotlin these are NucSeq and ProteinSeq, both of which implement Seq.

%use biokotlin
import biokotlin.seq.*

Creating a sequence¶

The top-level Seq function infers a nucleotide sequence. DNA is inferred here rather than RNA because of the T.

val seq = Seq("GCAGAT")
seq

GCAGAT

Complement the sequence:

seq.complement()

CGTCTA

Translate it to protein. The first call also prints the codon table ID, which defaults to 1.

seq.translate()

1
AD

Transcribe it to RNA:

seq.transcribe()

GCAGAU

Protein sequences¶

ProteinSeq is a distinct type, so the compiler will not let you mix it with nucleotide sequences. Count the glycines in a peptide with the AminoAcid enum rather than a bare string.

val proseq = ProteinSeq("GCAGAT")
val gCount = proseq.count(AminoAcid.G)
println("${AminoAcid.G.name3letter} count is $gCount")

Gly count is 2

The + operator concatenates two sequences of the same type:

val proseq2 = ProteinSeq("RMFGVE")
proseq2 + proseq

RMFGVEGCAGAT

DNA and RNA¶

DNA and RNA share one backing store and differ only in how they are viewed. Name the alphabet explicitly when there is no T or U to infer from.

val dna = NucSeq("GCTA")                  // inferred DNA
val rna = NucSeq("GCUA")                  // inferred RNA
val rnaSpecified = NucSeq("GCACCCCC", NUC.RNA)

println(dna.transcribe() == rna)          // true
println(dna)                              // GCTA
println(dna.repr())                       // NucSeqByte('GCTA',[A,C,G,T])

true
GCTA
NucSeqByte('GCTA',[A, C, G, T])

Working at scale¶

Sequences can be repeated with * and searched with count. Unambiguous DNA is stored at two bits per base, so even large sequences stay compact.

val bigSeq = seq * 100_000
bigSeq.count(Seq("TGC"))

Counting 6 bp palindromes is one of the operations in the benchmarks. Here it is on a one million base sequence.

val startTime = System.currentTimeMillis()
val palSeq = Seq("GATATCC") * 150_000
var totalCount = 0
var count = 0
for (i in 0..(palSeq.size() - 7)) {
    totalCount++
    val site = palSeq[i..(i + 5)]
    if (site == site.reverse_complement()) count++
}
println("count=$count, totalCount=$totalCount")
println("elapsed ${System.currentTimeMillis() - startTime}ms")

count=150000, totalCount=1049994
elapsed 242ms

Translation over the same sequence:

val startTime2 = System.currentTimeMillis()
val pro = palSeq.translate()
println(pro[0..2])
println("translated ${palSeq.size()} bases in ${System.currentTimeMillis() - startTime2}ms")

DIR
translated 1050000 bases in 31ms