Feature suggestion: support BinaryCIF file format

[marcom]

I think it would be good to support the BinaryCIF format. This might perhaps make for a good GSoC project.

"BinaryCIF is a data format for storing text based CIF files using a more efficient binary encoding. "

https://github.com/molstar/BinaryCIF

BinaryCIF is the replacement format for MMTF:

"As of July 2, 2024, RCSB PDB will no longer serve PDB data in the MMTF compression format. Users are strongly encouraged to switch to accessing the data files offered in the compressed BinaryCIF (BCIF) format."

https://www.rcsb.org/news/65a1af31c76ca3abcc925d0c

[jgreener64]

Yes I agree. RIP MMTF.

Hopefully this shouldn't be too hard, as once the file is decoded the existing mmCIF parser can be used.

[BradyAJohnston]

I'm still quite new to Julia, so there are lots of things to improve upon, but I've got a working .bcif parser that I've been using as a way to learn Julia over the past couple of weeks.

It's currently bare-bones, but extracting the _atom_site into a dictionary is ~10x speedup over the parsing of the .cif file, but this will also slow down once it is properly processed into the MolecularStructure.

.cif

julia> @benchmark read("7cgo.cif", MMCIFFormat)
BenchmarkTools.Trial: 5 samples with 1 evaluation per sample.
 Range (min … max):  1.119 s …    1.433 s  ┊ GC (min … max): 17.01% … 34.88%
 Time  (median):     1.298 s               ┊ GC (median):    27.95%
 Time  (mean ± σ):   1.279 s ± 114.868 ms  ┊ GC (mean ± σ):  27.10% ±  6.57%

  █                    █           █ █                     █  
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁█▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█ ▁
  1.12 s         Histogram: frequency by time         1.43 s <

 Memory estimate: 841.24 MiB, allocs estimate: 11335792.

.bcif

julia> @benchmark read("7cgo.bcif", BCIFFormat)
BenchmarkTools.Trial: 32 samples with 1 evaluation per sample.
 Range (min … max):  126.804 ms … 249.072 ms  ┊ GC (min … max):  0.00% … 46.85%
 Time  (median):     146.438 ms               ┊ GC (median):     7.43%
 Time  (mean ± σ):   159.367 ms ±  36.114 ms  ┊ GC (mean ± σ):  16.04% ± 14.23%

        ▂ ▅▂█                                                    
  ▅▁▁▁███████▅█▁▅▁▁▅▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▅▁█▁▅▁▁▅ ▁
  127 ms           Histogram: frequency by time          249 ms <

julia> read("7cgo.bcif", BCIFFormat)
Dict{String, Any} with 21 entries:
  "pdbx_PDB_model_num" => Int32[1, 1, 1, 1, 1, 1, 1, 1, 1, 1  …  1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
  "label_seq_id"       => Int32[1, 1, 1, 1, 1, 1, 1, 1, 2, 2  …  401, 401, 402, 402, 402, 402, 402, 402, 402, 402]
  "Cartn_z"            => [299.079, 300.42, 300.369, 301.365, 301.173, 300.544, 301.444, 301.39, 299.2, 299.004  …  152.812, 151.84, 156…
  "label_entity_id"    => Int32[1, 1, 1, 1, 1, 1, 1, 1, 1, 1  …  8, 8, 8, 8, 8, 8, 8, 8, 8, 8]
  "pdbx_PDB_ins_code"  => ["", "", "", "", "", "", "", "", "", ""  …  "", "", "", "", "", "", "", "", "", ""]
  "pdbx_formal_charge" => ["", "", "", "", "", "", "", "", "", ""  …  "", "", "", "", "", "", "", "", "", ""]
  "B_iso_or_equiv"     => [70.8, 70.8, 70.8, 70.8, 70.8, 70.8, 70.8, 70.8, 66.97, 66.97  …  122.17, 122.17, 127.53, 127.53, 127.53, 127.…
  "auth_comp_id"       => ["MET", "MET", "MET", "MET", "MET", "MET", "MET", "MET", "ILE", "ILE"  …  "ASN", "ASN", "LEU", "LEU", "LEU", "…
  "Cartn_y"            => [310.699, 310.555, 310.065, 309.584, 311.88, 312.989, 314.537, 314.755, 310.188, 309.826  …  341.887, 342.33, …
  "group_PDB"          => ["ATOM", "ATOM", "ATOM", "ATOM", "ATOM", "ATOM", "ATOM", "ATOM", "ATOM", "ATOM"  …  "ATOM", "ATOM", "ATOM", "A…
  "id"                 => Int32[1, 2, 3, 4, 5, 6, 7, 8, 9, 10  …  335713, 335714, 335715, 335716, 335717, 335718, 335719, 335720, 335721…
  "occupancy"          => [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0  …  1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
  "Cartn_x"            => [342.063, 341.513, 340.072, 339.535, 341.588, 340.772, 340.93, 342.702, 339.451, 338.054  …  337.511, 335.537,…
  "type_symbol"        => ["N", "C", "C", "O", "C", "C", "S", "C", "N", "C"  …  "O", "N", "N", "C", "C", "O", "C", "C", "C", "C"]
  "auth_seq_id"        => Int32[1, 1, 1, 1, 1, 1, 1, 1, 2, 2  …  401, 401, 402, 402, 402, 402, 402, 402, 402, 402]
  "auth_asym_id"       => ["A", "A", "A", "A", "A", "A", "A", "A", "A", "A"  …  "EG", "EG", "EG", "EG", "EG", "EG", "EG", "EG", "EG", "E…
  "auth_atom_id"       => ["N", "CA", "C", "O", "CB", "CG", "SD", "CE", "N", "CA"  …  "OD1", "ND2", "N", "CA", "C", "O", "CB", "CG", "CD…
  "label_asym_id"      => ["A", "A", "A", "A", "A", "A", "A", "A", "A", "A"  …  "KH", "KH", "KH", "KH", "KH", "KH", "KH", "KH", "KH", "K…
  "label_atom_id"      => ["N", "CA", "C", "O", "CB", "CG", "SD", "CE", "N", "CA"  …  "OD1", "ND2", "N", "CA", "C", "O", "CB", "CG", "CD…
  "label_alt_id"       => ["", "", "", "", "", "", "", "", "", ""  …  "", "", "", "", "", "", "", "", "", ""]
  "label_comp_id"      => ["MET", "MET", "MET", "MET", "MET", "MET", "MET", "MET", "ILE", "ILE"  …  "ASN", "ASN", "LEU", "LEU", "LEU", "…

As I'm still pretty fresh with Julia, I wanted to ask what the best approach would for turning this dictionary into the proper format?

With the .bcif format, each category of data is decoded all at once. To me it would make sense and be fastest to initialise a structure with empty / default values with the known number of atoms, then as each column of data is decode fill those directly into the structure.

From what I can tell looking at the rest of the BioStructures codebase, this differs to the usual atom-by-atom approach.

[jgreener64]

That sounds great @BradyAJohnston, it's good to have someone like you trying out Julia. If you have a Dict you should be able to call the function at https://github.com/BioJulia/BioStructures.jl/blob/master/src/mmcif.jl#L321 to get the structure (MMCIFDict just wraps a Dict).

This would be very welcome as a PR. There are definitely optimisations that could be made, but we don't currently have an "array of atoms" type internal representation. It's hard for us to know the shape of the hierarchical data structure without going through and reading it.

[BradyAJohnston]

I've been enjoying learning it and just needed a proper project to cut my teeth on. I couldn't just re-use the MolecularStructure function as it was expecting some of the values as strings when read from .cif, but everything is already parsed into numbers when coming from .bcif, but I got it working and returning a structure.

With reading from .bcif we can also multithread the decoding of columns, leading to further speed ups down to ~40 ms. I also tested on the largest structure in the PDB and it was ~300 ms vs 12 s for .cif parsing.

julia> @benchmark  read("7cgo.bcif", BCIFFormat)
BenchmarkTools.Trial: 106 samples with 1 evaluation per sample.
 Range (min … max):  30.472 ms … 95.469 ms  ┊ GC (min … max):  0.00% … 67.02%
 Time  (median):     43.610 ms              ┊ GC (median):    24.83%
 Time  (mean ± σ):   47.235 ms ± 15.639 ms  ┊ GC (mean ± σ):  31.65% ± 17.51%

  ▄      ▃▄ ▅▂ █  ▅                                            
  ██▅▆▁▁▅██▅██▆█▇▅█▇▅▁▃▅▃▃▁▃▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▃▁▃▃▃▁▁▁▃▁▁▁▅▃▅ ▃
  30.5 ms         Histogram: frequency by time        93.4 ms <

 Memory estimate: 177.05 MiB, allocs estimate: 2075542.

If I return it as a MolecularStructure this unfortunately adds a fair amount of processing time, bringing it up to ~600 ms.

julia> @benchmark  read("7cgo.bcif", BCIFFormat)
BenchmarkTools.Trial: 9 samples with 1 evaluation per sample.
 Range (min … max):  537.780 ms … 647.459 ms  ┊ GC (min … max):  5.81% … 20.76%
 Time  (median):     617.231 ms               ┊ GC (median):    16.81%
 Time  (mean ± σ):   604.015 ms ±  37.943 ms  ┊ GC (mean ± σ):  14.48% ±  5.54%

  █                 █ █     █                 █    █  █     █ █  
  █▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁█▁▁▁▁▁█▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁█▁▁▁▁█▁▁█▁▁▁▁▁█▁█ ▁
  538 ms           Histogram: frequency by time          647 ms <

 Memory estimate: 461.82 MiB, allocs estimate: 9341291.

julia> m = read("7cgo.bcif", BCIFFormat)
MolecularStructure with 1 models, 219 chains (0,1,2,3,4,5,6,7,8,9,A,B,C,D,E,F,G,H,I,J,K,L,M,N,O,P,Q,R,S,T,U,V,W,X,a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r,s,t,u,v,w,x,y,z,AA,AB,AC,AD,AE,AF,AG,AH,AI,AJ,AK,AL,AM,AN,AO,AP,AQ,AR,AS,AT,AU,AV,AW,AX,AY,AZ,BA,BB,BC,BD,BE,BF,BG,BH,BI,BJ,BK,BL,BM,BN,BO,BP,BQ,BR,BS,BT,BU,BV,BW,BX,BY,BZ,CA,CB,CC,CD,CE,CF,Ca,Cb,Cc,Cd,Ce,Cf,Cg,Ch,Ci,Cj,Ck,Cl,Cm,Cn,Co,Cp,Cq,Cr,Cs,Ct,Cu,Cv,Cw,Cx,Cy,Cz,DA,DB,DC,DD,DE,DF,DG,DH,DI,DJ,DK,DL,DM,DN,DO,DP,DQ,DR,DS,DT,DU,DV,DW,DX,DY,DZ,Da,Db,Dc,Dd,De,Df,Dg,Dh,Di,Dj,Dk,Dl,Dm,Dn,Do,Dp,Dq,Dr,Ds,Dt,Du,Dv,Dw,EA,EB,EC,ED,EE,EF,EG,Ea,Eb,Ec,Ed,Ee,Ef,Eg,Eh,GA,GB,GC,GD,GE,GF,GG,GH,GI,GJ,GK), 45082 residues, 335722 atoms

I'll clean up the code I have a bit and open a draft PR about it. I'd like to get some feedback and ensure everything is following a "Julia-esque" approach.

I'll admit that the idea for an "array of atoms" makes more sense to me than the current model in BioStructures. I'm a fan of biotite that just uses numpy arrays for each column of data.

I'm not so familiar with optimisations in Julia, but that feels like it would offer performance improvements over the current nesting, but I also might be way off there. As the data already exists, there could maybe be a struct that allows iterating through like the current interface while maintaining arrays underneath? A topic for another thread though.

[jgreener64]

Yes see some discussion at #48. Arrays under the hood would probably be faster. I expect the time taken to create the structure is as much due to unoptimised code though, I wrote it as part of my first Julia package. I could deal with 600 ms for 300k atoms, but I don't deal with large structures much.

You could also just read in the dictionary and use the arrays from there if you don't need the structural information, this is available for CIF files by calling MMCIFDict on the file path.