Dear rob234king,
I wrote a commented version of the script. To use the script you should have perl installed. place the files in the same folder as the script and on the command line write:
perl script.pl -file1 1.txt -file2 2.txt -cMlimit 4 -stepsize 2 -chr 1A
here, -file1 and -file2 refer to the options for your files 1.txt and 2.txt. -cMlimit is the maximum cM (be careful with big M) cM you want to compute to. -stepsize is the bin width of cM intervals, so if it is 2, than your bins will be 0-2, then 2-4 an so on. -chr is the chromosome name, in your case 1A. I commented the script and paste it here so everyone can use it. Please test it and tell me if there are mistakes etc. It might be useful to me later as well. Your files 1 and 2 need not to be ordered. Your result file will be in the same folder with name resultFile.txt.
I hope this helps,
Good luck with your research,
############################################CODE HERE###################################################
#!/usr/bin/perl
use strict;
use warnings;
use Getopt::Long;
#the global variables
my $stepCounter = 0;
my $stepSize;
my $chr;
my @chrList;
my $fileName1;
my $fileName2;
my $sampleNumber;
my $cMlimit;
#get the options
GetOptions("stepsize:f" => \$stepSize,"chr:s" => \$chr,"file1=s"=>\$fileName1,"file2=s"=>\$fileName2, "cMlimit=i"=>\$cMlimit);
$fileName1 = "./".$fileName1;
$fileName2 = "./".$fileName2;
#if chr is not defined it uses a standard list of human chromosome names
if (defined $chr) {
push(@chrList,$chr);
} else {
@chrList = (1..22,"X","Y");
}
$stepSize = defined $stepSize ? $stepSize : 1;
#open your file1 and file 2 which is in the same location of the script.
open(my $file1, '<',$fileName1) or die("cant open file 1!");
open(my $file2, '<',$fileName2) or die("cant open file 2!");
#skip the headers
skipHeader($file1);
skipHeader($file2);
#form a hash that will contain information for each scaffold the first key will be scaffold name, the entries will be name,chr,cM,sample1,sample2...
my %data = ();
while (<$file1>) {
chomp;
my @line = split(/\s+|\t+/,$_);
$data{$line[0]}{"name"} = $line[0];
for (my $i=1;$i<= $#line;$i++) {
$data{$line[0]}{"sample".$i} = $line[$i];
}
}
while (<$file2>) {
chomp;
my @line = split(/\s+|\t+/,$_);
$data{$line[0]}{"chr"} = $line[1];
$data{$line[0]}{"cM"} = $line[2];
}
#initiate a file handle which will contain your result.
open(my $resultFile, '>',"./resultFile.txt");
#print the first header..
print $resultFile "chr"."\t"."cM"."\t"."samples..\n";
#iterate over the hash
for (my $i = 0;$i<scalar(@chrList);$i++) {
$stepCounter = 0;
for (my $j = $stepCounter*$stepSize;$j<=$cMlimit; $j+=$stepSize) {
print $resultFile $chrList[$i]."\t".$stepCounter*$stepSize."\t";
for (my $k=1;$k <=$sampleNumber;$k++) {
my $sum = 0;
my $count = 0;
foreach my $key (keys %data) {
if($data{$key}{"chr"} eq $chrList[$i] && $data{$key}{"cM"}>=$stepSize*$stepCounter && $data{$key}{"cM"}<$stepSize*($stepCounter+1)) {
$sum += $data{$key}{"sample".$k};
$count++;
}
}
my $result;
if ($count != 0) {
$result = sprintf("%.1f",$sum/$count);
} else {
$result = 0;
}
if ($k == $sampleNumber) {
print $resultFile $result."\n";
} else {
print $resultFile $result."\t";
}
}
$stepCounter++;
}
}
sub skipHeader {
my $handle = $_[0];
if ($_[0] == $file1) {
my $line = <$handle>;
$sampleNumber = scalar(split(/\s+|\t+/,$line))-1;
} else {
<$handle>;
}
}
############################################CODE HERE###################################################
•
written
by
rob234king
