Mercurial > repos > dereeper > sniplay

comparison MDSbasedOnIBSmatrix.pl @ 13:734a3572c1d6 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author dereeper
date Tue, 08 Jan 2019 08:45:34 -0500
parents
children
comparison
equal deleted inserted replaced
12:88748d846a20 13:734a3572c1d6
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