convert long string to integer list

Question

at first, i would like to thank all the people here. i'm working with mathematica for half a year now and always found great advice and tricks.

now, i came up with a problem i couldn't find a solution for. from a file, i import a string with around 23 mio. characters which is reasonably fast (it acutually is a dna sequence). i then want to split this string into its characters and convert each of them to an integer number.

my first, naive approach was the following:

MaxMemoryUsed[]/1024^2//N
AbsoluteTiming[list=Characters[a]/.{"A"->1,"T"->2,"G"->3,"C"->4};]
ByteCount[list]/1024^2//N
MaxMemoryUsed[]/1024^2//N

which creates the following output:

134.286
{13.301206,Null}
702.261
1976.42

as you can see, the timing is quite good but the memory usage during the process is horrible. my next approach was reading from a stream and convert each character to an integer instantly:

MaxMemoryUsed[]/1024^2//N
string=StringReplace[string,{"A"->"1","T"->"2","G"->"3","C"->"4"}];
str=StringToStream[string];
AbsoluteTiming[list=Table[ToExpression[Read[str,Character]],{i,1,StringLength[string]}];]
Close[str];
ByteCount[list]/1024^2//N
MaxMemoryUsed[]/1024^2//N

which results in

134.281
{176.502736,Null}
702.252
813.208

as yo can see, this approach needs way less memory but takes considerably longer. my question now is, if there is another memory efficient way to do this which is faster than my second try?

thanks a lot!

Answer 1

string = StringJoin @@ RandomChoice[{"A", "T", "G", "C"}, 30000000];

This seems to take little memory

Boole[# === 65] + 2 Boole[# === 84] + Boole[# === 71] 3 + 
   Boole[# === 67] 4 & /@ ToCharacterCode[string]

but twice as long as

ToCharacterCode[
  StringReplace[
   string, {"A" -> "1", "T" -> "2", "G" -> "3", "C" -> "4"}]] - 48

Compiling the first solution makes it faster than the second

Compile[{{x, _Integer}}, \[Piecewise] {
    {1, x === 65},
    {2, x === 84},
    {3, x === 71},
    {4, x === 67}
   }, RuntimeAttributes -> Listable]@ToCharacterCode[string]

Answer 2

OK, we don't need command-line tools for this. We can do better, faster and with less memory.

Lets create a test file with 30 million letters:

Export["tmp/sequence.txt", StringJoin@RandomChoice[{"A", "C", "G", "T"}, 30000000]]

The trick is to make the conversion parallel on all data and to keep packet-arrays. First we define a parallel, compiled function which does the conversion for one character:

repl = Compile[{{n, _Integer, 0}},
  If[n == 65, 1,
   If[n == 67, 2,
    If[n == 71, 3,
     If[n == 84, 4, -1]
     ]
    ]
   ], CompilationTarget -> "C", RuntimeAttributes -> {Listable}, 
  Parallelization -> True
  ]

Yes, those Ifs look pretty awful, but they get compiled down while Switch or Which doesn't. See, that we don't read characters, we read the ASCII codes instead which is a nicely packed integer array:

MaxMemoryUsed[]/1024^2 // N
AbsoluteTiming[dna = repl[Developer`ToPackedArray[ReadList["tmp/sequence.txt", Byte]]];]
MaxMemoryUsed[]/1024^2 // N

Here this takes 2 seconds and less than 500 MB of memory.

Timings

I timed all methods on my machine. I really don't know how image_doctor got his method so fast. It was the slowest here (native Linux, SSD, Intel i7 extreme, 32 GB Ram). I didn't preload the data. If a method used a string, I used Import.

Currently fastest method

I thought, that it's Rojos Piecewise construct was faster, but as it turns out this was cosmic radiation since the ugly nested Ifs and the Piecewise solution produce exact the same compiled code. The full test code consists of a sequence of only 3 million digits

Export["sequence.txt", StringJoin@RandomChoice[{"A", "C", "G", "T"}, 3000000]];

which is called 100 times.

halirutan

f = Compile[{{n, _Integer, 0}}, 
  If[n == 65, 1, If[n == 67, 2, If[n == 71, 3, If[n == 84, 4, -1]]]], 
  CompilationTarget -> "C", RuntimeAttributes -> {Listable}, 
  Parallelization -> True];

Mean@Table[First@AbsoluteTiming[
    f[ToCharacterCode[Import["tmp/sequence.txt"]]];], {100}]

(* Out[2]= 0.181717 *)

Rojo compiled to "C" and parallel

f = Compile[{{x, _Integer}}, Piecewise[{{1, x === 65}, 
           {2, x === 84}, {3, x === 71}, {4, x === 67}}], 
       CompilationTarget -> "C", RuntimeAttributes -> {Listable},
   Parallelization -> True];

Mean@Table[First@AbsoluteTiming[f[
      ToCharacterCode[Import["tmp/sequence.txt"]]];], {100}]

(* Out[3]= 0.183226 *)

Note here, how mystical it is, that Piecewise (and Ifs with then branch (!)) get compiled correctly to "C", while Switch and Which stays in a MainEvaluate call.

Answer 3

I read 30 million characters stored in around 10,000 lines in a text file into a list of integers like this.

dna = Flatten[IntegerDigits /@ ReadList["!cat /tmp/dna.tmp | tr [ACGT] [1234]"]];
//AbsoluteTiming

{2.766828, Null}

Memory usage was 0.93GB for 29.6 mega characters.

The 0.93GB is effectively the space required for the 30 million integers. I guess this storage amount could be compressed if required as the underlying alphabet only has 4 symbols.

Answer 4

On my PC, using a lookup table is slightly faster than @halirutan's If or Piecewise method.

Export["tmp/sequence.txt", StringJoin@RandomChoice[{"A", "C", "G", "T"}, 30000000]]

lookup = ConstantArray[0, 256];
lookup[[ToCharacterCode["ACGT"]]] = Range[4];

f = Compile[{{x, _Integer}, {lut, _Integer, 1}}, lut[[x]],
   CompilationTarget -> "C", RuntimeAttributes -> {Listable}, 
   Parallelization -> True];

AbsoluteTiming[f[ToCharacterCode[Import["tmp/sequence.txt"]], lookup];]

Takes 3.1198962s, compared to 3.2079074 with the Piecewise code in @halirutan's answer.

If you don't care what the function returns for characters other than A,C,G,T, this is even faster:

f = Compile[{{x, _Integer, 1}},
   UnitStep[x - 65] + UnitStep[x - 84] + UnitStep[x - 71] + UnitStep[x - 67]];
AbsoluteTiming[f[ToCharacterCode[Import["sequence.txt"]]];]

Takes 2.9723774s

(Probably because UnitStep can somehow use SIMD instructions and If and Part cannot)

Answer 5

Though not quite as fast as Rojo's second method this is a little shorter and uses less memory on my machine:

ToCharacterCode @ string /. {65 -> 1, 84 -> 2, 71 -> 3, 67 -> 4}

You could also read your data as bytes to bypass ToCharacterCode like this:

BinaryReadList["gen.txt"] /. {65 -> 1, 84 -> 2, 71 -> 3, 67 -> 4}

However, his new compiled form is as fast and uses much less memory.

Comparison

First generating test data:

string = StringJoin @@ RandomChoice[{"A", "T", "G", "C"}, 30000000];

DumpSave["bigGen.txt", string];

Quit[], then:

ToCharacterCode[StringReplace[
     ReadList["bigGen.txt", String],
     {"A" -> "1", "T" -> "2", "G" -> "3", "C" -> "4"}]] - 48 // 
  Timing // First

MaxMemoryUsed[]

2.075

850227296

Quit[] again, and:

BinaryReadList["bigGen.txt"] /. {65 -> 1, 84 -> 2, 71 -> 3, 67 -> 4} // Timing // First

MaxMemoryUsed[] 

2.387

614962712

Now with Rojo's compiled piecewise form:

f = Compile[{{x, _Integer}},
   \[Piecewise]{{1, x === 65}, {2, x === 84}, {3, x === 71}, {4, x === 67}}];

f /@ BinaryReadList["bigGen.txt"] // Timing // First

MaxMemoryUsed[] 

2.262

254977528

This should be even better in Version 8 with the enhanced Compile.