annotate 2.4/script/Bam2pair.pl @ 13:e3609c8714fb draft

Uploaded
author plus91-technologies-pvt-ltd
date Fri, 30 May 2014 03:37:55 -0400
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
13
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1 #!/usr/bin/perl
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
2 #Author Steven Hart, PhD
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
3 #11-15-2012
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
4 #Convert and filter BAM files into merged bed
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
5 #Output should be
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
6 #chrA StartA EndA chrB StartB EndB Gene_id #supportingReads StrandA StrandB
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
7 #chr9 1000 5000 chr9 3000 3800 bedpe_example2 100 - +
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
8
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
9 use Cwd;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
10 use File::Basename;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
11 #Usage
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
12 sub usage(){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
13 print "Usage: perl Bam2Pair.pl -b <BAM> -o <outfile>\n
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
14 -isize [10000]\t\tThe insert size to be considered discordant\n
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
15 -winsize [10000]\tThe distance between mate pairs to be considered the same\n
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
16 -min [1]\t\tThe minimum number of reads required to support an SV event\n
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
17 -prefix need a random prefix so files with the same name don't get created\n\n"
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
18 ;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
19 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
20 $bedtools=`which intersectBed`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
21 $samtools=`which samtools`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
22
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
23 if(!defined($bedtools)){die "BEDtools must be installed\n";}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
24 if(!defined($samtools)){die "Samtools must be installed\n";}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
25 use Getopt::Long;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
26 #Declare variables
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
27 GetOptions(
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
28 'b=s' => \$BAM_FILE, #path to bam
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
29 'out=s' => \$outfile, #path to output
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
30 'java:s' => \$java ,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
31 'chrom:s' => \$chrom ,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
32 'isize=i' => \$isize,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
33 'winsize=i' => \$winsize,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
34 'prefix=s' => \$prefix,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
35 'min=i' => \$minSupport,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
36 'blacklist:s' => \$new_blacklist,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
37 'q=s' => \$qual,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
38 'v' => \$verbose
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
39 );
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
40 #if(defined($picard_path)){$picard_path=$picard_path} else {die "Must specify a path to PICARD so that files can be sorted and indexed properly\n"};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
41 if(!defined($BAM_FILE)){die "Must specify a BAM file!\n".usage();}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
42 if(!defined($outfile)){die "Must specify an out filename!\n".usage();}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
43 if(!defined($java)){$java=$java;}else{$java=`which java`}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
44 if(!defined($qual)){$qual=20}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
45 if($new_blacklist){$new_blacklist=" -L $new_blacklist"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
46
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
47
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
48 $Filter_BAM=$BAM_FILE;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
49
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
50 @bam=split("/",$Filter_BAM);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
51 $Filter_BAM=@bam[@bam-1];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
52 $Filter_BAM=~s/.bam/$prefix.bam/;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
53 $Filter_sam=$Filter_BAM;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
54 $Filter_sam=~s/.bam/.sam/;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
55
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
56
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
57
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
58
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
59 print "\nLooking for Discordant read pairs (and Unmated reads) without soft-clips\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
60
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
61 #$command=join("","samtools view -h -q 20 -f 1 -F 1804 ",$BAM_FILE," ",$chrom," ",$new_blacklist," | awk -F\'\\t\' \'{if (\$9 > ", $isize, " || \$9 < -",$isize," || \$9 == 0 || \$1 ~ /^@/) print \$0}' > ",$Filter_sam);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
62
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
63
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
64 #Change command to allow reads where mate is unmapped & remove qual
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
65 $command=join("","samtools view -h -f 1 -F 1800 -q $qual ",$BAM_FILE," ",$chrom," ",$new_blacklist," | awk -F\'\\t\' \'{if (\$9 > ", $isize, " || \$9 < -",$isize," || \$9 == 0 || \$1 ~ /^@/) print \$0}' > ",$Filter_sam);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
66
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
67 print "$command\n" if ($verbose);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
68 system($command);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
69 $path = dirname(__FILE__);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
70 $Filter_cluster=$Filter_sam;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
71 $Filter_cluster=~s/.sam/.cluster/;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
72 $command=join("",$path,"/ReadCluster.pl -i=$Filter_sam -o=$Filter_cluster -m $minSupport");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
73 if($verbose){print "\n$command\n"};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
74
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
75 system($command);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
76
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
77 ##################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
78 #Now there are 2 SAM files of filtered reads
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
79 #.filter.cluster.inter.sam
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
80 #.filter.cluster.intra.sam
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
81 $result_pe=join("",$Filter_cluster,".out");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
82 $command=join("","cat ",$Filter_cluster,".int\*|perl -ane 'next if(\@F[0]=~/^\@/);if(\@F[6]!~/=/){print join(\"\\t\",\$F[11],\@F[2],\@F[3],\@F[6],\@F[7],\"\\n\")}else{print join(\"\\t\",\$F[11],\@F[2],\@F[3],\@F[2],\@F[7],\"\\n\")}' >",$result_pe);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
83 if($verbose){print $command."\n"};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
84 system($command);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
85 #my ($sample, $chrstart, $start, $chrend, $end)
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
86 $command=join("","cat ",$result_pe," | ",$path,"/cluster.pair.pl ",$winsize," |awk '(\$6 >",$minSupport,")' >> ", $outfile);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
87 if($verbose){print $command."\n"};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
88 system($command);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
89 $filt1=join("",$Filter_cluster,".inter.sam");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
90 $filt2=join("",$Filter_cluster,".intra.sam");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
91
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
92
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
93 unlink($Filter_sam,$filt1,$filt2,$result_pe);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
94
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
95 #########################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
96 #Now determine if left or righ clipping surrogate
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
97 print "\nBam2pair.pl Done\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
98