Parse secondary structure annotations from mmCIF (#65)

Unless `run_dssp` or `run_stride`, this will parse any secondary
structure annotations in the mmCIF file.

---

Co-authored by: Joe Greener <jgreener@hotmail.co.uk>

Author: timholy

docs/src/documentation.md

@@ -286,7 +286,7 @@ The operators currently supported are:
 | Operators                   | Acts on                             |
 | :-------------------------- | :---------------------------------- |
 | `=`, `>`, `<` `>=`, `<=`    | Properties with real values.        |
-| `and`, `or`, `not`          | Selections subsets.                 | 
+| `and`, `or`, `not`          | Selections subsets.                 |
 
 The keywords supported are:
 
@@ -541,13 +541,31 @@ So if you wanted the graph of chain contacts in a protein complex you could give
 
 ## Assigning secondary structure
 
-The secondary structure code of a residue or atom can be accessed after assigning the secondary structure using [DSSP](https://github.com/PDB-REDO/dssp) or [STRIDE](https://webclu.bio.wzw.tum.de/stride).
+Any secondary structure assignment at a residue is accessed via [`sscode`](@ref):
+```julia
+sscode(res)
+sscode(at)
+```
+Or for the whole structure:
+```julia
+sscode.(collectresidues(struc))
+```
+`sscode` may return `'-'`, indicating either that the residue has no recognized secondary structure or that it has not yet been assigned.
+The secondary structure code of a residue can be changed using [`sscode!`](@ref).
+
+mmCIF files with secondary structure [annotations](https://mmcif.wwpdb.org/dictionaries/mmcif_pdbx_v50.dic/Items/_struct_conf_type.id.html) will, by default, parse and assign secondary structure.
+
+Secondary structure can also be assigned using [DSSP](https://github.com/PDB-REDO/dssp) or [STRIDE](https://webclu.bio.wzw.tum.de/stride).
+This functionality is provided through extensions, so you must load the corresponding package library (`DSSP_jll` or `STRIDE_jll`) before this functionality will work.
+
 To assign secondary structure when reading the structure:
 ```julia
 # Assign secondary structure using DSSP
+using DSSP_jll
 read("/path/to/pdb/file.pdb", PDBFormat, run_dssp=true)
 
 # Assign secondary structure using STRIDE
+using STRIDE_jll
 read("/path/to/pdb/file.pdb", PDBFormat, run_stride=true)
 ```
 [`rundssp!`](@ref), [`runstride!`](@ref), [`rundssp`](@ref) and [`runstride`](@ref) can also be used to assign secondary structure to a [`MolecularStructure`](@ref) or [`Model`](@ref):
@@ -557,21 +575,12 @@ runstride!(struc)
 ```
 The assignment process may fail if the structure is too large, since we use an intermediate PDB file where the atom serial cannot exceed 99999 and the chain ID must be a single character.
 To get access to the secondary structure code of a residue or atom as a `Char`:
-```julia
-sscode(res)
-sscode(at)
-```
-Or for the whole structure:
-```julia
-sscode.(collectresidues(struc))
-```
-The secondary structure code of a residue can be changed using [`sscode!`](@ref).
 [`rundssp`](@ref) and [`runstride`](@ref) can also be run directly on structure files:
 ```julia
 rundssp("/path/to/pdb/file.pdb", "out.dssp") # Also works with mmCIF files
 runstride("/path/to/pdb/file.pdb", "out.stride")
 ```
-See the documentation for [DSSP_jll](https://docs.juliahub.com/General/DSSP_jll/stable/autodocs) and [STRIDE_jll](https://docs.juliahub.com/General/STRIDE_jll/stable/autodocs), and also [ProteinSecondaryStructures.jl](https://github.com/m3g/ProteinSecondaryStructures.jl), for other ways to run these programs.
+See the documentation for [DSSP_jll](https://docs.juliahub.com/General/DSSP_jll/stable/autodocs) and [STRIDE_jll](https://docs.juliahub.com/General/STRIDE_jll/stable/autodocs), and also [ProteinSecondaryStructures.jl](https://github.com/BioJulia/ProteinSecondaryStructures.jl), for other ways to run these programs.
 
 ## Downloading PDB files
 

src/mmcif.jl

@@ -41,6 +41,24 @@ const mmciforder = Dict(
     ]
 )
 
+const cif_sscodes = Dict(
+    "HELX_LH_AL_P" => 'H',
+    "HELX_P" => 'H',
+    "STRN" => 'E',
+    "HELX_LH_3T_P" => 'G',
+    "HELX_LH_PI_P" => 'I',
+    "HELX_LH_PP_P" => 'P',
+    "TURN_OT_P" => 'T',
+    "TURN_P" => 'T',
+    "TURN_TY1P_P" => 'T',
+    "TURN_TY2P_P" => 'T',
+    "TURN_TY3P_P" => 'T',
+    "TURN_TY1_P" => 'T',
+    "TURN_TY2_P" => 'T',
+    "TURN_TY3_P" => 'T',
+    "BEND" => 'S',
+)
+
 """
     MMCIFDict(filepath; gzip=false)
     MMCIFDict(io; gzip=false)
@@ -150,17 +168,23 @@ function tokenizecif(f::IO)
     return tokens
 end
 
-# Get tokens from a mmCIF file corresponding to atom site records only
+# Get tokens from a mmCIF file corresponding to atom_site or struct_conf records only
 # This will fail if there is only a single atom record in the file
 #   and it is not in the loop format
 function tokenizecifstructure(f::IO)
     tokens = String[]
     reading = false
     in_keys = true
+    category_groups = ["_atom_site.", "_struct_conf."]
     for line in eachline(f)
-        if (!reading && !startswith(line, "_atom_site.")) || startswith(line, "#")
+        startswith(line, "#") && continue
+        iscat = any(k -> startswith(line, k), category_groups)
+        if !reading && !iscat
             continue
-        elseif startswith(line, "_atom_site.")
+        elseif iscat
+            if !reading
+                push!(tokens, "loop_")
+            end
             reading = true
             push!(tokens, rstrip(line))
         elseif startswith(line, ";")
@@ -175,14 +199,14 @@ function tokenizecifstructure(f::IO)
             end
             push!(tokens, join(token_buffer, "\n"))
         elseif (!in_keys && startswith(line, "_")) || startswith(line, "loop_") ||
-                    startswith(line, "LOOP_")
-            break
+            startswith(line, "LOOP_")
+            reading = false
+            in_keys = true
         else
             in_keys = false
             append!(tokens, splitline(line))
         end
     end
-    length(tokens) > 0 ? pushfirst!(tokens, "loop_") : nothing
     return tokens
 end
 
@@ -322,6 +346,30 @@ function MolecularStructure(mmcif_dict::MMCIFDict;
         fixlists!(struc)
     end
 
+    # Secondary structure assignment
+    if !(run_dssp | run_stride) && haskey(mmcif_dict, "_struct_conf.conf_type_id")
+        for (i, id) in pairs(mmcif_dict["_struct_conf.conf_type_id"])
+            chainid = get(mmcif_dict, "_struct_conf.beg_auth_asym_id", nothing)[i]
+            chainidend = get(mmcif_dict, "_struct_conf.end_auth_asym_id", nothing)[i]
+            if chainid === nothing
+                chainid = mmcif_dict["_struct_conf.beg_label_asym_id"][i]
+                chainidend = mmcif_dict["_struct_conf.end_label_asym_id"][i]
+            end
+            chainidend == chainid || continue   # mismatch in chain id
+            haskey(chains(struc), chainid) || continue
+            chain = struc[chainid]
+            bg = tryparse(Int, mmcif_dict["_struct_conf.beg_auth_seq_id"][i])
+            nd = tryparse(Int, mmcif_dict["_struct_conf.end_auth_seq_id"][i])
+            (bg === nothing || nd === nothing) && continue
+            code = get(cif_sscodes, id, '-')
+            for j in bg:nd
+                r = get(residues(chain), string(j), nothing)
+                isa(r, Residue) || continue
+                r.ss_code = code
+            end
+        end
+    end
+
     if run_dssp && run_stride
         throw(ArgumentError("run_dssp and run_stride cannot both be true"))
     end

test/runtests.jl

@@ -1787,6 +1787,7 @@ end
         @test dic["_pdbx_database_status.recvd_initial_deposition_date"] == ["1991-11-08"]
         @test dic["_audit_author.name"] == ["Mueller, C.W.", "Schulz, G.E."]
         @test length(dic["_atom_site.group_PDB"]) == 3816
+        @test length(dic["_struct_conf.conf_type_id"]) == 18
         dic["_pdbx_database_status.recvd_initial_deposition_date"] = ["changed"]
         @test dic["_pdbx_database_status.recvd_initial_deposition_date"] == ["changed"]
         @test length(keys(dic)) == 610
@@ -1942,10 +1943,14 @@ end
 
     open(testfilepath("mmCIF", "1AKE.cif")) do f
         tokens = tokenizecifstructure(f)
-        @test length(tokens) == 80158
-        @test tokens[1] == "loop_"
-        @test tokens[10] == "_atom_site.label_seq_id"
-        @test tokens[80089] == "77.69"
+        loopidxs = findall(==("loop_"), tokens)
+        @test length(loopidxs) == 2
+        idxsc = loopidxs[1]
+        @test tokens[idxsc + 1] == "_struct_conf.conf_type_id"
+        idxas = loopidxs[2]
+        @test length(tokens) - idxas + 1 == 80158
+        @test tokens[idxas + 9] == "_atom_site.label_seq_id"
+        @test tokens[idxas + 80088] == "77.69"
     end
 
     # Test parsing empty file
@@ -1996,6 +2001,8 @@ end
     res = collect(chs[1])
     @test length(res) == 456
     @test resid(res[20]) == "20"
+    @test sscode(res[11]) == '-'
+    @test sscode(res[12]) == 'H'
     ats = collect(res[20])
     @test length(ats) == 8
     @test atomname.(ats) == ["N", "CA", "C", "O", "CB", "CG1", "CG2", "CD1"]