LqMGcNYx

1#!/usr/bin/perl -w
2use strict;
3
4# Given properly formatted herbarium records, this script will produce two
5# datasets of equal size for two time periods.  From these, for each species,
6# the mean value of each time period is calculated.  The script outputs a
7# tab-delimited table of columns for species, mean value of the first time 
8# period, and the mean value of the second time period.  This table can be used
9# to perform classical statistical tests in R.  Justification for the technique 
10# is provided below.
11#
12# Chris Grassa (7/13/2010)
13
14#   "When working with herbarium specimens, which were usually not sampled in a 
15# sytematic way, special consideration has to be given to potential biases
16# arising out of different sampling efforts.  To account for different sampling
17# efforts in the two time periods, we applied a random subsampling procedure
18# similiar to Warren et al. (2001) so that equal numbers of records were 
19# obtained in each 100 m band..."
20#   -Bergamini et al. (2009)
21
22#   "To equalize recorder effort, we subsampled the 1995-99 data by randomly
23# selecting the 1970-82 number of records from the 1995-99 data, subsampling
24# separately for each 100-km Ordnance Survey grid square to retain the broad
25# geographical distribution of 1970-82 records."
26#   -Warren et al. (2001)
27
28# open up our gbif/climate dataset.  the in and out filenames are passed as arguments when the script is executed.
29my $in = shift;
30my $out = shift;
31open IN, "<", $in;
32open OUT, ">", $out;
33
34# bin boundry dates are inclusive.
35# these dates make up the "before climate change" date range.
36my $bin1begin = 1900;
37my $bin1end = 1970;
38# these dates make up the "after climate change" date range.
39my $bin2begin = 1985;
40my $bin2end = 2010;
41
42# these hashes will contain our data from the two time periods.
43my %bin1;
44my %bin2;
45
46# for each record in our dataset.
47while ($in = <IN>)
48{
49    my $fileline = "$in";
50    chomp $fileline;
51    # extract pertinent data from pur records.  we could choose to analyze some variable other than elevation (latitude?) if we wanted to.
52    my @line = split /\t/,$fileline;
53    my $species = $line[0];
54    my $year = $line[1];
55    my $ele = $line[4];
56    
57    # reformat elevation to a five digit integer, with leading zeros if need be.
58    $ele = sprintf("%05d", "$ele");
59    # here we remove the last two digits of the elevation.  this number classifies the record into a 100 meter "band."
60    my $band = "$ele";
61    chop $band;
62    chop $band;
63    
64    # the dataset will be separated into two time period bins (before and after climate change).
65    # each bin will contain its records assigned to elevation bands.  we will use these bands when subsampling later.
66    # two hashes of arrays are used to map this data structure.
67    
68    if (($year >= $bin1begin) and ($year <= $bin1end))
69    # does the record belong in the "before(%bin1)" time period?
70    {
71        # if so, take the record($fileline) and push it onto the @array keyed by its $band.
72        push(@{ $bin1{"$band"}}, "$fileline");
73    }
74    elsif (($year >= $bin2begin) and ($year <= $bin2end))
75    # does the record belong in the "after(%bin2)" time period?
76    {
77        # if so, take the record($fileline) and push it onto the @array keyed by its $band.
78        push(@{ $bin2{"$band"}}, "$fileline");
79    }    
80}
81# now we have our dataset seperated into two time periods, each of which is divided by elevation.
82# this will allow us to subsample from the larger bin while simultaneously avoiding bias from sampling at higher elevations in later years. 
83
84#these arrays will hold our before and after binned datasets of equal size.
85my @before;
86my @after;
87
88#this array will be used to temporarily hold the dataset subsampled from the larger bin for each band.
89my @subsample;
90
91# there aren't many records in the upper bands, but i go through all bands that are earthly possible (mount everest is a little less than 8900 meters above sea level), to be conservative.
92my $n = 0;
93while ($n<90)
94{
95    # tempororary arrays hold the data extracted from the hashes for each band while we're doing comparisons and subsampling. 
96    my @band1;
97    my @band2;
98    # here i generate the three digit key corresponding to the elevation bands.
99    my $nnn = sprintf("%03d", "$n");
100    
101    # are there records in the current band for each bin?
102    if (exists $bin1{"$nnn"})
103    {
104        # if so, take them out, and put them in the temporary arrays.
105        @band1 = @{ $bin1{"$nnn"}};
106    }
107    if (exists $bin2{"$nnn"})
108    {
109        @band2 = @{ $bin2{"$nnn"}};
110    }
111    
112    # formally evaluate the size of each array.
113    my $count1 = scalar @band1;
114    my $count2 = scalar @band2;
115    # then compare, as we subsample from the larger of the two.
116    if ($count1<$count2)
117    {
118        # we can keep the smaller of the two elevation bands, so we just append it to our final array.
119        @before = (@before, @band1);
120        # this subroutine mathematically permutes the larger array in place.
121        FISHER_YATES_SHUFFLE( \@band2 );
122        # now subsample by taking a slice from the shuffled larger array that is equal in size to the smaller array. 
123        @subsample = @band2[1..$count1];
124        # append the subsample to the final array.
125        @after = (@after,@subsample);
126    }
127    
128    elsif ($count1>$count2)
129    {
130        @after = (@after, @band2);
131        FISHER_YATES_SHUFFLE( \@band1 );
132        @subsample = @band1[1..$count2];
133        @before = (@before,@subsample);
134    }
135    # no subsampling is needed if the elevation bands are of equal size in both bins.
136    else
137    {
138        @before = (@before, @band1);
139        @after = (@after, @band2);
140    }
141    # increment to the next band.
142    ++$n;
143}
144# now we have two datasets of equal size with approximately equal distributions of elevation.
145# what is their size?  
146my $count = scalar @before;
147################################################################################
148# uncomment this block to create a file containing all the subsampled data
149################################################################################
150#my $out2 = "after_subsampling";
151#open OUT2, ">", $out2;
152#foreach(@before)
153#{
154#    print OUT2 "$_\n";
155#}
156#foreach(@after)
157#{
158#    print OUT2 "$_\n";
159#}
160#`head -n $count after_subsampling | cut -f 5 >tmp1`;
161#`tail -n $count after_subsampling| cut -f 5 >tmp2`;
162#`paste tmp1 tmp2 >after_subsampling`;
163#`rm tmp1 tmp2`;
164#
165################################################################################
166
167
168# we will be comparing the mean elevevation for each species.  i again chose to use hashes of arrays.
169my %species_before;
170my %species_after;
171# i keep a running list of each species encountered
172my @specieslist;
173# for all the record lines
174for(1..$count)
175{
176    # take a record from each dataset.
177    my $beforeline= pop @before;
178    my $afterline= pop @after;
179    # split the fields.
180    my @before_data = split /\t/,$beforeline;
181    my @after_data = split /\t/,$afterline;
182    # extract the species name,
183    my $before_species_key = $before_data[0];
184    my $after_species_key = $after_data[0];
185    # and the elevation datum.
186    my $before_species_value = $before_data[4];
187    my $after_species_value = $after_data[4];
188    # push the elevation datum($before_species_value) onto an @array keyed by the species it belongs to($before_species_key) in its proper time period hash(%species_before).
189    push(@{ $species_before{"$before_species_key"}}, "$before_species_value");
190    push(@{ $species_after{"$after_species_key"}}, "$after_species_value");
191    # push both species names onto the running list of those encountered.
192    push @specieslist, ("$before_species_key","$after_species_key");        
193}
194
195
196# get an array of unique species names.  i didn't write the next four lines, but i did test them.
197my @unique_specieslist;
198# declare a hash.
199my %saw;
200# make sure it's undefined.
201undef %saw;
202# our unique species list will consist of the first encounter of each species as we iterate through the list of all species encountered in the previous block.
203@unique_specieslist = grep(!$saw{$_}++, @specieslist);
204
205# for each species in our dataset, we'll calculate its mean elevation before climate change, and after.
206foreach (@unique_specieslist)
207{
208    # i like to name the current string.
209    my $species = "$_";
210    # these arrays will temporarily hold the elevations extracted from the hashes for the current species.
211    my @beforevalues;
212    my @aftervalues;
213    # these strings will hold the mean elevation.
214    my $beforemean;
215    my $aftermean;
216    # if the species is represented for the time period,
217    if (exists $species_before{"$species"})
218    {
219        # extract its elevation data from the hash.
220        @beforevalues = @{ $species_before{"$species"}};
221    }
222    if (exists $species_after{"$species"})
223    {
224        @aftervalues = @{ $species_after{"$species"}};
225    }
226    # how many records do we have in each time period?
227    my $count_before = scalar @beforevalues;
228    my $count_after = scalar @aftervalues;
229    # we will analyze species with at least 12 records in each time period.
230    if (($count_before>0) and ($count_after>0))
231    {
232        # get the mean elevation before and after climate change.  we pass the array containing them to a subroutine.
233        $beforemean = MEAN(\@beforevalues);
234        $aftermean = MEAN(\@aftervalues);
235        # print the species name and both means. 
236        print OUT "$species\t$beforemean\t$aftermean\n";
237    }
238
239}
240
241# FISHER_YATES_SHUFFLE( \@array ) : generate a random permutation of @array in place.
242    # this code is from the Perl Cookbook ch 4.17; the algorithm's literal description was read from Wikipedia.
243# the fisher yates shuffle is unbiased, so that every permutation is equally likely.
244sub FISHER_YATES_SHUFFLE
245{
246    # take the passed array.
247    my $array = shift;
248    my $i;
249    # for i, from the array size down to one,
250    for ($i = @$array; --$i; )
251    {
252        # obtain a random integer, j that is 0<=j<=i.
253        my $j = int rand ($i+1);
254        # equivalent swaps do nothing.  wikipedia implies that this line may be waste.
255        next if $i == $j;
256        # exchange array element i for j.
257        @$array[$i,$j] = @$array[$j,$i];
258    }
259}
260
261# calculate the mean of an array of numbers.
262sub MEAN
263{
264    # make sure the subroutine is passed a defined array.
265    @_ == 1 or die ('Sub usage: $mean = mean(\@array);'); 
266    # take the passed array.
267    my ($array_ref) = @_; 
268    my $sum;
269    # get its size.
270    my $count = scalar @$array_ref; 
271    # sum the values contained in the array
272    foreach (@$array_ref)
273    {
274        $sum += $_;
275    } 
276    # divide the array's sum by its size to return its mean.
277    return $sum / $count; 
278}