annotate MDSplot/MDSbasedOnIBSmatrix.pl @ 0:3e19d0dfcf3e draft

Uploaded
author dereeper
date Mon, 23 Mar 2015 05:57:27 -0400
parents
children 6bf69b40365c
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
1 #!/usr/bin/perl
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
2
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
3 use strict;
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
4 use Switch;
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
5 use Getopt::Long;
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
6 use Bio::SeqIO;
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
7
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
8 my $PLINK_EXE= "plink";
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
9
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
10 my $usage = qq~Usage:$0 <args> [<opts>]
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
11 where <args> are:
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
12 -i, --in <input>
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
13 -o, --out <output>
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
14 ~;
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
15 $usage .= "\n";
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
16
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
17 my ($in,$out);
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
18
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
19
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
20 GetOptions(
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
21 "in=s" => \$in,
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
22 "out=s" => \$out
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
23 );
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
24
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
25 die $usage
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
26 if ( !$in || !$out);
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
27
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
28
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
29 my $plink_command = $PLINK_EXE . " --file $in --noweb --cluster --matrix --mds-plot 2 --out $out >>$in.plink.log 2>&1";
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
30 system($plink_command);
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
31
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
32 my $awk_cmd = "awk \{\'print \$1\'\} $in.ped";
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
33 my $inds = `$awk_cmd`;
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
34 my @individuals = split("\n",$inds);
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
35
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
36 my %populations;
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
37 if (-e "$in.individual_info.txt")
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
38 {
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
39 open(my $I,"$in.individual_info.txt");
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
40 while(<$I>)
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
41 {
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
42 my $line = $_;
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
43 $line =~s/\n//g;
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
44 $line =~s/\r//g;
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
45 my ($ind,$pop) = split(/;/,$line);
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
46 $populations{$ind} = $pop;
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
47 }
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
48 close($I);
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
49 }
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
50
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
51 open(my $OUT,">$out.mds_plot.txt");
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
52 my $go = 0;
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
53 open(my $O,"$out.mds");
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
54 while(<$O>)
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
55 {
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
56 if ($go)
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
57 {
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
58 my $line = $_;
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
59 $line =~s/\n//g;
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
60 $line =~s/\r//g;
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
61 my @i = split(/\s+/,$line);
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
62 if ($line =~/^ /)
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
63 {
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
64 my $ind = $i[1];
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
65 my $pop = "Pop1";
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
66 if ($populations{$ind})
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
67 {
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
68 $pop = $populations{$ind};
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
69 }
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
70 print $OUT "$pop $ind ".$i[4]." ".$i[5]."\n";
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
71 }
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
72 if ($line =~/^\w/)
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
73 {
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
74 my $ind = $i[0];
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
75 my $pop = "Pop1";
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
76 if ($populations{$ind})
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
77 {
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
78 $pop = $populations{$ind};
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
79 }
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
80 print $OUT "$pop $ind ".$i[3]." ".$i[4]."\n";
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
81 }
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
82
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
83 }
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
84 if (/C1/){$go = 1;}
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
85 }
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
86 close($O);
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
87 close($OUT);
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
88
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
89
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
90 my $j = 0;
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
91 open(my $IBS,">$out.ibs_matrix.txt");
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
92 print $IBS "Individuals " . join("\t",@individuals)."\n";
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
93 open(my $O2,"$out.mibs");
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
94 while(<$O2>)
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
95 {
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
96 my $line = $_;
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
97 $line =~s/\n//g;
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
98 $line =~s/\r//g;
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
99 my @i = split(/\s+/,$line);
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
100 print $IBS $individuals[$j]. " ". join("\t",@i)."\n";
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
101 $j++;
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
102 }
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
103 close($O2);
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
104 close($IBS);
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
105
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
106
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
107
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
108
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
109
3e19d0dfcf3e Uploaded
dereeper
parents:
diff changeset
110