|
13
|
1 #!/usr/bin/perl
|
|
|
2
|
|
|
3 use strict;
|
|
|
4 use Getopt::Long;
|
|
|
5 use Bio::SeqIO;
|
|
|
6
|
|
|
7 my $PLINK_EXE= "plink";
|
|
|
8
|
|
|
9 my $usage = qq~Usage:$0 <args> [<opts>]
|
|
|
10 where <args> are:
|
|
|
11 -i, --in <input>
|
|
|
12 -o, --out <output>
|
|
|
13 ~;
|
|
|
14 $usage .= "\n";
|
|
|
15
|
|
|
16 my ($in,$out);
|
|
|
17
|
|
|
18
|
|
|
19 GetOptions(
|
|
|
20 "in=s" => \$in,
|
|
|
21 "out=s" => \$out
|
|
|
22 );
|
|
|
23
|
|
|
24 die $usage
|
|
|
25 if ( !$in || !$out);
|
|
|
26
|
|
|
27
|
|
|
28 my $plink_command = $PLINK_EXE . " --vcf $in --cluster --allow-extra-chr --matrix --mds-plot 3 --out $out >>$in.plink.log 2>&1";
|
|
|
29 system($plink_command);
|
|
|
30
|
|
|
31 my @individuals = ();
|
|
|
32
|
|
|
33 my %populations;
|
|
|
34 if (-e "$in.individual_info.txt")
|
|
|
35 {
|
|
|
36 open(my $I,"$in.individual_info.txt");
|
|
|
37 while(<$I>)
|
|
|
38 {
|
|
|
39 my $line = $_;
|
|
|
40 $line =~s/\n//g;
|
|
|
41 $line =~s/\r//g;
|
|
|
42 my ($ind,$pop) = split(/;/,$line);
|
|
|
43 $populations{$ind} = $pop;
|
|
|
44 }
|
|
|
45 close($I);
|
|
|
46 }
|
|
|
47
|
|
|
48
|
|
|
49 my $line_ind = `grep CHROM $in`;
|
|
|
50 $line_ind =~s/\n//g;$line_ind =~s/\r//g;
|
|
|
51 my @tab = split(/\t/,$line_ind);
|
|
|
52 for (my $i = 9; $i <= $#tab; $i++){
|
|
|
53 push(@individuals,$tab[$i]);
|
|
|
54 }
|
|
|
55
|
|
|
56 open(my $OUT,">$out.mds_plot.txt");
|
|
|
57 my $go = 0;
|
|
|
58 print $OUT "Pop sample val1 val2 val3\n";
|
|
|
59 open(my $O,"$out.mds");
|
|
|
60 my $numline = 0;
|
|
|
61 while(<$O>)
|
|
|
62 {
|
|
|
63 if ($go)
|
|
|
64 {
|
|
|
65 my $line = $_;
|
|
|
66 $line =~s/\n//g;
|
|
|
67 $line =~s/\r//g;
|
|
|
68 my @i = split(/\s+/,$line);
|
|
|
69 if ($line =~/^ /)
|
|
|
70 {
|
|
|
71 #my $ind = $i[1];
|
|
|
72 my $ind = $individuals[$numline];
|
|
|
73 my $pop = "Pop1";
|
|
|
74 #if ($ind=~/^d/){$pop="Pop2";}
|
|
|
75 if ($populations{$ind})
|
|
|
76 {
|
|
|
77 $pop = $populations{$ind};
|
|
|
78 }
|
|
|
79 print $OUT "$pop $ind ".$i[4]." ".$i[5]." ".$i[6]."\n";
|
|
|
80 }
|
|
|
81 if ($line =~/^\w/)
|
|
|
82 {
|
|
|
83 #my $ind = $i[0];
|
|
|
84 my $ind = $individuals[$numline];
|
|
|
85 my $pop = "Pop1";
|
|
|
86 #if ($ind=~/^d/){$pop="Pop2";}
|
|
|
87 if ($populations{$ind})
|
|
|
88 {
|
|
|
89 $pop = $populations{$ind};
|
|
|
90 }
|
|
|
91 print $OUT "$pop $ind ".$i[3]." ".$i[4]." ".$i[5]."\n";
|
|
|
92 }
|
|
|
93 $numline++;
|
|
|
94 }
|
|
|
95 if (/C1/){$go = 1;}
|
|
|
96 }
|
|
|
97 close($O);
|
|
|
98 close($OUT);
|
|
|
99
|
|
|
100
|
|
|
101 my $j = 0;
|
|
|
102 open(my $IBS,">$out.ibs_matrix.txt");
|
|
|
103 print $IBS "Individuals " . join("\t",@individuals)."\n";
|
|
|
104 open(my $O2,"$out.mibs");
|
|
|
105 while(<$O2>)
|
|
|
106 {
|
|
|
107 my $line = $_;
|
|
|
108 $line =~s/\n//g;
|
|
|
109 $line =~s/\r//g;
|
|
|
110 my @i = split(/\s+/,$line);
|
|
|
111 print $IBS $individuals[$j]. " ". join("\t",@i)."\n";
|
|
|
112 $j++;
|
|
|
113 }
|
|
|
114 close($O2);
|
|
|
115 close($IBS);
|
|
|
116
|
|
|
117
|
|
|
118
|
|
|
119
|
|
|
120
|
|
|
121
|
