comparison ngsap-vc/varscan/varscan_somatic.pl @ 3:0d10255b5434 draft default tip

Uploaded

3:0d10255b5434

#!/usr/bin/perl


use strict;
use Cwd;

die qq(
Bad numbr of inputs

) if(!@ARGV);

my $options ="";
my $normal="";
my $command="";
my $tumor="";
my $output="";
my $working_dir = cwd();
my $snp = "$working_dir/output.snp.vcf";
my $indels = "$working_dir/output.indel.vcf";

foreach my $input (@ARGV) 
{
	my @tmp = split "::", $input;
	if($tmp[0] eq "COMMAND") 
	{
		$command = $tmp[1];
	} 
	if($tmp[0] eq "NORMAL") 
	{
		$normal = $tmp[1];
	} 
	elsif($tmp[0] eq "TUMOR") 
	{
		$tumor = $tmp[1];
	}
	elsif($tmp[0] eq "OPTION") 
	{
		$options = "$options ${tmp[1]}";
	}
	elsif($tmp[0] eq "OUTPUT") 
	{
		$output = $tmp[1];
	}
	
	else 
	{
		die("Unknown Input: $input\n");
	}
}

system ("$command $normal $tumor $options ");
system("grep -v '^\#' $indels | grep -v '^chrom position' >> $snp");

my @chr_ord = chromosome_order($tumor);

vs2vcf($snp, $output,\@chr_ord);


sub vs2vcf 
{

	#
	# G l o b a l     v a r i a b l e s 
	#
	my $version = '0.1';

	#
	# Read in file
	#
	my $input = shift;
	my $output = shift;
	my $chr_ord = shift;
	open(IN, $input) or die "Can't open $input': $!\n";
	open(OUT, ">$output") or die "Can't create $output': $!\n";
	my %output;

	while ( <IN> )
	{
		if ( /^#/ )
		{
			print OUT;
			next;
		}
		chomp;
		my $line = $_;

		my @flds = split ( "\t", $line );
		my $ref = $flds[3];
		my $alt = $flds[4];
		#
		# Deletion of bases
		#
		if ( $alt =~ /^\-/ )
		{
			($flds[3], $flds[4]) = ($ref.substr($alt,1), $ref);
		}

		#
		# Insertion of bases
		#
		if ( $alt =~ /^\+/ )
		{
			$flds[4] = $ref.substr($alt,1);
		}
		print OUT join( "\t", @flds),"\n" unless defined $chr_ord;
		$output{$flds[0]}{$flds[1]} = join( "\t", @flds)."\n" if defined $chr_ord;
	}
	close(IN);
	# if chromosome order given return in sorted order
	if(defined $chr_ord) 
	{
		for my $chrom (@{ $chr_ord }) 
		{
			for my $pos (sort {$a<=>$b} keys %{ $output{$chrom} }) 
			{
				print OUT $output{$chrom}{$pos};
			}
		}
	}
	close(OUT);
}


sub chromosome_order 
{
	my $input = shift;
	# calculate flagstats
	my $COMM = "samtools view -H $input | grep '^\@SQ'";
	my @SQ = `$COMM`;
	chomp @SQ;
	for(my $i = 0; $i <= $#SQ; $i++) 
	{
		$SQ[$i] =~ s/^\@SQ\tSN:(.*?)\tLN:\d+$/$1/;
	} 
	return(@SQ);
}