df_wide <- data.frame(
  ID = rep(1:2, each=10),
  JN = rep(1:2, each=5),
  Frame = rep(1:5, 4),
  System = rep(1:2, 10),
  RKX = rep(1:10+rnorm(10,mean=1,sd=0.5),2),
  RKY = rep(1:10+rnorm(10,mean=1,sd=0.5),2),
  RKZ = rep(1:10+rnorm(10,mean=1,sd=0.5), 2),
  LHeX = rep(1:10-rnorm(10,mean=1,sd=0.5),2),
  LHeY = rep(1:10-rnorm(10,mean=1,sd=0.5),2),
  LHeZ = rep(1:10-rnorm(10,mean=1,sd=0.5),2))

head(df_wide, 2)
  ID JN Frame System      RKX      RKY      RKZ        LHeX       LHeY      LHeZ
1  1  1     1      1 1.332827 2.068720 2.295742 -0.02336031 -0.3011227 -1.212326
2  1  1     2      2 3.570076 3.306799 3.136177  2.08828231  1.9226740  2.106496

ID JN Frame System joint         X         Y         Z
1   1  1     1      1    RK  1.440103  2.221676  1.621871
2   1  1     1      1   LHe  3.537940  3.060948  2.856955

df_3D <- df_wide %>%
 gather(keys, values, -ID, -JN, -Frame, -System)%>%
  extract(keys, c("X", "Y", "Z", "joint"), "(X$) (Y$) (Z$) ([A-Z].$)")%>%
  spread(X, values)

df_wide %>%
    gather(keys, values, -ID, -JN, -Frame, -System) %>%
    extract(keys, c("Joint", "Coord"), "(.*)(X|Y|Z)$") %>%
    spread(Coord, values)

ID JN Frame System Joint          X          Y           Z
1   1  1     1      1   LHe  0.1344259 -0.2927277  0.05375166
2   1  1     1      1    RK  1.8083539  2.4053498  2.32899399
3   1  1     2      2   LHe  1.1777492  1.1780538  0.96549849
4   1  1     2      2    RK  3.2254236  2.4100235  2.79816371

library(tidyr)
library(stringr)

df2  <- df_wide %>%
  # leave the other columns
  gather( jointid, position, -ID, -JN, -Frame, -System ) %>%
  # insert a seperator to make it easier to split the X/Y/Z from the joint name
  mutate(jointid = str_replace( jointid, "X|Y|Z", ";\0")) %>%
  # split the joint name and the dimension apart
  tidyr::separate(jointid, c('joint', 'dim'), sep = ";" ) %>%
  # spread the joint and position apart into 3 columns
  spread(dim, position)