|
13
|
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 }
|
