Convert a bunch of individual STL files to a large mesh collection sff file.

stl2sff.py

## Convert all mitochondria STL files in current directory to SFF/HFF files
## Usage: python stl2sff.py FOLDERNAME
## Assumes there are multiple mitochondria per tomogram and up to one ER per tomogram.
## Performs both per-segment and global annotations based on assumption of MitoGFP MEF cells treated with thapsigargin (sometimes)
## Outputs HFF format sff files.
## Author: Benjamin Barad <[email protected]> with guidance from Paul Korir

import os
import sys
import glob
import numpy as np


## I could definitely clean up these imports! But these are the ones I ended up using.
from sfftk.readers import stlreader
from sfftk.formats.stl import STLSegment
from sfftk.notes.modify import SimpleNote, ExternalReference
import sfftkrw as sff


os.chdir(sys.argv[1])
stlfiles = glob.glob("*.stl")
basenames = set([i.split("_")[0] for i in stlfiles])

## External references are things I do not totally understand - why one API for the segment references vs another for the global external references?
## It might be possible to only do one.
## External references trigger an http request when declared, so I do them here and then reference them for each entry to speed things up.
immref = ExternalReference(resource='ncit', url="http://purl.obolibrary.org/obo/NCIT_C13331", accession="NCIT_C13331")
ommref = ExternalReference(resource='ncit', url="http://purl.obolibrary.org/obo/NCIT_C13330", accession="NCIT_C13330")
erref = ExternalReference(resource='ncit', url="http://purl.obolibrary.org/obo/NCIT_C13230", accession="NCIT_C13230")
mefref = sff.SFFExternalReference(resource='efo', url="http://www.ebi.ac.uk/efo/EFO_0004040", accession="EFO_0004040")
thapsref = sff.SFFExternalReference(resource='chebi', url="https://www.ebi.ac.uk/chebi/chebiOntology.do?chebiId=CHEBI:9516", accession="CHEBI_9516")



for basename in basenames:
    ## All the files we might need for each tomogram!
    imms = [i for i in stlfiles if basename+"_IMM" in i]
    omms = [i for i in stlfiles if basename+"_OMM" in i]
    ers = [i for i in stlfiles if basename+"_ER" in i]
    ## Condense down to only 1 or 0 ER - for simpler coding we used multiple copies of ER (1 per mito) during mito-oriented quantifications.
    if len(ers) > 0:
        er_bool = True
        er = ers[0]
    else:
        er_bool = False
    
    ## This is just how we reference treated vs untreated, but it is a useful way to grab details for the description.
    if basename[0]=="T":
        treatment = "Thapsigargin-treated"
    else:
        treatment = "Vehicle-treated"
    
    ## Straightforward to include name and details during instantiation. 
    ## Also need to instantiate as a mesh_list descriptor - for now, you only get 1 descriptor per sff so no obvious way to mix meshes with voxels with particle positions.
    seg = sff.SFFSegmentation(name="Seg {}".format(basename), details=f"Segmentation of mitochondria and potentially ER in tomogram {basename} of a {treatment} mouse embryonic fibroblast expressing mitochondria-targeted GFP.", primary_descriptor="mesh_list")
    ## Global references for the SFF
    seg.global_external_references = sff.SFFGlobalExternalReferenceList()
    seg.global_external_references.append(mefref)
    if treatment == "Thapsigargin-treated":
        seg.global_external_references.append(thapsref)
        
        
    ## Step two of annotation - all the software we used! Since each one is unique I didn't bother with a loop and just did these by hand.
    seg.software_list = sff.SFFSoftwareList()

    seg.software_list.append(
        sff.SFFSoftware(
            name='Tomosegmemtv',
            version='Apr2020',
            processing_details='We used tomosegmemtv to enhance the membranes of tomograms for further segmentations.'
        )
    )   
    seg.software_list.append(
        sff.SFFSoftware(
            name='AMIRA',
            version='2020.1',
            processing_details='Tomosegmemtv "Saliency" output was loaded into amira and enhanced membranes were separated and cleaned up with the segmentation panel and the magic wand and paintbrush tools. They were then exported as 3D segmentation volumes.'
        )
    )       
    seg.software_list.append(
        sff.SFFSoftware(
            name='Surface Morphometrics Pipeline',
            version='v1.0beta',
            processing_details='Voxel segmentations were used to generate high quality meshes using the screened poisson method in the surface morphometrics pipeline. These were then used for quantifications using the built in tools in the pipeline'
        )
    )   
    seg.software_list.append(
        sff.SFFSoftware(
            name="sfftk",
            version="0.8.3",
            processing_details="Used the python bindings to perform format conversions to sff"
        )
    )
    
    ## Done with global annotation - time to load segments!
    seg.segment_list = sff.SFFSegmentList()
  
    ## Each of these 3 loops ends up looking basically identical. 
    ## Load data, add annotations, add color. Just for IMM, then OMM, then ER.
    for immname in imms:
        mito_number = int(immname.split("_")[3])
        print(immname)
        ### This is the best way I found to load in stl data - there may be a better API eventually.
        stl_seg = stlreader.get_data(immname)
        stl_segment = STLSegment(*stl_seg[0])
        sff_segment = stl_segment.convert()
        ### Simple note interface is what I found best for annotating segment/meshes. 
        ## This was advised by Paul Korir and I found it much easier than directly modifying the biological annotation.
        note=SimpleNote()
        note.name=f"{basename}_mito_{mito_number}_IMM"
        note.description=f"Screened poisson generated mesh depicting IMM (inner mitochondrial membrane) of mitochondria #{mito_number} of tomogram {basename}"
        note.external_references = [immref]
        sff_segment = note.add_to_segment(sff_segment)
        ## We use HEX for our standard colors in the lab so I converted by hand here.
        sff_segment.colour= sff.SFFRGBA(
            red=237/255,
            green=141/255,
            blue=52/255,
            alpha=1
        )
        seg.segment_list.append(sff_segment)

    
    for ommname in omms:
        mito_number = int(ommname.split("_")[3])
        print(ommname)
        stl_seg = stlreader.get_data(ommname)
        stl_segment = STLSegment(*stl_seg[0])
        sff_segment = stl_segment.convert()
        note=SimpleNote()
        note.name=f"{basename}_mito_{mito_number}_OMM"
        note.description=f"Screened poisson generated mesh depicting OMM (outer mitochondrial membrane) of mitochondria #{mito_number} of tomogram {basename}"
        note.external_references = [ommref]
        sff_segment = note.add_to_segment(sff_segment)        
        sff_segment.colour= sff.SFFRGBA(
            red=168/255,
            green=83/255,
            blue=160/255,
            alpha=1
        )
        seg.segment_list.append(sff_segment)
    
    if er_bool:
        stl_seg = stlreader.get_data(er)
        stl_segment = STLSegment(*stl_seg[0])
        sff_segment = stl_segment.convert()
        note=SimpleNote()
        note.name=f"{basename}_ER"
        note.description=f"Screened poisson generated mesh depicting ER (endoplasmic reticulum) of tomogram {basename}"
        note.external_references = [erref]
        sff_segment = note.add_to_segment(sff_segment) 
        seg.segment_list.append(sff_segment)
        sff_segment.colour= sff.SFFRGBA(
            red=155/255,
            green=216/255,
            blue=213/255,
            alpha=1
        )






    # Since the surfaces are nm instead of angstrom scale they need a transformation. 
    # 1 nm = 10Å so that leads to the 10!
    # Sometimes amira outputs stuff that is flipped in Z - you could put a -10 in the final array to fix that.
    # Our surfaces match with our voxel segs - we are prioritizing that for this submission.
    tx = np.array([[10,0,0,0], [0,10,0,0], [0,0,10,0]])
    transform = sff.SFFTransformationMatrix.from_array(tx)  
    seg.transform_list = sff.SFFTransformList()
    seg.transform_list.append(transform)


    ## Output is super simple
    seg.export("{}.hff".format(basename))