Mercurial > repos > dereeper > sniplay

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

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