Ytdyklly

Goal of the program

The goal of the program is to translate a nucleic acid sequence into its corresponding amino acid sequence. The nucleic sequences
have to be formatted in a specific format called fasta. There is an existing implementation of this program in C here: emboss transeq

Fasta format

A fasta file looks like this:

>sequenceId comment

nucleic sequence

for example:

>Seq1 [organism=Carpodacus mexicanus] [clone=6b] actin (act) mRNA, partial cds

CCTTTATCTAATCTTTGGAGCATGAGCTGGCATAGTTGGAACCGCCCTCAGCCTCCTCATCCGTGCAGAA

CCCAGTCCTGTACCAACACCTCTTCTGATTCTTCGGCCATCCAGAAGTCTATATCCTCATTTTAC



>Seq2 [organism=uncultured bacillus sp.] [isolate=A2] corticotropin (CT) gene

GGTAGGTACCGCCCTAAGNCTCCTAATCCGAGCAGAACTANGCCAACCCGGAGCCCTTCTGGGAGACGAC

AATCAACATAAAA

A nucleic sequence is a string composed of the letters A, C, T, G, U, N

Expected output

A combinaison of 3 nucleic acid, named codon, gives a specif amino acid, for exemple GCT is the code for Alanine, symbolised by the letter A

With the Seq1 define above:

codon      CCT TTA TCT AAT CTT TGG AGC ATG ...

amino acid  P   L   S   N   L   W   S   M  ...

The expected output is:

>Seq1_1 [organism=Carpodacus mexicanus] [clone=6b] actin (act) mRNA, partial cds

PLSNLWSMSWHSWNRPQPPHPCRTQSCTNTSSDSSAIQKSISSFY

>Seq2_1 [organism=uncultured bacillus sp.] [isolate=A2] corticotropin (CT) gene

GRYRPKXPNPSRTXPTRSPSGRRQST*X

Specific rules

The program takes 4 parameters:

• 
  clean : if true, write STOP codon as X instead of *
  


• 
  trim : if true, remove all X and * chars from the right end of the amino - acid sequence


• 
  alternative: different way to compute reverse frame


• 
  frame: a string value in ["1", "2", "3", "F", "-1", "-2", "-3", "R", "6" ]
  

The frame define the position in the nucleic sequence to start from:

-frame 1

|-frame 2

||-frame 3

|||

CCTTTATCTAATCTTTGGAGCATGAGCTGGCATAGTTGGAACCGCCCTCAGCCTCCTCATCCGTGCAGAA

CCCAGTCCTGTACCAACACCTCTTCTGATTCTTCGGCCATCCAGAAGTCTATATCCTCATTTTAC

                                                              |||

                                                              ||-frame -1

                                                              |-frame -2

                                                              -frame -3





frame "F" = "1", "2", "3"

frame "R" = "-1", "-2", "-3"

frame "6" = "1", "2", "3", "-1", "-2", "-3"

In the output file, we add the frame used to the sequenceId: sequenceId = sequenceId_frame

For example if the program is used with frame=6, the expected output is

>Seq1_1 [organism=Carpodacus mexicanus] [clone=6b] actin (act) mRNA, partial cds

PLSNLWSMSWHSWNRPQPPHPCRTQSCTNTSSDSSAIQKSISSFY

>Seq1_2 [organism=Carpodacus mexicanus] [clone=6b] actin (act) mRNA, partial cds

LYLIFGA*AGIVGTALSLLIRAEPSPVPTPLLILRPSRSLYPHFT

>Seq1_3 [organism=Carpodacus mexicanus] [clone=6b] actin (act) mRNA, partial cds

FI*SLEHELA*LEPPSASSSVQNPVLYQHLF*FFGHPEVYILILX

>Seq1_4 [organism=Carpodacus mexicanus] [clone=6b] actin (act) mRNA, partial cds

VK*GYRLLDGRRIRRGVGTGLGSARMRRLRAVPTMPAHAPKIR*R

>Seq1_5 [organism=Carpodacus mexicanus] [clone=6b] actin (act) mRNA, partial cds

KMRI*TSGWPKNQKRCWYRTGFCTDEEAEGGSNYASSCSKD*IKX

>Seq1_6 [organism=Carpodacus mexicanus] [clone=6b] actin (act) mRNA, partial cds

*NEDIDFWMAEESEEVLVQDWVLHG*GG*GRFQLCQLMLQRLDKG

>Seq2_1 [organism=uncultured bacillus sp.] [isolate=A2] corticotropin (CT) gene

GRYRPKXPNPSRTXPTRSPSGRRQST*X

>Seq2_2 [organism=uncultured bacillus sp.] [isolate=A2] corticotropin (CT) gene

VGTALXLLIRAELXQPGALLGDDNQHKX

>Seq2_3 [organism=uncultured bacillus sp.] [isolate=A2] corticotropin (CT) gene

*VPP*XS*SEQNXANPEPFWETTINIK

>Seq2_4 [organism=uncultured bacillus sp.] [isolate=A2] corticotropin (CT) gene

LC*LSSPRRAPGWXSSARIRXLRAVPT

>Seq2_5 [organism=uncultured bacillus sp.] [isolate=A2] corticotropin (CT) gene

FMLIVVSQKGSGLA*FCSD*EX*GGTYX

>Seq2_6 [organism=uncultured bacillus sp.] [isolate=A2] corticotropin (CT) gene

FYVDCRLPEGLRVGXVLLGLGXLGRYLP

Improvements

• Performance improvements


• Is the code easy to read / understand ?


• Is the code idiomatic ?

( full project with tests + flags handling is avaible on github: gotranseq)

code :

package transeq



import (

    "bufio"

    "bytes"

    "context"

    "encoding/binary"

    "fmt"

    "io"

    "runtime"

    "sync"

)



var (

    letterCode = map[byte]uint8{

        'A': aCode,

        'C': cCode,

        'T': tCode,

        'G': gCode,

        'N': nCode,

        'U': uCode,

    }

    standard = map[string]byte{

        "TTT": 'F',

        "TCT": 'S',

        "TAT": 'Y',

        "TGT": 'C',

        "TTC": 'F',

        "TCC": 'S',

        "TAC": 'Y',

        "TGC": 'C',

        "TTA": 'L',

        "TCA": 'S',

        "TAA": '*',

        "TGA": '*',

        "TTG": 'L',

        "TCG": 'S',

        "TAG": '*',

        "TGG": 'W',

        "CTT": 'L',

        "CCT": 'P',

        "CAT": 'H',

        "CGT": 'R',

        "CTC": 'L',

        "CCC": 'P',

        "CAC": 'H',

        "CGC": 'R',

        "CTA": 'L',

        "CCA": 'P',

        "CAA": 'Q',

        "CGA": 'R',

        "CTG": 'L',

        "CCG": 'P',

        "CAG": 'Q',

        "CGG": 'R',

        "ATT": 'I',

        "ACT": 'T',

        "AAT": 'N',

        "AGT": 'S',

        "ATC": 'I',

        "ACC": 'T',

        "AAC": 'N',

        "AGC": 'S',

        "ATA": 'I',

        "ACA": 'T',

        "AAA": 'K',

        "AGA": 'R',

        "ATG": 'M',

        "ACG": 'T',

        "AAG": 'K',

        "AGG": 'R',

        "GTT": 'V',

        "GCT": 'A',

        "GAT": 'D',

        "GGT": 'G',

        "GTC": 'V',

        "GCC": 'A',

        "GAC": 'D',

        "GGC": 'G',

        "GTA": 'V',

        "GCA": 'A',

        "GAA": 'E',

        "GGA": 'G',

        "GTG": 'V',

        "GCG": 'A',

        "GAG": 'E',

        "GGG": 'G',

    }

)



const (

    nCode = uint8(0)

    aCode = uint8(1)

    cCode = uint8(2)

    tCode = uint8(3)

    uCode = uint8(3)

    gCode = uint8(4)



    stopByte = '*'

    unknown  = 'X'

    // Length of the array to store code/bytes

    // uses gCode because it's the biggest uint8 of all codes

    arrayCodeSize = (uint32(gCode) | uint32(gCode)<<8 | uint32(gCode)<<16) + 1

)



func createCodeArray(clean bool) byte {



    resultMap := map[uint32]byte{}

    twoLetterMap := map[string]byte{}



    tmpCode := make(uint8, 4)



    for codon, aaCode := range standard {

        // generate 3 letter code

        for i := 0; i < 3; i++ {

            tmpCode[i] = letterCode[codon[i]]

        }

        // each codon is represented by an unique uint32:

        // each possible nucleotide is represented by an uint8 (255 possibility)

        // the three first bytes are the the code for each nucleotide

        // last byte is unused ( eq to uint8(0) )

        // example:

        // codon 'ACG' ==> uint32(aCode) | uint32(cCode)<<8 | uint32(gCode)<<16

        uint32Code := uint32(tmpCode[0]) | uint32(tmpCode[1])<<8 | uint32(tmpCode[2])<<16

        resultMap[uint32Code] = aaCode



        // generate 2 letter code

        codes, ok := twoLetterMap[codon[0:2]]

        if !ok {

            twoLetterMap[codon[0:2]] = byte{aaCode}

        } else {

            twoLetterMap[codon[0:2]] = append(codes, aaCode)

        }

    }

    for twoLetterCodon, codes := range twoLetterMap {

        uniqueAA := true

        for i := 0; i < len(codes); i++ {



            if codes[i] != codes[0] {

                uniqueAA = false

            }

        }

        if uniqueAA {

            first := letterCode[twoLetterCodon[0]]

            second := letterCode[twoLetterCodon[1]]



            uint32Code := uint32(first) | uint32(second)<<8

            resultMap[uint32Code] = codes[0]

        }

    }

    // if clean is specified, we want to replace all '*' by 'X' in the output

    // sequence, so replace all occurrences of '*' directly in the ref map

    if clean {

        for k, v := range resultMap {

            if v == stopByte {

                resultMap[k] = unknown

            }

        }

    }

    r := make(byte, arrayCodeSize)

    for k, v := range resultMap {

        r[k] = v

    }

    return r

}



func computeFrames(frameName string) (frames int, reverse bool, err error) {



    frames = make(int, 6)

    reverse = false



    switch frameName {

    case "1":

        frames[0] = 1

    case "2":

        frames[1] = 1

    case "3":

        frames[2] = 1

    case "F":

        for i := 0; i < 3; i++ {

            frames[i] = 1

        }

    case "-1":

        frames[3] = 1

        reverse = true

    case "-2":

        frames[4] = 1

        reverse = true

    case "-3":

        frames[5] = 1

        reverse = true

    case "R":

        for i := 3; i < 6; i++ {

            frames[i] = 1

        }

        reverse = true

    case "6":

        for i := range frames {

            frames[i] = 1

        }

        reverse = true

    default:

        err = fmt.Errorf("wrong value for -f | --frame parameter: %s", frameName)

    }

    return frames, reverse, err

}



type writer struct {

    buf            *bytes.Buffer

    currentLineLen int

    bytesToTrim    int

}



func (w *writer) addByte(b byte) {

    w.buf.WriteByte(b)

    w.currentLineLen++

    if b == stopByte || b == unknown {

        w.bytesToTrim++

    } else {

        w.bytesToTrim = 0

    }

}



func (w *writer) addUnknown() {

    w.buf.WriteByte(unknown)

    w.currentLineLen++

    w.bytesToTrim++

}



func (w *writer) newLine() {

    w.buf.WriteByte('n')

    w.currentLineLen = 0

    w.bytesToTrim++

}



const (

    // size of the buffer for writing to file

    maxBufferSize = 1024 * 1024 * 30

    // max line size for sequence

    maxLineSize = 60

    // suffixes ta add to sequence id for each frame

    suffixes = "123456"

)



// Translate read a fata file, translate each sequence to the corresponding prot sequence in the specified frame

func Translate(inputSequence io.Reader, out io.Writer, frame string, clean, trim, alternative bool) error {



    arrayCode := createCodeArray(clean)

    framesToGenerate, reverse, err := computeFrames(frame)

    if err != nil {

        return err

    }



    fnaSequences := make(chan encodedSequence, 10)

    errs := make(chan error, 1)



    var wg sync.WaitGroup

    ctx, cancel := context.WithCancel(context.Background())

    defer cancel()



    for nWorker := 0; nWorker < runtime.NumCPU(); nWorker++ {



        wg.Add(1)



        go func() {



            defer wg.Done()



            startPosition := make(int, 3)

            w := &writer{

                buf:            bytes.NewBuffer(nil),

                bytesToTrim:    0,

                currentLineLen: 0,

            }



            for sequence := range fnaSequences {



                select {

                case <-ctx.Done():

                    return

                default:

                }



                frameIndex := 0

                startPosition[0], startPosition[1], startPosition[2] = 0, 1, 2



                idSize := int(binary.LittleEndian.Uint32(sequence[0:4]))

                nuclSeqLength := len(sequence) - idSize



            Translate:

                for _, startPos := range startPosition {



                    if framesToGenerate[frameIndex] == 0 {

                        frameIndex++

                        continue

                    }



                    // sequence id should look like

                    // >sequenceID_<frame> comment

                    idEnd := bytes.IndexByte(sequence[4:idSize], ' ')

                    if idEnd != -1 {

                        w.buf.Write(sequence[4 : 4+idEnd])

                        w.buf.WriteByte('_')

                        w.buf.WriteByte(suffixes[frameIndex])

                        w.buf.Write(sequence[4+idEnd : idSize])

                    } else {

                        w.buf.Write(sequence[4:idSize])

                        w.buf.WriteByte('_')

                        w.buf.WriteByte(suffixes[frameIndex])

                    }

                    w.newLine()



                    // if in trim mode, nb of bytes to trim (nb of successive 'X', '*' and 'n'

                    // from right end of the sequence)

                    w.bytesToTrim = 0

                    w.currentLineLen = 0



                    // read the sequence 3 letters at a time, starting at a specific position

                    // corresponding to the frame

                    for pos := startPos + 2 + idSize; pos < len(sequence); pos += 3 {



                        if w.currentLineLen == maxLineSize {

                            w.newLine()

                        }

                        // create an uint32 from the codon, to retrieve the corresponding

                        // AA from the map

                        codonCode := uint32(sequence[pos-2]) | uint32(sequence[pos-1])<<8 | uint32(sequence[pos])<<16



                        b := arrayCode[codonCode]

                        if b != byte(0) {

                            w.addByte(b)

                        } else {

                            w.addUnknown()

                        }

                    }



                    // the last codon is only 2 nucleotid long, try to guess

                    // the corresponding AA

                    if (nuclSeqLength-startPos)%3 == 2 {



                        if w.currentLineLen == maxLineSize {

                            w.newLine()

                        }

                        codonCode := uint32(sequence[len(sequence)-2]) | uint32(sequence[len(sequence)-1])<<8



                        b := arrayCode[codonCode]

                        if b != byte(0) {

                            w.addByte(b)

                        } else {

                            w.addUnknown()

                        }

                    }



                    // the last codon is only 1 nucleotid long, no way to guess

                    // the corresponding AA

                    if (nuclSeqLength-startPos)%3 == 1 {

                        if w.currentLineLen == maxLineSize {

                            w.newLine()

                        }

                        w.addUnknown()

                    }



                    if trim && w.bytesToTrim > 0 {

                        // remove the last bytesToTrim bytes of the buffer

                        // as they are 'X', '*' or 'n'

                        w.buf.Truncate(w.buf.Len() - w.bytesToTrim)

                        w.currentLineLen -= w.bytesToTrim

                    }



                    if w.currentLineLen != 0 {

                        w.newLine()

                    }

                    frameIndex++

                }



                if reverse && frameIndex < 6 {



                    // get the complementary sequence.

                    // Basically, switch

                    //   A <-> T

                    //   C <-> G

                    // N is not modified

                    for i, n := range sequence[idSize:] {



                        switch n {

                        case aCode:

                            sequence[i+idSize] = tCode

                        case tCode:

                            // handle both tCode and uCode

                            sequence[i+idSize] = aCode

                        case cCode:

                            sequence[i+idSize] = gCode

                        case gCode:

                            sequence[i+idSize] = cCode

                        default:

                            //case N -> leave it

                        }

                    }

                    // reverse the sequence

                    for i, j := idSize, len(sequence)-1; i < j; i, j = i+1, j-1 {

                        sequence[i], sequence[j] = sequence[j], sequence[i]

                    }



                    if !alternative {

                        // Staden convention: Frame -1 is the reverse-complement of the sequence

                        // having the same codon phase as frame 1. Frame -2 is the same phase as

                        // frame 2. Frame -3 is the same phase as frame 3

                        //

                        // use the matrix to keep track of the forward frame as it depends on the

                        // length of the sequence

                        switch nuclSeqLength % 3 {

                        case 0:

                            startPosition[0], startPosition[1], startPosition[2] = 0, 2, 1

                        case 1:

                            startPosition[0], startPosition[1], startPosition[2] = 1, 0, 2

                        case 2:

                            startPosition[0], startPosition[1], startPosition[2] = 2, 1, 0

                        }

                    }

                    // run the same loop, but with the reverse-complemented sequence

                    goto Translate

                }



                if w.buf.Len() > maxBufferSize {

                    _, err := out.Write(w.buf.Bytes())

                    if err != nil {

                        select {

                        case errs <- fmt.Errorf("fail to write to output file: %v", err):

                        default:

                        }

                        cancel()

                        return

                    }

                    w.buf.Reset()

                }

                pool.Put(sequence)

            }



            if w.buf.Len() > 0 {

                _, err := out.Write(w.buf.Bytes())

                if err != nil {

                    select {

                    case errs <- fmt.Errorf("fail to write to output file: %v", err):

                    default:

                    }

                    cancel()

                    return

                }

            }

        }()

    }

    readSequenceFromFasta(ctx, inputSequence, fnaSequences)



    wg.Wait()

    select {

    case err, ok := <-errs:

        if ok {

            return err

        }

    default:

    }

    return nil

}



func readSequenceFromFasta(ctx context.Context, inputSequence io.Reader, fnaSequences chan encodedSequence) {



    feeder := &fastaChannelFeeder{

        idBuffer:       bytes.NewBuffer(nil),

        commentBuffer:  bytes.NewBuffer(nil),

        sequenceBuffer: bytes.NewBuffer(nil),

        fastaChan:      fnaSequences,

    }

    // fasta format is:

    //

    // >sequenceID some comments on sequence

    // ACAGGCAGAGACACGACAGACGACGACACAGGAGCAGACAGCAGCAGACGACCACATATT

    // TTTGCGGTCACATGACGACTTCGGCAGCGA

    //

    // see https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=BlastHelp

    // section 1 for details

    scanner := bufio.NewScanner(inputSequence)

Loop:

    for scanner.Scan() {



        line := scanner.Bytes()

        if len(line) == 0 {

            continue

        }

        if line[0] == '>' {



            if feeder.idBuffer.Len() > 0 {

                select {

                case <-ctx.Done():

                    break Loop

                default:

                }

                feeder.sendFasta()

            }

            feeder.reset()



            // parse the ID of the sequence. ID is formatted like this:

            // >sequenceID comments

            seqID := bytes.SplitN(line, byte{' '}, 2)

            feeder.idBuffer.Write(seqID[0])



            if len(seqID) > 1 {

                feeder.commentBuffer.WriteByte(' ')

                feeder.commentBuffer.Write(seqID[1])

            }

        } else {

            // if the line doesn't start with '>', then it's a part of the

            // nucleotide sequence, so write it to the buffer

            feeder.sequenceBuffer.Write(line)

        }

    }

    // don't forget to push last sequence

    select {

    case <-ctx.Done():

    default:

        feeder.sendFasta()

    }

    close(fnaSequences)

}



// a type to hold an encoded fasta sequence

//

//  s[0:4] stores the size of the sequence id + the size of the comment as an uint32 (little endian)

//  s[4:idSize] stores the sequence id, and the comment id there is one

//  s[idSize:] stores the nucl sequence

type encodedSequence byte



var pool = sync.Pool{

    New: func() interface{} {

        return make(encodedSequence, 512)

    },

}



func getSizedSlice(idSize, requiredSize int) encodedSequence {

    s := pool.Get().(encodedSequence)

    binary.LittleEndian.PutUint32(s[0:4], uint32(idSize))



    for len(s) < requiredSize {

        s = append(s, byte(0))

    }

    return s[0:requiredSize]

}



func (f *fastaChannelFeeder) sendFasta() {



    idSize := 4 + f.idBuffer.Len() + f.commentBuffer.Len()

    requiredSize := idSize + f.sequenceBuffer.Len()



    s := getSizedSlice(idSize, requiredSize)



    if f.commentBuffer.Len() > 0 {

        copy(s[idSize-f.commentBuffer.Len():idSize], f.commentBuffer.Bytes())

    }



    copy(s[4:4+f.idBuffer.Len()], f.idBuffer.Bytes())



    // convert the sequence of bytes to an array of uint8 codes,

    // so a codon (3 nucleotides | 3 bytes ) can be represented

    // as an uint32

    for i, b := range f.sequenceBuffer.Bytes() {



        switch b {

        case 'A':

            s[i+idSize] = aCode

        case 'C':

            s[i+idSize] = cCode

        case 'G':

            s[i+idSize] = gCode

        case 'T', 'U':

            s[i+idSize] = tCode

        case 'N':

            s[i+idSize] = nCode

        default:

            fmt.Printf("WARNING: invalid char in sequence %s: %s, ignoring", s[4:4+idSize], string(b))

        }

    }

    f.fastaChan <- s

}



type fastaChannelFeeder struct {

    idBuffer, commentBuffer, sequenceBuffer *bytes.Buffer

    fastaChan                               chan encodedSequence

}



func (f *fastaChannelFeeder) reset() {

    f.idBuffer.Reset()

    f.sequenceBuffer.Reset()

    f.commentBuffer.Reset()

}