Mercurial > repos > dereeper > sniplay

comparison MDSplot/MDSbasedOnIBSmatrix.pl.org @ 3:345f88a8f483 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author dereeper
date Fri, 10 Jul 2015 10:38:43 -0400
parents 420b57c3c185
children
comparison
equal deleted inserted replaced
2:feb40a9a8eae 3:345f88a8f483
1 #!/usr/bin/perl
2
3 use strict;
4 use Switch;
5 use Getopt::Long;
6 use Bio::SeqIO;
7
8 my $PLINK_EXE= "/apps/www/sniplay.cirad.fr/tools/plink/plink-1.07-x86_64/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
37 open(my $OUT,">$out.mds_plot.txt");
38 my $go = 0;
39 open(my $O,"$out.mds");
40 while(<$O>)
41 {
42 if ($go)
43 {
44 my $line = $_;
45 $line =~s/\n//g;
46 $line =~s/\r//g;
47 my @i = split(/\s+/,$line);
48 my $ind = $i[1];
49 print $OUT "$ind ".$i[4]." ".$i[5]."\n";
50 }
51 if (/C1/){$go = 1;}
52 }
53 close($O);
54 close($OUT);
55
56
57 my $j = 0;
58 open(my $IBS,">$out.ibs_matrix.txt");
59 print $IBS "Individuals " . join("\t",@individuals)."\n";
60 open(my $O2,"$out.mibs");
61 while(<$O2>)
62 {
63 my $line = $_;
64 $line =~s/\n//g;
65 $line =~s/\r//g;
66 my @i = split(/\s+/,$line);
67 print $IBS $individuals[$j]. " ". join("\t",@i)."\n";
68 $j++;
69 }
70 close($O2);
71 close($IBS);
72
73
74
75
76
77