annotate MDSplot/MDSbasedOnIBSmatrix.pl @ 15:31c23d943c29 draft default tip

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