Bjagla at 19:34, 18 October 2007

2007-10-18T19:34:47Z

New page

{{BJNav}}
#!/usr/bin/perl
## Create PED files for PLINK
# C:\SAEC\scripts>perl makePED.pl -mapfile mapping_info.txt -snpfile snp_ids.txt -dir ../data -phenofile phenotypes.txt -outfile makePED.out

use Getopt::Long;
use Data::Dumper;
use strict;

my $VERSION = '1_00';
my $mappingFile = 'C:/SAEC/SJS Delivery from GSK/PGX40001_Illumina1M/Documents/Locus_Annotation_Files/illumina1M.map';
# Check if programm called correctly

my ($snpfp, $mapfp, $datadir, $phefp, $outfp);
my $result = GetOptions ("snpfile=s" => \$snpfp,
"mapfile=s" => \$mapfp,
"dir=s" => \$datadir,
"phenofile=s" => \$phefp,
"outfile=s" => \$outfp);
if($snpfp eq "" || $mapfp eq "" || $datadir eq "" || $phefp eq "" || $outfp eq "" ){
print <<EOF;
usage: perl makePED.pl
-snpfile <SNP ids (1st column)>
-mapfile <individual numbers (2nd column)>
-dir <directory with snp data files>
-phenofile <phenotype file (1st col Individual ID, 2nd col sex, 3rd col phenotype)>
-outfile <output file>\n
EOF
exit(-1);
}

open(OUTPUT2, '>', $outfp) || die("Could not open file $outfp!");
print2("#perl makePED.pl -snpfile $snpfp -mapfile $mapfp -dir $datadir -phenofile $phefp -outfile $outfp");
print2("#Version: $VERSION\n");

#read in information about mapping between Illumina individuals and GSK individuals
# file has been created from "PGX40001_GSK_SJS_B137_29Aug2007_DNAReport.xls"
# we are only interested in the second column
print2("#Reading Mapfile $mapfp\n");
open(DAT, $mapfp) || die("Could not open file $mapfp!");
my @map_data=<DAT>;
my %indi; #individual IDs
foreach my $m (@map_data) {
if ($m =~ !/^#/) {#ignore comment lines
$m =~ s/[\r\n]//g; #remove end of line characters
my @a = split(/\t/, $m);
$indi{lc $a[1]} = 1;
}
}
close(DAT);
if (!defined %indi) {
print2("#%indi not defined");
die("%indi not defined");
}
#print Dumper(\%indi);

#read in information about which SNPs should be there
print2("#Reading SNP file $mappingFile\n");
open(DAT, $mappingFile) || die("Could not open file $mappingFile!");
my @snp_data=<DAT>;
my %snps; #SNP names
foreach my $m (@snp_data) {
if ($m =~ /^#/) { #ignore comment lines
next;
}
$m =~ s/[\r\n]//g; #remove end of line characters
my @a = split(/\t/, $m);
$snps{$a[1]} = " 0 0"; #missing value
}
close(DAT);

if (!defined %snps) {
print2("#%snps not defined");
die("%snps not defined");
}

#read in information phenotypes
print2("#Reading Phenotype file $phefp\n");
open(DAT, $phefp) || die("Could not open file $phefp!");
my @phe_data=<DAT>;
my %phes; #phenotypic data plus SNP data

#print Dumper(\@phe_data);
#initialize, we are only interested in things that have been declared in $mapfp
foreach my $key (keys %indi) {
$phes{$key, "S"} = ""; # sex
$phes{$key, "P"} = ""; # phenotype
}
foreach my $m (@phe_data) {
#print $m ;
if ($m =~ /^#/) { #ignore comment lines
next;
}
$m =~ s/[\r\n]//g; #remove end of line characters
my @a = split(/\t/, $m);
#print "D1:" . Dumper(\@a);
if (exists $phes{$a[0], "S"}) {
$phes{$a[0], "S"} = $a[1];
$phes{$a[0], "P"} = $a[2];
}else {
print2("#Unknown phenotype " . $a[0] . "\n");
}

}
close(DAT);
#print Dumper(\%phes);
#exit(-1);

foreach my $key (keys %indi) {
if ($phes{$key, "S"} eq "") {
print2("#No Information for Individual $key\n");
next;
}

if (!open(DAT, $datadir . "/" . $key . ".txt")) {
print2("#No data file for $key\n");
next;
}
#initialize SNP hash
foreach my $m (keys %snps) {
$snps{$m} = " 0 0"; #missing value
}
print2( "#working on $key ");
# read the data
my $co=0;
while (<DAT>) {
#print "$_";
if (/^#/){next;} #ignore comment lines
s/[\r\n]//g; #remove end of line characters
my @line = split(/\t/);
if (!exists $snps{$line[0]}) {
print2("#Unknown SNP in Individual SNP file " . $line[0] . "\n");
next;
}
$co ++;
$line[1] =~ s/-/0/;
$line[2] =~ s/-/0/;
if ($line[1] !~ /[DIACGT0]/ || $line[2] !~ /[DIACGT0]/) {
print ("#$line[0] $line[1] $line[2]\n");
}
$snps{$line[0]} = " " . $line[1] . " " . $line[2];
#print Dumper(\@line) . "\n";
}
close(DAT);
print2("found $co SNPs\n");

#print out data
my $sex = '0';
my $phenotype = "-9";
if (uc $phes{$key, "S"} eq 'M') {
$sex = '1';
}
if (uc $phes{$key, "S"} eq 'F') {
$sex = '2';
}
if (uc $phes{$key, "P"} eq 'CO') {
$phenotype = '1';
}
if (uc $phes{$key, "P"} eq 'CA') {
$phenotype = '2';
}
my $co2=0;
foreach my $m (keys %snps) {
if ($snps{$m} eq " 0 0") {
$co2++;
}
}
print2("#found no info for $co2 SNPs\n");
print2("$key 1 0 0 $sex $phenotype ");
foreach my $m (sort keys %snps) {
#print this only to the outputfile
# print( OUTPUT2 "#$m");
print( OUTPUT2 $snps{$m});
}
print2 ("\n");

}

if (!defined %snps) {
print2("#%snps not defined");
die("%snps not defined");
}

#===============================

sub print2 {

my ($str) = @_;
print $str ;
print(OUTPUT2 $str);
}

MakePED - Revision history

Bjagla at 19:34, 18 October 2007