Mercurial > repos > dereeper > sniplay
view MDSplot/MDSbasedOnIBSmatrix.pl @ 15:31c23d943c29 draft default tip
Uploaded
author | dereeper |
---|---|
date | Tue, 08 Jan 2019 08:47:56 -0500 |
parents | 6bf69b40365c |
children |
line wrap: on
line source
#!/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);