Mercurial > repos > dereeper > sniplay

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

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:47:56 -0500
parents 6bf69b40365c
children
comparison
equal deleted inserted replaced
14:d15869b3731a 15:31c23d943c29
23 23
24 die $usage 24 die $usage
25 if ( !$in || !$out); 25 if ( !$in || !$out);
26 26
27 27
28 my $plink_command = $PLINK_EXE . " --file $in --noweb --cluster --matrix --mds-plot 2 --out $out >>$in.plink.log 2>&1"; 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); 29 system($plink_command);
30 30
31 my $awk_cmd = "awk \{\'print \$1\'\} $in.ped"; 31 my @individuals = ();
32 my $inds = `$awk_cmd`;
33 my @individuals = split("\n",$inds);
34 32
35 my %populations; 33 my %populations;
36 if (-e "$in.individual_info.txt") 34 if (-e "$in.individual_info.txt")
37 { 35 {
38 open(my $I,"$in.individual_info.txt"); 36 open(my $I,"$in.individual_info.txt");
45 $populations{$ind} = $pop; 43 $populations{$ind} = $pop;
46 } 44 }
47 close($I); 45 close($I);
48 } 46 }
49 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
50 open(my $OUT,">$out.mds_plot.txt"); 56 open(my $OUT,">$out.mds_plot.txt");
51 my $go = 0; 57 my $go = 0;
58 print $OUT "Pop sample val1 val2 val3\n";
52 open(my $O,"$out.mds"); 59 open(my $O,"$out.mds");
60 my $numline = 0;
53 while(<$O>) 61 while(<$O>)
54 { 62 {
55 if ($go) 63 if ($go)
56 { 64 {
57 my $line = $_; 65 my $line = $_;
58 $line =~s/\n//g; 66 $line =~s/\n//g;
59 $line =~s/\r//g; 67 $line =~s/\r//g;
60 my @i = split(/\s+/,$line); 68 my @i = split(/\s+/,$line);
61 if ($line =~/^ /) 69 if ($line =~/^ /)
62 { 70 {
63 my $ind = $i[1]; 71 #my $ind = $i[1];
72 my $ind = $individuals[$numline];
64 my $pop = "Pop1"; 73 my $pop = "Pop1";
65 #if ($ind=~/^d/){$pop="Pop2";} 74 #if ($ind=~/^d/){$pop="Pop2";}
66 if ($populations{$ind}) 75 if ($populations{$ind})
67 { 76 {
68 $pop = $populations{$ind}; 77 $pop = $populations{$ind};
69 } 78 }
70 print $OUT "$pop $ind ".$i[4]." ".$i[5]."\n"; 79 print $OUT "$pop $ind ".$i[4]." ".$i[5]." ".$i[6]."\n";
71 } 80 }
72 if ($line =~/^\w/) 81 if ($line =~/^\w/)
73 { 82 {
74 my $ind = $i[0]; 83 #my $ind = $i[0];
84 my $ind = $individuals[$numline];
75 my $pop = "Pop1"; 85 my $pop = "Pop1";
76 #if ($ind=~/^d/){$pop="Pop2";} 86 #if ($ind=~/^d/){$pop="Pop2";}
77 if ($populations{$ind}) 87 if ($populations{$ind})
78 { 88 {
79 $pop = $populations{$ind}; 89 $pop = $populations{$ind};
80 } 90 }
81 print $OUT "$pop $ind ".$i[3]." ".$i[4]."\n"; 91 print $OUT "$pop $ind ".$i[3]." ".$i[4]." ".$i[5]."\n";
82 } 92 }
83 93 $numline++;
84 } 94 }
85 if (/C1/){$go = 1;} 95 if (/C1/){$go = 1;}
86 } 96 }
87 close($O); 97 close($O);
88 close($OUT); 98 close($OUT);