| Previous changeset 12:88748d846a20 (2018-08-14) Next changeset 14:d15869b3731a (2019-01-08) |
|
Commit message:
Uploaded |
|
added:
MDSbasedOnIBSmatrix.pl |
| b |
| diff -r 88748d846a20 -r 734a3572c1d6 MDSbasedOnIBSmatrix.pl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/MDSbasedOnIBSmatrix.pl Tue Jan 08 08:45:34 2019 -0500 |
| [ |
| @@ -0,0 +1,121 @@ +#!/usr/bin/perl + +use strict; +use Getopt::Long; +use Bio::SeqIO; + +my $PLINK_EXE= "plink"; + +my $usage = qq~Usage:$0 <args> [<opts>] +where <args> are: + -i, --in <input> + -o, --out <output> +~; +$usage .= "\n"; + +my ($in,$out); + + +GetOptions( + "in=s" => \$in, + "out=s" => \$out +); + +die $usage + if ( !$in || !$out); + + +my $plink_command = $PLINK_EXE . " --vcf $in --cluster --allow-extra-chr --matrix --mds-plot 3 --out $out >>$in.plink.log 2>&1"; +system($plink_command); + +my @individuals = (); + +my %populations; +if (-e "$in.individual_info.txt") +{ + open(my $I,"$in.individual_info.txt"); + while(<$I>) + { + my $line = $_; + $line =~s/\n//g; + $line =~s/\r//g; + my ($ind,$pop) = split(/;/,$line); + $populations{$ind} = $pop; + } + close($I); +} + + +my $line_ind = `grep CHROM $in`; +$line_ind =~s/\n//g;$line_ind =~s/\r//g; +my @tab = split(/\t/,$line_ind); +for (my $i = 9; $i <= $#tab; $i++){ + push(@individuals,$tab[$i]); +} + +open(my $OUT,">$out.mds_plot.txt"); +my $go = 0; +print $OUT "Pop sample val1 val2 val3\n"; +open(my $O,"$out.mds"); +my $numline = 0; +while(<$O>) +{ + if ($go) + { + my $line = $_; + $line =~s/\n//g; + $line =~s/\r//g; + my @i = split(/\s+/,$line); + if ($line =~/^ /) + { + #my $ind = $i[1]; + my $ind = $individuals[$numline]; + my $pop = "Pop1"; + #if ($ind=~/^d/){$pop="Pop2";} + if ($populations{$ind}) + { + $pop = $populations{$ind}; + } + print $OUT "$pop $ind ".$i[4]." ".$i[5]." ".$i[6]."\n"; + } + if ($line =~/^\w/) + { + #my $ind = $i[0]; + my $ind = $individuals[$numline]; + my $pop = "Pop1"; + #if ($ind=~/^d/){$pop="Pop2";} + if ($populations{$ind}) + { + $pop = $populations{$ind}; + } + print $OUT "$pop $ind ".$i[3]." ".$i[4]." ".$i[5]."\n"; + } + $numline++; + } + if (/C1/){$go = 1;} +} +close($O); +close($OUT); + + +my $j = 0; +open(my $IBS,">$out.ibs_matrix.txt"); +print $IBS "Individuals " . join("\t",@individuals)."\n"; +open(my $O2,"$out.mibs"); +while(<$O2>) +{ + my $line = $_; + $line =~s/\n//g; + $line =~s/\r//g; + my @i = split(/\s+/,$line); + print $IBS $individuals[$j]. " ". join("\t",@i)."\n"; + $j++; +} +close($O2); +close($IBS); + + + + + + |