diff 2.4/script/Annotate_SoftSearch.pl @ 18:1163c16cb3c0 draft

Uploaded
author plus91-technologies-pvt-ltd
date Mon, 02 Jun 2014 07:35:53 -0400
parents e3609c8714fb
children
line wrap: on
line diff
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/2.4/script/Annotate_SoftSearch.pl	Mon Jun 02 07:35:53 2014 -0400
@@ -0,0 +1,40 @@
+#!/usr/bin/perl
+open(VCF,"$ARGV[0]")||die "Usage: <VCF> <Annotation.bed>\n\n\t\t The annotation BED should be of exons\n";
+
+$bedtools=`which intersectBed`;
+if(!$bedtools){die "Requires Bedtools in path\n\n"}
+if(!$ARGV[1]){die "Usage: <VCF> <Annotation.bed>\n\n";}
+
+while (<VCF>){
+	if($_=~/^#/){print;next}
+	chomp;
+	@data=split(/\t/,$_);
+	#Get left pair information
+	$chr1=$data[0];
+        $pos1a=$data[1]-1;
+        $pos1b=$data[1];
+	#Get right pair information
+	$data[4]=~s/[ACTGactghr\[\]]//g;#$data[4]=~s/hr/chr/;
+	@pos2=split(/:/,$data[4]);
+	$chr2="chr";
+	$chr2.=$pos2[0];
+	$pos2a=$pos2[1]-1;
+	$pos2b=$pos2[1];
+	#Now get left side annotations
+	#
+	#print "LEFT=get_anno($chr1,$pos1a,$pos1b)\n";
+	$left_gene=get_anno($chr1,$pos1a,$pos1b);
+        #print "RIGHT=get_anno($chr2,$pos2a,$pos2b)\n";
+        $right_gene=get_anno($chr2,$pos2a,$pos2b);
+	print "$_\t$left_gene\t$right_gene\n";
+}
+
+close VCF;
+
+sub get_anno(){
+	my ($chr,$pos1,$pos2)=@_;
+ 	$result=`perl -e 'if(($chr)&&($pos1)&&($pos2)){print join(\"\\t\",$chr,$pos1,$pos2,\"\\n\")}else {print STDERR "Not all variables defined: chr,pos1,pos2=$chr,$pos1,$pos2\n$_\n"}'|intersectBed -a $ARGV[1] -b stdin|cut -f4|head -1`;
+	$result=~chomp;$result=~s/\n//;
+	if(!$result){$result="NA"};
+	return $result;
+}