Your description of what you need to do is pretty much the same, but the missing step is that you need to get the glyph closure from only the icon glyphs and not the a-z glyphs (seen here: https://github.com/rsheeter/subset-gf-icons/blob/main/src/subset_gf_icons/subset_gf_icons.py#L106). This can be done in harfbuzz by generating a subset plan and then extracting the set of included glyphs from it (via: hb_subset_plan_create_or_fail and hb_subset_plan_old_to_new_glyph_mapping)
The result of that closure would then be used to form the final subset in combination with the a-z gids and no layout closure.