File: gkStoreEncode.C

package info (click to toggle)
canu 1.7.1+dfsg-1~bpo9+1
  • links: PTS, VCS
  • area: main
  • in suites: stretch-backports
  • size: 7,680 kB
  • sloc: cpp: 66,708; perl: 13,682; ansic: 4,020; makefile: 627; sh: 472; python: 39
file content (157 lines) | stat: -rw-r--r-- 4,213 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157

/******************************************************************************
 *
 *  This file is part of canu, a software program that assembles whole-genome
 *  sequencing reads into contigs.
 *
 *  This software is based on:
 *    'Celera Assembler' (http://wgs-assembler.sourceforge.net)
 *    the 'kmer package' (http://kmer.sourceforge.net)
 *  both originally distributed by Applera Corporation under the GNU General
 *  Public License, version 2.
 *
 *  Canu branched from Celera Assembler at its revision 4587.
 *  Canu branched from the kmer project at its revision 1994.
 *
 *  Modifications by:
 *
 *    Brian P. Walenz beginning on 2015-DEC-03
 *      are a 'United States Government Work', and
 *      are released in the public domain
 *
 *  File 'README.licenses' in the root directory of this distribution contains
 *  full conditions and disclaimers for each license.
 */

#include "gkStore.H"


//  Encode seq as 2-bit bases.  Doesn't touch qlt.
uint32
gkReadData::gkReadData_encode2bit(uint8 *&chunk, char *seq, uint32 seqLen) {

  //  Scan the read, if there are non-acgt, return length 0; this cannot encode it.

  for (uint32 ii=0; ii<seqLen; ii++) {
    char  base = seq[ii];

    if ((base != 'a') && (base != 'A') &&
        (base != 'c') && (base != 'C') &&
        (base != 'g') && (base != 'G') &&
        (base != 't') && (base != 'T'))
      return(0);
  }

  uint8  acgt[256] = { 0 };

  acgt['a'] = acgt['A'] = 0x00;
  acgt['c'] = acgt['C'] = 0x01;
  acgt['g'] = acgt['G'] = 0x02;
  acgt['t'] = acgt['T'] = 0x03;

  uint32 chunkLen = 0;

  chunk    = new uint8 [ seqLen / 4 + 1];

  for (uint32 ii=0; ii<seqLen; ) {
    uint8  byte = 0;

    if (ii + 4 < seqLen) {
      byte  = acgt[seq[ii++]];  byte <<= 2;
      byte |= acgt[seq[ii++]];  byte <<= 2;
      byte |= acgt[seq[ii++]];  byte <<= 2;
      byte |= acgt[seq[ii++]];
    }

    else {
      if (ii < seqLen)  byte |= acgt[seq[ii++]];   byte <<= 2;  //  The if here is redundant, but pretty.
      if (ii < seqLen)  byte |= acgt[seq[ii++]];   byte <<= 2;  //  Yes, everything shifts, not a mistake to leave out the braces.
      if (ii < seqLen)  byte |= acgt[seq[ii++]];   byte <<= 2;
      if (ii < seqLen)  byte |= acgt[seq[ii++]];
    }

    chunk[chunkLen++] = byte;
  }

  return(chunkLen);
}



bool
gkReadData::gkReadData_decode2bit(uint8 *chunk, uint32 chunkLen, char *seq, uint32 seqLen) {

  if (chunkLen == 0)
    return(false);

  uint32   chunkPos = 0;

  char     acgt[4] = { 'A', 'C', 'G', 'T' };

  for (uint32 ii=0; ii<seqLen; ) {
    assert(chunkPos < chunkLen);

    uint8  byte = chunk[chunkPos++];

    if (ii + 4 < seqLen) {
      seq[ii++] = acgt[((byte >> 6) & 0x03)];
      seq[ii++] = acgt[((byte >> 4) & 0x03)];
      seq[ii++] = acgt[((byte >> 2) & 0x03)];
      seq[ii++] = acgt[((byte >> 0) & 0x03)];
    }

    else {
      if (ii < seqLen)  seq[ii++] = acgt[((byte >> 6) & 0x03)];  //  This if is also redundant, and also pretty.
      if (ii < seqLen)  seq[ii++] = acgt[((byte >> 4) & 0x03)];
      if (ii < seqLen)  seq[ii++] = acgt[((byte >> 2) & 0x03)];
      if (ii < seqLen)  seq[ii++] = acgt[((byte >> 0) & 0x03)];
    }
  }

  seq[seqLen] = 0;

  return(true);
}



//  Encode seq as 3-bases-in-7-bits.  Doesn't touch qlt.
uint32
gkReadData::gkReadData_encode3bit(uint8 *&UNUSED(chunk), char *UNUSED(seq), uint32 UNUSED(seqLen)) {
  return(0);
}

bool
gkReadData::gkReadData_decode3bit(uint8 *UNUSED(chunk), uint32 UNUSED(chunkLen), char *UNUSED(seq), uint32 UNUSED(seqLen)) {
  return(false);
}





//  Encode qualities as 4 bit integers.  Doesn't touch seq.
uint32
gkReadData::gkReadData_encode4bit(uint8 *&UNUSED(chunk), uint8 *UNUSED(qlt), uint32 UNUSED(qltLen)) {
  return(0);
}

bool
gkReadData::gkReadData_decode4bit(uint8 *UNUSED(chunk), uint32 UNUSED(chunkLen), uint8 *UNUSED(qlt), uint32 UNUSED(qltLen)) {
  return(false);
}





//  Encode qualities as 5 bit integers.  Doesn't touch seq.
uint32
gkReadData::gkReadData_encode5bit(uint8 *&UNUSED(chunk), uint8 *UNUSED(qlt), uint32 UNUSED(qltLen)) {
  return(0);
}

bool
gkReadData::gkReadData_decode5bit(uint8 *UNUSED(chunk), uint32 UNUSED(chunkLen), uint8 *UNUSED(qlt), uint32 UNUSED(qltLen)) {
  return(false);
}