I want to reduce memory storage of chars using Bitwise Operations,
for Example
input: {"ACGT"}
output: {"0xE4"}
where we represtnt the number in binary then into Hexadecimal
if A=00, C=01, G=10, T=11
so ACGT = 0xE4 = 11100100b
I cant Figure the whole way so here is what I did So Far
enum NucleicAc {
A = 0,
C = 1,
G = 2,
T = 3,
} ;
struct _DNA {
char * seq ;
};
DNAString DSCreate(char * mseq) {
DNAString dna = malloc(sizeof(struct _DNA));
if (!dna){
exit(-1);
}
const int length = strlen(mseq);
// left shift will create the base , in the case of ACGT --> 0000,0000
int base = 0 << length * 2;
//loop and set bits
int i = 0 ;
int k = 0 ; // the counter where i want to modify the current bit
//simple looping gor the string
for ( ; i < length ; i ) {
switch (*(mseq i)) {
case 'A': // 00
k ;
break;
case 'C': // 0 1
modifyBit(&base, k, 1);
k ;
break;
case 'G': //10
k ;
modifyBit(&base,k , 1);
break;
case 'T': // 11
modifyBit(&base, k, 1);
k ;
modifyBit(&base, k,1);
break;
default:
break;
} //end of switch
k ;
}//end of for
char * generatedSeq ;
//convert the base to hex ??
return dna;
}
void bin(unsigned n){
unsigned i;
for (i = 1 << 7; i > 0; i = i / 2){
(n & i) ? printf("1") : printf("0");
}
}
and if we print the base, the value is 11100100b as expected,
how to store the hexadecimal representation as String to the char *mseq in the struct ? any direction or exact solutions or is there any better approach?
and also later i want to get the letter using only the index for Example
DSGet(dsStruct , '0')--> will return 'A' hence dsStruct contains the "ACGT" before Encode?
CodePudding user response：
Assuming you really do want to encode your dna string into a hex string and that you want to read the input string from left to right, but output hex chars right to left, here's a simple, but slightly slow implementation.
First, your DNAString needs to keep track whether there are really an even or odd number of acid sequences in the list. This will make additional appendages easier.
struct DNAString
{
char* seq;
bool odd; // if odd bit is set, then the front char is already allocated and holds one acid
};
And now let's introduce a little helper function to convert ACGT into 0,1,2,3.
char acid_to_value(char c)
{
switch (c)
{
case 'A': return 0;
case 'C': return 1;
case 'G': return 2;
case 'T': return 3;
}
// ASSERT(false)
return 0;
}
Then the core implementation is to keep "prepending" new hexchars onto the string your are building. If the string is already of an odd length, the code will just "fixup" the front character by converting it from hex to integer, then shifting the new acid value into it, then converting it back to a hex char
extern char fixup(char previous, char acid);
{
char hexchars[16] = { '0','1','2','3','4','5','6','7','8','9','A','B','C','D','E','F' };
char tmp[2] = { previous, '\0' };
unsigned long asnumber = strtol(tmp, nullptr, 16);
asnumber = asnumber & 0x3; // last two bits
asnumber = asnumber | (acid_to_value(acid) << 2);
return hexchars[asnumber];
}
void prepend_nucleic_acid_to_hexstring(struct DNAString* dnastring, char acid)
{
if (dnastring->odd)
{
// find the first char in the string and fix it up hexwise
dnastring->seq[0] = fixup(dnastring->seq[0], acid);
dnastring->odd = false;
}
else
{
size_t currentlength = dnastring->seq ? strlen(dnastring->seq) : 0;
const char* currentstring = dnastring->seq ? dnastring->seq : "";
char* newseq = (char*)calloc(currentlength 2, sizeof(char)); // 1 for new char and 1 for null char
newseq[0] = acid_to_value(acid) '0'; // prepend the next hex char
strcpy(newseq 1, currentstring); // copy the old string into the new string space
free(dnastring->seq);
dnastring->seq = newseq;
dnastring->odd = true;
}
}
Then your DNACreate function is real simple:
struct DNAString DSCreate(const char* mseq)
{
DNAString dnastring = { 0 };
while (*mseq)
{
prepend_nucleic_acid_to_hexstring(&dnastring, *mseq);
mseq ;
}
return dnastring;
}
I don't claim this approach to be efficient since he literally keeps reallocating memory for each char. But it does enable you to have flexability to invoke the prepend function later for additional sequencing.
And then to test it:
int main()
{
struct DNAString dnastring = DSCreate("ACGT");
printf("0x%s\n", dnastring.seq);
return 0;
}
CodePudding user response：
This should do what you want. I thought the compression was a pretty cool idea so I wrote this real quick. As mentioned by @kaylum, hex encoding is just a way to read the underlying data in memory, which is always just bits. So, you only need to worry about that on print statements.
Let me know if this works or you have any questions about what I did.
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
typedef struct {
unsigned char *bits;
unsigned long length; // use this to store the number of letters encoded, for safety
} DNA;
typedef enum {
A = 0,
C = 1,
G = 2,
T = 3
} NucleicAc;
This returns the base at a given index with some bounds checking
char base_at_index(DNA *dna, unsigned long index) {
if (index >= dna->length) {
fputs("error: base_at_index: index out of range", stderr);
exit(EXIT_FAILURE);
}
// offset is index / 4, this gives us the correct byte
// shift amount is index % 4 to give us the correct 2 bits within the byte.
// This must then be multiplied by 2 because
// each base takes 2 bits to encode
// then we have to bitwise-and this value with
// 3 (0000 0011 in binary) to retrieve the bits we want.
// so, the formula we need is
// (dna->bits[index / 4] >> (2 * (index % 4))) & 3
switch((dna->bits[index / 4] >> 2 * (index % 4)) & 3) {
case A: return 'A';
case C: return 'C';
case G: return 'G';
case T: return 'T';
default:
fputs("error: base_at_index: invalid encoding", stderr);
exit(EXIT_FAILURE);
}
}
This encodes a string of bases to bytes
/* you can fit four 2-bit DNA codes in each byte (unsigned char).
len is the maximum number of characters to read. result must be at least len bytes long
*/
void encode_dna(unsigned char *result, char *sequence, unsigned long len) {
// keep track of what byte we are on in the result
unsigned result_index = 0;
// our shift for writing to the correct position in the byte
unsigned shift = 0;
// first clear result or else bitwise operations will produce errors
// this could be removed if you were certain result parameter was zero-filled
memset(result, 0, len);
// iterate through characters of the sequence
while(*sequence) {
switch (*sequence) {
// do nothing for 'A' since it is just zero
case 'A': break;
case 'C':
// we are doing a bitwise or with the current byte
// and C (1) shifted to the appropriate position within
// the byte, and then assigning the byte with the result
result[result_index] |= C << shift;
break;
case 'G':
result[result_index] |= G << shift;
break;
case 'T':
result[result_index] |= T << shift;
break;
default:
fputs("error: encode_dna: invalid base pair", stderr);
exit(EXIT_FAILURE);
}
// increase shift amount by 2 to the next 2-bit slot in the byte
shift = 2;
// on every 4th iteration, reset our shift to zero since the byte is now full
// and move to the next byte in our result buffer
if (shift == 8) {
shift = 0;
result_index ;
}
// advance sequence to next nucleotide character
sequence ;
}
}
And here's a test
int main(int argc, char **argv) {
// allocate some storage for encoded DNA
unsigned char encoded_dna[32];
const unsigned long sample_length = 15;
// encode the given sample sequence
encode_dna(encoded_dna, "ACGTAGTCGTCATAG", sample_length);
// hh here means half of half word, which is a byte
// capital X for capitalized hex output
// here we print some bytes
printf("0x%hhX\n", encoded_dna[0]); // should output 0xE4
printf("0x%hhX\n", encoded_dna[1]); // should be 0x78
printf("0x%hhX\n", encoded_dna[2]); // should be 0x1E
printf("0x%hhX\n", encoded_dna[3]); // should be 0x23
DNA test_dna; // allocate a sample DNA structure
test_dna.bits = encoded_dna;
test_dna.length = sample_length; // length of the sample sequence above
// test some indices and see if the results are correct
printf("test_dna index 4: %c\n", base_at_index(&test_dna, 4));
printf("test_dna index 7: %c\n", base_at_index(&test_dna, 7));
printf("test_dna index 12: %c\n", base_at_index(&test_dna, 12));
return 0;
}
Output:
0xE4
0x78
0x1E
0x23
test_dna index 4: A
test_dna index 7: C
test_dna index 12: T
CodePudding user response：
There are several ways you can approach encoding the sequence. Your enum is fine, but for your encoded sequence, a struct that captures the bytes as unsigned char, the original sequence length and the encoded size in bytes will allow an easy decoding. You will get 4-to-1 compression (plus 1 byte if you sequence isn't evenly divisible by 4). The enum and struct could be:
enum { A, C, G, T };
typedef struct {
unsigned char *seq;
size_t len, size;
} encoded;
To map characters in your string to the encoded values a simple function that returns the enum value matching the character is all you need (don't forget to handle any errors)
/* convert character to encoded value */
unsigned char getencval (const char c)
{
if (c == 'A')
return A;
else if (c == 'C')
return C;
else if (c == 'G')
return G;
else if (c == 'T')
return T;
/* exit on anything other than A, C, G, T */
fprintf (stderr, "error: invalid sequence character '%c'\n", c);
exit (EXIT_FAILURE);
}
To encode the original sequence, you will populate the encoded struct with the original length (len) and number of bytes (size) needed to hold the encoded string. Don't forget that any 1-character of the next 4-characters will require another byte of storage. You can use a simple add-and-divide to account for any partial 4-character ending portion of the sequence, e.g.
/* encode sequence of characters as 2-bit pairs (4-characters per-byte)
* returns encoded struct with allocated .seq member, on failure the .seq
* member is NULL. User is resposible for freeing .seq member when done.
*/
encoded encode_seq (const char *seq)
{
size_t len = strlen(seq),
size = (len 3) / 4; /* integer division intentional */
encoded enc = { .seq = calloc (1, size), /* encoded sequence struct */
.len = len,
.size = size };
if (!enc.seq) { /* validate allication */
perror ("calloc-enc.seq");
return enc;
}
/* loop over each char (i) with byte index (ndx)
* shifting each 2-bit pair by (shift * 2) amount.
*/
for (int i = 0, ndx = 0, shift = 4; seq[i] && seq[i] != '\n'; i ) {
if (!shift--) /* decrement shift, reset if 0 */
shift = 3;
if (i && i % 4 == 0) /* after each 4th char, increment ndx */
ndx = 1;
/* shift each encoded value (multiply by 2 for shift of 6, 4, 2, 0) */
enc.seq[ndx] |= getencval (seq[i]) << shift * 2;
}
return enc; /* return encoded struct with allocated .seq member */
}
To get the original sequence back from your encoded struct, the use of a lookup table (shown below with the full code) makes it a breeze. You simply loop over all stored byte values appending the corresponding strings from the lookup table until the final byte. For the final byte, you need to determine if it is a partial string and, if so, how many characters remain to copy. (that's why you store the original sequence length in your struct). Then simply use strncat to append that many characters from the final byte, e.g.
/* decodes encoded sequence. Allocates storage for decoded sequence
* and loops over each encoded byte using lookup-table to obtain
* original 4-character string from byte value. User is responsible
* for freeing returned string when done. Returns NULL on allocation
* failure.
*/
char *decode_seq (encoded *eseq)
{
char *seq = malloc (eseq->len 1); /* allocate storage for sequence */
size_t i = 0, offset = 0, remain;
if (!seq) { /* validate allocation */
perror ("malloc-seq");
return NULL;
}
/* loop appending strings from lookup table for all but last byte */
for (; i < eseq->size - 1; i ) {
memcpy (seq offset, lookup[eseq->seq[i]], 4);
offset = 4; /* increment offset by 4 */
}
/* determine the number of characters in last byte */
remain = eseq->len - (eseq->size - 1) * 4;
memcpy (seq offset, lookup[eseq->seq[i]], remain);
seq[offset remain] = 0; /* nul-terminate seq */
return seq; /* return allocated sequence */
}
Adding the lookup table and putting all the pieces together, one way to approach this problem is:
(edit: lookup table reordered to match your byte-value encoding, optimized decode_seq() to not scan for end-of-string on copy)
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
enum { A, C, G, T };
typedef struct {
unsigned char *seq;
size_t len, size;
} encoded;
const char lookup[][5] = {
"AAAA","CAAA","GAAA","TAAA","ACAA","CCAA","GCAA","TCAA",
"AGAA","CGAA","GGAA","TGAA","ATAA","CTAA","GTAA","TTAA",
"AACA","CACA","GACA","TACA","ACCA","CCCA","GCCA","TCCA",
"AGCA","CGCA","GGCA","TGCA","ATCA","CTCA","GTCA","TTCA",
"AAGA","CAGA","GAGA","TAGA","ACGA","CCGA","GCGA","TCGA",
"AGGA","CGGA","GGGA","TGGA","ATGA","CTGA","GTGA","TTGA",
"AATA","CATA","GATA","TATA","ACTA","CCTA","GCTA","TCTA",
"AGTA","CGTA","GGTA","TGTA","ATTA","CTTA","GTTA","TTTA",
"AAAC","CAAC","GAAC","TAAC","ACAC","CCAC","GCAC","TCAC",
"AGAC","CGAC","GGAC","TGAC","ATAC","CTAC","GTAC","TTAC",
"AACC","CACC","GACC","TACC","ACCC","CCCC","GCCC","TCCC",
"AGCC","CGCC","GGCC","TGCC","ATCC","CTCC","GTCC","TTCC",
"AAGC","CAGC","GAGC","TAGC","ACGC","CCGC","GCGC","TCGC",
"AGGC","CGGC","GGGC","TGGC","ATGC","CTGC","GTGC","TTGC",
"AATC","CATC","GATC","TATC","ACTC","CCTC","GCTC","TCTC",
"AGTC","CGTC","GGTC","TGTC","ATTC","CTTC","GTTC","TTTC",
"AAAG","CAAG","GAAG","TAAG","ACAG","CCAG","GCAG","TCAG",
"AGAG","CGAG","GGAG","TGAG","ATAG","CTAG","GTAG","TTAG",
"AACG","CACG","GACG","TACG","ACCG","CCCG","GCCG","TCCG",
"AGCG","CGCG","GGCG","TGCG","ATCG","CTCG","GTCG","TTCG",
"AAGG","CAGG","GAGG","TAGG","ACGG","CCGG","GCGG","TCGG",
"AGGG","CGGG","GGGG","TGGG","ATGG","CTGG","GTGG","TTGG",
"AATG","CATG","GATG","TATG","ACTG","CCTG","GCTG","TCTG",
"AGTG","CGTG","GGTG","TGTG","ATTG","CTTG","GTTG","TTTG",
"AAAT","CAAT","GAAT","TAAT","ACAT","CCAT","GCAT","TCAT",
"AGAT","CGAT","GGAT","TGAT","ATAT","CTAT","GTAT","TTAT",
"AACT","CACT","GACT","TACT","ACCT","CCCT","GCCT","TCCT",
"AGCT","CGCT","GGCT","TGCT","ATCT","CTCT","GTCT","TTCT",
"AAGT","CAGT","GAGT","TAGT","ACGT","CCGT","GCGT","TCGT",
"AGGT","CGGT","GGGT","TGGT","ATGT","CTGT","GTGT","TTGT",
"AATT","CATT","GATT","TATT","ACTT","CCTT","GCTT","TCTT",
"AGTT","CGTT","GGTT","TGTT","ATTT","CTTT","GTTT","TTTT"};
/* convert character to encoded value */
unsigned char getencval (const char c)
{
if (c == 'A')
return A;
else if (c == 'C')
return C;
else if (c == 'G')
return G;
else if (c == 'T')
return T;
/* exit on anything other than A, C, G, T */
fprintf (stderr, "error: invalid sequence character '%c'\n", c);
exit (EXIT_FAILURE);
}
/* encode sequence of characters as 2-bit pairs (4-characters per-byte)
* returns encoded struct with allocated .seq member, on failure the .seq
* member is NULL. User is resposible for freeing .seq member when done.
*/
encoded encode_seq (const char *seq)
{
size_t len = strlen(seq),
size = (len 3) / 4; /* integer division intentional */
encoded enc = { .seq = calloc (1, size), /* encoded sequence struct */
.len = len,
.size = size };
if (!enc.seq) { /* validate allication */
perror ("calloc-enc.seq");
return enc;
}
/* loop over each char (i) with byte index (ndx)
* shifting each 2-bit pair by (shift * 2) amount.
*/
for (int i = 0, ndx = 0, shift = 0; seq[i] && seq[i] != '\n'; i , shift ) {
if (shift == 4) /* reset to 0 */
shift = 0;
if (i && i % 4 == 0) /* after each 4th char, increment ndx */
ndx = 1;
/* shift each encoded value (multiply by 2 for shift of 0, 2, 4, 6) */
enc.seq[ndx] |= getencval (seq[i]) << shift * 2;
}
return enc; /* return encoded struct with allocated .seq member */
}
/* decodes encoded sequence. Allocates storage for decoded sequence
* and loops over each encoded byte using lookup-table to obtain
* original 4-character string from byte value. User is responsible
* for freeing returned string when done. Returns NULL on allocation
* failure.
*/
char *decode_seq (encoded *eseq)
{
char *seq = malloc (eseq->len 1); /* allocate storage for sequence */
size_t i = 0, offset = 0, remain;
if (!seq) { /* validate allocation */
perror ("malloc-seq");
return NULL;
}
/* loop appending strings from lookup table for all but last byte */
for (; i < eseq->size - 1; i ) {
memcpy (seq offset, lookup[eseq->seq[i]], 4);
offset = 4; /* increment offset by 4 */
}
/* determine the number of characters in last byte */
remain = eseq->len - (eseq->size - 1) * 4;
memcpy (seq offset, lookup[eseq->seq[i]], remain);
seq[offset remain] = 0; /* nul-terminate seq */
return seq; /* return allocated sequence */
}
/* short example program that takes string to encode as 1st argument
* using "ACGT" if no argument is provided by default
*/
int main (int argc, char **argv) {
char *seq = NULL;
encoded enc = encode_seq(argc > 1 ? argv[1] : "ACGT");
if (!enc.seq) /* validate encoded allocation */
return 1;
/* output original string, length and encoded size */
printf ("encoded str : %s\nencoded len : %zu\nencoded size : %zu\n",
argc > 1 ? argv[1] : "ACGT", enc.len, enc.size);
/* loop outputting byte-values of encoded string */
fputs ("encoded seq :", stdout);
for (size_t i = 0; i < enc.size; i )
printf (" 0xx", enc.seq[i]);
putchar ('\n');
seq = decode_seq (&enc); /* decode seq from byte values */
printf ("decoded seq : %s\n", seq); /* output decoded string */
free (seq); /* don't forget to free what you allocated */
free (enc.seq);
}
In most cases a lookup-table provides a great deal of efficiency advantage compared to computing and building each 4-character string during decoding. This is enhanced by the lookup table staying resident in cache for most cases.
The length of the DNA sequence you can encode and decode is limited only by the amount of virtual memory you have available.
Example Use/Output
The program takes the sequence to encode and decode as the first argument (default "ACGT"). So the default output is:
$ ./bin/dnaencodedecode
encoded str : ACGT
encoded len : 4
encoded size : 1
encoded seq : 0xe4
decoded seq : ACGT
4-byte encoded in 1-byte. Note the byte value of 0x1b and not 0xe4 due to the table ordering.
A longer example:
./bin/dnaencodedecode ACGTGGGTCAGACTTA
encoded str : ACGTGGGTCAGACTTA
encoded len : 16
encoded size : 4
encoded seq : 0xe4 0xea 0x21 0x3d
decoded seq : ACGTGGGTCAGACTTA
16-character encoded in 4-bytes.
Finally, what of a sequence that isn't divisible by 4 so you have a partial number of characters in the last encoded byte? That is handled as well, e.g.
$ ./bin/dnaencodedecode ACGTGGGTCAGACTTAG
encoded str : ACGTGGGTCAGACTTAG
encoded len : 17
encoded size : 5
encoded seq : 0xe4 0xea 0x21 0x3d 0x02
decoded seq : ACGTGGGTCAGACTTAG
17 characters encoded in 5-bytes. (not the pure 4-to-1 compression, but as the sequence size increases, the significance of any partial group of characters in the last byte becomes negligible)
As far a perfomance, for a sequence of 100,000 characters and output of the the byte values and strings replaced with a simple loop that compares the decoded seq to the original argv[1] it only takes a few thousandth of a second (on an old i7 Gen2 laptop with SSD) to encode and decode and validate, e.g.
$
time ./bin/dnaencodedecodet2big $(< dat/dnaseq100k.txt)
encoded len : 100000
encoded size : 25000
all tests passed
real 0m0.014s
user 0m0.012s
sys 0m0.003s
There are a lot of ways to do this, but given your description, this was what came to my mind that you were trying to accomplish. There is a lot here, so take your time going through it.
Look things over (the code is commented), and let me know if you have further questions. Just drop a comment below.