Mercurial > repos > plus91-technologies-pvt-ltd > softsearch

comparison 2.4/script/Annotate_SoftSearch.pl @ 13:e3609c8714fb draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author plus91-technologies-pvt-ltd
date Fri, 30 May 2014 03:37:55 -0400
parents
children
comparison
equal deleted inserted replaced
12:980fce54e60a 13:e3609c8714fb
1 #!/usr/bin/perl
2 open(VCF,"$ARGV[0]")||die "Usage: <VCF> <Annotation.bed>\n\n\t\t The annotation BED should be of exons\n";
3
4 $bedtools=`which intersectBed`;
5 if(!$bedtools){die "Requires Bedtools in path\n\n"}
6 if(!$ARGV[1]){die "Usage: <VCF> <Annotation.bed>\n\n";}
7
8 while (<VCF>){
9 if($_=~/^#/){print;next}
10 chomp;
11 @data=split(/\t/,$_);
12 #Get left pair information
13 $chr1=$data[0];
14 $pos1a=$data[1]-1;
15 $pos1b=$data[1];
16 #Get right pair information
17 $data[4]=~s/[ACTGactghr\[\]]//g;#$data[4]=~s/hr/chr/;
18 @pos2=split(/:/,$data[4]);
19 $chr2="chr";
20 $chr2.=$pos2[0];
21 $pos2a=$pos2[1]-1;
22 $pos2b=$pos2[1];
23 #Now get left side annotations
24 #
25 #print "LEFT=get_anno($chr1,$pos1a,$pos1b)\n";
26 $left_gene=get_anno($chr1,$pos1a,$pos1b);
27 #print "RIGHT=get_anno($chr2,$pos2a,$pos2b)\n";
28 $right_gene=get_anno($chr2,$pos2a,$pos2b);
29 print "$_\t$left_gene\t$right_gene\n";
30 }
31
32 close VCF;
33
34 sub get_anno(){
35 my ($chr,$pos1,$pos2)=@_;
36 $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`;
37 $result=~chomp;$result=~s/\n//;
38 if(!$result){$result="NA"};
39 return $result;
40 }