# HG changeset patch
# User dereeper
# Date 1546955134 18000
# Node ID 734a3572c1d6fe80f9057d28113584dc9b318a8a
# Parent  88748d846a208a4caafb04165888889713b45505
Uploaded

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);
+
+
+	
+	
+
+