Tuesday, January 1, 2013

Processing .tps morphometric files in R

This has absolutely nothing to do with phylogenetics, but today (New Year's day 2013!) I've been starting the year with some R programming to pull data out of the raw output files from a program tpsDIG2 (by Jim Rohlf) that we use to get data from landmarks on digital skeletal x-rays of lizards.

We are just getting linear measurements at the moment (no fancy geometric morphometrics), so our data analysis is pretty simple. Here's what a skeletal x-ray of a lizard with the landmark coordinates overlain looks like (click for larger version):
(There is also a scale bar, not shown.)

The output from tpsDIG looks as follows (some lines omitted):
LM=43
726.00000 1580.00000
773.00000 1457.00000
679.00000 1456.00000
840.00000 1493.00000
840.00000 1463.00000
...
IMAGE=KMW_122_adj.tif
SCALE=0.125544

The first line gives the number of landmarks; subsequent lines give the horizontal & vertical pixel positions for each landmark; and the final lines give the image information & scale. Obviously, in R we want to read in the file, translate it to meaningful (i.e., non-pixel) scale, and perhaps compute our desired linear measurements. We might also want to be able to visually inspect our data for errors.

I have written two functions tps.process (to pull data out of the file and translate the scale); and tps.postProcess, which computes the linear measures we want from our input files. The first function is thus "general" (suitable for many people getting linear measures from .tps files); whereas the latter is specific to our data, although it could probably be modified and re-purposed easily. The code is posted here.

Here's a quick example. Note that I have also included a mode tps.process(...,type="anole") that visualizes the landmark data (for error checking) in a sensible way specifically for our x-rays.

> list.files()
[1] "KMW_004.tps" "KMW_122.tps" "KMW_317.tps"
[4] "tps.process.R"
> files<-list.files(pattern=".tps")
> files
[1] "KMW_004.tps" "KMW_122.tps" "KMW_317.tps"
> source("tps.process.R")
> X<-tps.process(files,type="anole")
OK? Press <Enter> to continue...
OK? Press <Enter> to continue...
OK? Press <Enter> to continue...
One of the plots above is created for each input .tps file so that we can easily review our input data for errors. (You can see it even looks kind of like a lizard!)

> # Here's our data rescaled
> X
$KMW_004adj
         V1        V2
1  28.338576 75.445244
2  33.807424 59.287284
3  21.751100 59.784452
4  49.343924 62.643168
5  47.355252 59.287284
6   6.214600 60.157328
...
$KMW_122_adj
         V1        V2
1  29.126208 80.724792
2  35.026776 65.282880
...
> # now let's automatically post-process
> LL<-tps.postProcess(X)
> LL
                rJL      lJL       JW    rMETC    lMETC
KMW_004adj  17.05837 16.98986 12.06657 3.900868 3.890955
KMW_122_adj 16.53086 16.64819 11.80180 3.766320 3.694494
KMW_317adj  15.80653 16.18932 11.41486 3.283738 3.391124
               rRAD     lRAD      rULN      lULN     rHUM
KMW_004adj  8.728803 8.672875  9.743256  9.175747 11.57585
KMW_122_adj 9.094794 9.109513 10.081895 10.313728 12.43139
KMW_317adj  8.239660 8.533229  9.663571  9.765426 12.13156
               lHUM     PECT     PELV     rFEM     lFEM
KMW_004adj  12.01011 8.082803 6.840578 16.44463 16.42960
KMW_122_adj 12.19394 7.910268 6.872889 17.09659 16.83975
KMW_317adj  11.85941 7.864227 6.990424 15.13729 15.12956
               rTIB     lTIB     rFIB     lFIB   rMETT1
KMW_004adj  13.59564 13.38724 13.41548 13.33405 5.133726
KMW_122_adj 14.58960 14.54415 14.47245 14.17035 5.280316
KMW_317adj  12.61329 12.73256 12.48540 12.61020 4.745144
             rMETT2   lMETT1   lMETT2
KMW_004adj  7.990535 4.991837 8.342392
KMW_122_adj 8.799730 5.176312 8.956841
KMW_317adj  7.743424 4.383258 7.899812

Now we have left & right measures for all the characters in our x-ray, above. Cool.

4 comments:

  1. This comment has been removed by the author.

    ReplyDelete
  2. Really cool Liam!
    Just to share its existence, package geomorph has the function readland.tps(file) to read tps files as well and package Momocs allows the placement of landmarks in case someone would like to try and place the landmarks in R.
    Ricardo

    ReplyDelete
  3. Hi Liam
    Thank you for the great blog, I am trying to import tps file, like you show in the above output from tpsDIG. Could you please tell how to import tsp file into R? infromation under the data is look like Image=D0015.JPG ID=Speciemen1 Scale=0.00458 This format for all individuls (79). I need to transfer all data to R so can do further analysis
    Thank you
    Abdullah

    ReplyDelete
  4. Hi again
    I tried to apply R code you posted, but its seems difficult as it it long and I don know which is related to filename and specific comments on it. Could you please put some details. Thanks
    Abdullah

    ReplyDelete