annotate 2.4/script/SoftSearch.multi.pl @ 16:8eb7d93f7e58 draft

Uploaded
author plus91-technologies-pvt-ltd
date Sat, 31 May 2014 11:23:36 -0400
parents e3609c8714fb
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
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
3 ####
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
4 #### Usage: SoftSearch.pl [-lqrmsd] -b <BAM> -f <Genome.fa> -sam <samtools path> -bed <bedtools path>
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
5 #### Created 1-30-2012 by Steven Hart, PhD
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
6 #### hart.steven@mayo.edu
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
7 #### Required bedtools & samtools to be in path
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
8
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
9 #use FindBin;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
10 #use lib "$FindBin::Bin/lib";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
11 use lib "/data2/bsi/reference/softsearch/lib" ;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
12 use Getopt::Long;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
13 use strict;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
14 use warnings;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
15 use Data::Dumper;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
16 use LevD;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
17 use File::Basename;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
18 use List::Util qw(min max);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
19
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
20 my (@INPUT_BAM,$INPUT_FASTA,$OUTPUT_FILE,$minSoft,$minSoftReads,$dist_To_Soft,$bedtools,$samtools);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
21 my ($minRP, $temp_output, $num_sd, $MapQ, $chrom, $unmated_pairs, $minBQ, $pair_only, $disable_RP_only);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
22 my ($levD_local_threshold, $levD_distl_threshold,$pe_upper_limit,$high_qual,$sv_only,$blacklist,$genome_file,$verbose);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
23
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
24 my $cmd = "";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
25
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{1,}' => \@INPUT_BAM,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
29 'f=s' => \$INPUT_FASTA,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
30 'o:s' => \$OUTPUT_FILE,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
31 'm:i' => \$minRP,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
32 'l:i' => \$minSoft,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
33 'r:i' => \$minSoftReads,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
34 't:i' => \$temp_output,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
35 's:s' => \$num_sd,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
36 'd:i' => \$dist_To_Soft,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
37 'q:i' => \$MapQ,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
38 'c:s' => \$chrom,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
39 'u:s' => \$unmated_pairs,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
40 'x:s' => \$minBQ,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
41 'p' => \$pair_only,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
42 'g' => \$disable_RP_only, #ignore softclips
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
43 'j:s' => \$levD_local_threshold,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
44 'k:s' => \$levD_distl_threshold,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
45 'a:s' => \$pe_upper_limit,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
46 'e:s' => \$high_qual,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
47 'L' => \$sv_only,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
48 'v' => \$verbose,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
49 'blacklist:s' => \$blacklist,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
50 'genome_file:s' => \$genome_file,
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
51 "help|h|?" => \&usage);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
52 #print "Using @INPUT_BAM as INPUT_BAM\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
53 unless($sv_only){$sv_only=""};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
54 my $INPUT_BAM=$INPUT_BAM[0];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
55 #print "MY NEW INPUT BAM=$INPUT_BAM[0]\n\n";die;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
56 if(defined($INPUT_BAM)){$INPUT_BAM=$INPUT_BAM} else {print usage();die "Where is the BAM file?\n\n"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
57 if(defined($INPUT_FASTA)){$INPUT_FASTA=$INPUT_FASTA} else {print usage();die "Where is the fasta file?\n\n"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
58 my ($fn,$pathname) = fileparse($INPUT_BAM,".bam");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
59 #my $index=`ls $pathname/$fn*bai|head -1`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
60 #my $index =`ls \${INPUT_BAM%.bam}*bai`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
61 #print "INDEX=$index\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
62 #if(!$index){die "\n\nERROR: you need index your BAM file\n\n"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
63 my $index="";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
64 ### get current time
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
65 print "Start Time : " . &spGetCurDateTime() . "\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
66 my $now = time;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
67
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
68 #if(defined($OUTPUT_FILE)){$OUTPUT_FILE=$OUTPUT_FILE} else {$OUTPUT_FILE="output.vcf"; print "\nNo outfile specified. Using output.vcf as default\n\n"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
69 if(defined($minSoft)){$minSoft=$minSoft} else {$minSoft=5}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
70 if(defined($minRP)){$minRP=$minRP} else {$minRP=5}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
71 if(defined($minSoftReads)){$minSoftReads=$minSoftReads} else {$minSoftReads=5}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
72 if(defined($dist_To_Soft)){$dist_To_Soft=$dist_To_Soft} else {$dist_To_Soft=300}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
73 if(defined($num_sd)){$num_sd=$num_sd} else {$num_sd=6}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
74 if(defined($MapQ)){$MapQ=$MapQ} else {$MapQ=20}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
75
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
76 unless (defined $pe_upper_limit) { $pe_upper_limit = 10000; }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
77 unless (defined $levD_local_threshold) { $levD_local_threshold = 0.05; }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
78 unless (defined $levD_distl_threshold) { $levD_distl_threshold = 0.05; }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
79 #Get sample name if available
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
80 my $SAMPLE_NAME="";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
81 my $OUTNAME ="";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
82 $SAMPLE_NAME=`samtools view -f2 -H $INPUT_BAM|awk '{if(\$1~/^\@RG/){sub("ID:","",\$2);name=\$2;print name}}'|head -1`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
83 $SAMPLE_NAME=~s/\n//g;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
84 if (!$OUTPUT_FILE){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
85 if($SAMPLE_NAME ne ""){$OUTNAME=$SAMPLE_NAME.".vcf"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
86 else {$OUTNAME="output.vcf"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
87 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
88 else{$OUTNAME=$OUTPUT_FILE}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
89
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
90 print "Writing results to $OUTNAME\n";
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 ##Make sure if submitting on SGE, to prepned the "chr". Not all referecne FAST files require "chr", so we shouldn't force the issue.
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
94 if(!defined($chrom)){$chrom=""}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
95 if(!defined($unmated_pairs)){$unmated_pairs=0}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
96
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
97 my $badQualValue=chr($MapQ);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
98 if(defined($minBQ)){ $badQualValue=chr($minBQ); }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
99
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
100 if($badQualValue eq "#"){$badQualValue="\#"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
101
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
102 # adding and cheking for samtools and bedtools in the PATh
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
103 ## check for bedtools and samtools in the path
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
104 $bedtools=`which intersectBed` ;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
105 if(!defined($bedtools)){die "\nError:\n\tno bedtools. Please install bedtools and add to the path\n";}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
106 #$samtools=`samtools 2>&1`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
107 $samtools=`which samtools`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
108 if($samtools !~ /(samtools)/i){die "\nError:\n\tno samtools. Please install samtools and add to the path\n";}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
109
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
110 print "Usage = SoftSearch.pl -l $minSoft -q $MapQ -r $minSoftReads -d $dist_To_Soft -m $minRP -s $num_sd -c $chrom -b @INPUT_BAM -f $INPUT_FASTA -o $OUTNAME \n\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
111 sub usage {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
112 print "\nusage: SoftSearch.pl [-cqlrmsd] -b <BAM> -f <Genome.fa> \n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
113 print "\t-q\t\tMinimum mapping quality [20]\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
114 print "\t-l\t\tMinimum length of soft-clipped segment [5]\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
115 print "\t-r\t\tMinimum depth of soft-clipped reads at position [5]\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
116 print "\t-m\t\tMinimum number of discordant read pairs [5]\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
117 print "\t-s\t\tNumber of sd away from mean to be considered discordant [6]\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
118 print "\t-u\t\tNumber of unmated pairs [0]\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
119 print "\t-d\t\tMax distance between soft-clipped segments and discordant read pairs [300]\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
120 print "\t-o\t\tOutput file name [output.vcf]\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
121 print "\t-t\t\tPrint temp files for debugging [no|yes]\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
122 print "\t-c\t\tuse only this chrom or chr:pos1-pos2\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
123 print "\t-p\t\tuse paired-end mode only \n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
124 print "\t-g\t\tEnable paired-only seach. This will look for discordant read pairs even without soft clips.\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
125 print "\t-a\t\tset the minimum distance for a discordant read pair without soft-clipping info [10000]\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
126 print "\t-L\t\tFlag to print out even small deletions (low quality)\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
127 print "\t-e\t\tdisable strict quality filtering of base qualities in soft-clipped reads [no]\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
128 print "\t-blacklist\tareas of the genome to skip calling. Requires -genome_file\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
129 print "\t-genome_file\ttab seperated value of chromosome name and length. Only used with -blacklist option\n\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
130
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
131 exit 1;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
132 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
133
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
134
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
135 #############################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
136 # create temporary variable name
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
137 #############################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
138 srand (time ^ $$ ^ unpack "%L*", `ps axww | gzip -f`);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
139 our $random_name=join "", map { ("a".."z")[rand 26] } 1..8;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
140
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
141 #############################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
142 ## create green list
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
143 ##############################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
144 #
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
145 my $new_blacklist="";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
146 if($blacklist){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
147 if(!$genome_file){die "if using a blacklist, you must also specify the location of a genome_file
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
148 The format of the genome_file should be
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
149 chrom size
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
150 chr1 249250621
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
151 chr2 243199373
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
152 ...
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
153
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
154 If using hg19, you can ge the genome file by
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
155 mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -e \"select chrom, size from hg19.chromInfo\" > hg19.genome";}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
156
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
157 $cmd=join("","complementBed -i $blacklist -g $genome_file >",$random_name,".bed") ;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
158 system ($cmd);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
159 $new_blacklist=join(""," -L ",$random_name,".bed ");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
160 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
161
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
162 if($verbose){print "CMD=$cmd\nBlacklist is $new_blacklist\n";}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
163
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
164
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
165
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
166
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
167
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
168 #############################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
169 # Calcualte insert size distribution of properly mated reads
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
170 #############################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
171
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
172 #Change for compatability with other operating systems
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
173 #my $metrics=`samtools view -q $MapQ -f2 $INPUT_BAM $chrom|cut -f9|head -10000|awk '{if (\$1<0){\$1=-\$1}else{\$1=\$1} sum+=\$1; sumsq+=\$1*\$1} END {print sum/NR, sqrt(sumsq/NR - (sum/NR)**2)}'`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
174 #print "samtools view -q $MapQ -f2 $INPUT_BAM $chrom|cut -f9|head -10000|awk '{if (\$1<0){\$1=-\$1}else{\$1=\$1} sum+=\$1; sumsq+=\$1*\$1} END {print sum/NR, sqrt(sumsq/NR - (sum/NR)^2)}'\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
175 my $metrics=`samtools view -q $MapQ -f2 $INPUT_BAM $chrom|cut -f9|head -10000|awk '{if (\$1<0){\$1=-\$1}else{\$1=\$1} sum+=\$1; sumsq+=\$1*\$1} END {print sum/NR, sqrt(sumsq/NR - (sum/NR)^2)}'`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
176 #my ($mean,$stdev)=split(/ /,$metrics);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
177 my ($mean,$stdev)=split(/\s/,$metrics);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
178 $stdev=~s/\n//;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
179
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
180 #print "MEAN=$mean\tSTDEV=$stdev\n\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
181
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
182 my $upper_limit=int($mean+($num_sd*$stdev));
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
183 my $lower_limit=int($mean-($num_sd*$stdev));
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
184 die if (!$mean);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
185 print qq{The mean insert size is $mean +/- $stdev (sd)
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
186 The upper limit = $upper_limit
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
187 The lower limit = $lower_limit\n
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
188 };
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
189 if($lower_limit<0){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
190 print "Warning!! Given this insert size distribution, we can not call small indels. No other data will be affected\n\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
191 $lower_limit=1;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
192 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
193 my $tmp_name=join ("",$random_name,".tmp.bam");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
194 my $random_file_sc = "";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
195 my $command = "";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
196
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
197 #############################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
198 # Make sam file that has soft clipped reads
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
199 #############################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
200 #give file a name
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
201 if(!defined($pair_only)){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
202 foreach my $bam(@INPUT_BAM){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
203 $random_file_sc=join ("",$random_name,".sc.sam");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
204 $command=join ("","samtools view -q $MapQ -F 1024 $bam $chrom $new_blacklist| awk '{OFS=\"\\t\"}{c=0;if(\$6~/S/){++c};if(c == 1){print}}' | perl -ane '\$TR=(\@F[10]=~tr/\#//);if(\$TR<2){print}' >> ", $random_file_sc);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
205 print "Making SAM file of soft-clipped reads from $bam\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
206 if($verbose){ print "$command\n";}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
207 system("$command");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
208 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
209 #############################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
210 # Find areas that have deep enough soft-clip coverage
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
211 print "Identifying soft-clipped regions that are at least $minSoft bp long iin $random_file_sc\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
212 open (FILE,"$random_file_sc")||die "Can't open soft-clipped sam file $random_file_sc\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
213
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
214 my $tmpfile=join("",$random_file_sc,".sc.passfilter");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
215 open (OUT,">$tmpfile")||die "Can't write files here!\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
216
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
217 while(<FILE>){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
218 @_ = split(/\t/, $_);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
219 #### parse CIGAR string and create a hash of array of each operation
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
220 my @CIGAR = split(/([0-9]+[SMIDNHXP])/, $_[5]);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
221 my $hash;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
222 map { push(@{$hash->{$2}}, $1) if (/(\d+)([SMIDNHXP])/) } @CIGAR;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
223
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
224 #for ($i=0; $i<=$#softclip_pos; $i++) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
225 foreach my $softclip (@{$hash->{S}}) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
226 #if ($CIGAR[$softclip_pos[$i]] > $minSoft){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
227 if ($softclip > $minSoft){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
228 ###############Make sure base qualities don't have more than 2 bad marks
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
229 my $qual=$_[10];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
230 my $TR=($qual=~tr/$badQualValue//);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
231 if($badQualValue eq "#"){ $TR=($qual=~tr/\#//); }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
232 #Skip the soft clip if there is more than 2 bad qual values
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
233 #next if($TR > 2);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
234 # if (!$high_qual){next if($TR > 2);}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
235 print OUT;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
236 last;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
237 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
238 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
239 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
240 close FILE;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
241 close OUT;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
242
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
243 $command=join(" ","mv",$tmpfile,$random_file_sc);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
244 if($verbose){ print "$command\n";}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
245 system("$command");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
246 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
247
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
248 #########################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
249 #Stack up SoftClips
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
250 #########################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
251 my $random_file=join("",$random_name,".sc.direction.bed");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
252 if(!defined($pair_only)){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
253 open (FILE,"$random_file_sc")|| die "Can't open sam file\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
254 #$random_file=join("",$random_name,".sc.direction");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
255
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
256 print "Calling sides of soft-clips from $random_file_sc\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
257 #\nTMPOUT=$random_file\tINPUT=$random_file_sc\n\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
258 open (TMPOUT,">$random_file")|| die "Can't create tmp file\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
259
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
260 while (<FILE>){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
261 @_ = split(/\t/, $_);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
262 #### parse CIGAR string and create a hash of array of each operation
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
263 my @CIGAR = split(/([0-9]+[SMIDNHXP])/, $_[5]);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
264 my $hash;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
265 map { push(@{$hash->{$2}}, $1) if (/(\d+)([SMIDNHXP])/) } @CIGAR;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
266
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
267 #### next if softclips on each end
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
268 next if ($_[5] =~ /^[0-9]+S.*S$/);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
269
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
270 #### next softclip occurs in the middle
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
271 next if ($_[5] =~ /^[0-9]+[^S][0-9].*S.+$/);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
272
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
273 my $softclip = $hash->{S}[0];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
274
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
275 my $end1 = 0;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
276 my $end2 = 0;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
277 my $softBases = "";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
278 my $right_corrected="";my $left_corrected="";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
279
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
280 if ($softclip > $minSoft) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
281
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
282 ####If the soft clip occurs at end of read and its on the minus strand, then it's a right clip
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
283 if ($_[5] =~ /^.*S$/) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
284 $end1=$_[3]+length($_[9])-$softclip-1;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
285 $end2=$end1+1;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
286 next if ($end1<0);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
287 #RIGHT clip on Minus
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
288 $softBases=substr($_[9], length($_[9])-$softclip, length($_[9]));
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
289 #Right clips don't always get clipped correctly, so fix that
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
290 # Check to see if sc base matches ref
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
291 $right_corrected=baseCheck($_[2],$end2,"right",$softBases);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
292 print TMPOUT "$right_corrected\n"
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
293
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
294 } else {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
295 #### Begins with S (left clip)
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
296 $end1=$_[3]-$softclip;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
297 next if ($end1<0);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
298
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
299 $softBases=substr($_[9], 0,$softclip);#print "TMP=$softBases\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
300 $left_corrected=baseCheck($_[2],$end1,"left",$softBases);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
301 if(!$left_corrected){print "baseCheck($_[2],$end1,left,$softBases)\n";next}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
302 print TMPOUT "$left_corrected\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
303 #print "\nSEQ=$_[9]\t\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
304
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
305 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
306 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
307 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
308 close FILE;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
309 close TMPOUT;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
310 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
311 sub baseCheck{
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
312 my ($chrom,$pos,$direction,$softBases)=@_;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
313 #skip if position is less than 0, which is caused by MT DNA
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
314 return if ($pos<0);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
315 my $exit="";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
316
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
317 while(!$exit){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
318 if($direction=~/right/){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
319 my $refBase=getSeq($chrom,$pos,$INPUT_FASTA);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
320 my $softBase=substr($softBases,0,1);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
321 if ($softBase !~ /$refBase/){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
322 my $value=join("\t",$chrom,$pos,$pos+1,join("|",$softBases,$direction));
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
323 $exit="STOP";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
324 return $value;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
325 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
326 else{
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
327 $pos=$pos+1;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
328 $softBases=substr($softBases, 1,length($softBases));
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
329 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
330 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
331 else{
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
332 my $refBase=getSeq($chrom,$pos+1,$INPUT_FASTA);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
333 my $softBase=substr($softBases,-1,1);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
334 if ($softBase !~ /$refBase/){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
335 $pos=$pos-1+length($softBases);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
336 my $value=join("\t",$chrom,$pos-1,$pos,join("|",$softBases,$direction));
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
337 $exit="STOP";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
338 return $value;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
339 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
340 else{
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
341 $pos=$pos-1;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
342 $softBases=substr($softBases, 0, -1);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
343 #print "Trying again $softBases\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
344 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
345
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
346 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
347
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
348 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
349 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
350 #Remove SAM files to conserve space
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
351 unlink($random_file_sc);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
352
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
353
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
354
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
355 ###
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
356 #
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
357 ######################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
358 # Transform Read pair groups into softclip equivalents
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
359 ######################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
360 #
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
361 #
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
362 #
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
363 my $v="";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
364 #if($disable_RP_only){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
365 print "Running Bam2pair.pl\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
366 print "Looking for discordant read pairs without requiring soft-clipping information\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
367 use FindBin qw($Bin);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
368 my $path=$Bin;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
369 # print"\n\nPATH=$path\n\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
370 if($verbose){$v="-v"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
371 foreach my $random_file_disc(@INPUT_BAM){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
372 my $tmp_out=join("",$random_name,"RP.out");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
373 $command=join("","perl ",$path,"/Bam2pair.pl -b $random_file_disc -o $tmp_out -isize $pe_upper_limit -winsize $dist_To_Soft -min $minRP -chrom $chrom -prefix $random_name -q $MapQ -blacklist $random_name.bed $v");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
374 if($verbose){ print "$command\n"};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
375 system("$command");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
376 $command=join("","perl -ane '\$end1=\@F[1];\$end2=\@F[3];print join(\"\\t\",\@F[0..1],\$end1,\"unknown|left\");print \"\\n\";print join(\"\\t\",\@F[2..3],\$end2,\"unknown|left\");print \"\\n\"' ", $tmp_out," >> ",$random_file);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
377 if($verbose){print "$command\n"};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
378 system($command);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
379 unlink($tmp_out);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
380 #}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
381 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
382
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
383
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
384 ######################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
385 unlink("$random_file","$tmp_name","$random_file","$index","$random_name","$new_blacklist") if (-z $random_file || ! -e $random_file) ;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
386 if (-z $random_file || ! -e $random_file){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
387 print "Softclipped file is empty($random_file).\nNo soft clipping found using desired paramters\n\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
388 open (OUT,">$OUTNAME")||die "Can't write files here!\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
389 &print_header();
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
390 close OUT;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
391 exit;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
392 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
393
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
394
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
395 #############################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
396 # Make sure there are enough soft-clippped supporting reads
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
397 #############################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
398 my $outfile=join("",$random_file,".sc.merge.bed");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
399 #sortbed -i .sc.direction | mergeBed -nms -d 25 -i stdin > .sc.merge.bed
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
400 $command=join(" ","sortBed -i",$random_file," | mergeBed -nms -i stdin","|grep \";\"","|awk '{OFS=\"\t\"}(NF==4)'",">",$outfile);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
401
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
402 #print "$command\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
403 system("$command");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
404
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
405 if (-z $outfile || ! -e $outfile){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
406 unlink("$tmp_name","$random_file","$outfile","$index","$random_name","$new_blacklist");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
407 print "mergeBed file is empty.\nNo strucutral variants found\n\n" ;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
408 open (OUT,">$OUTNAME")||die "Can't write files here!\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
409 &print_header();
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
410 close OUT;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
411 exit;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
412 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
413
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
414 print "completed mergeBed\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
415
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
416 ###############################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
417 # If left and right are on the same line, make into 2 lines
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
418 ###############################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
419 open (INFILE,$outfile)||die "couldn't open temp file : $. \n\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
420 my $tmp2=join("",$random_name,".sc.fixed.merge.bed");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
421 #print "INFILE=$outfile\tOUTFILE=$tmp2\n\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
422 #INPUT FORMAT=chr9\t131467\t131473\tATGCTTATTAAAA|left;TTATTAAAAGCATA|left
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
423 open (OUTFILE,">$tmp2")||die "couldn't create temp file : $. \n\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
424 while(<INFILE>){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
425 chomp $_;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
426 my $l = $_;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
427
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
428 my @a = split(/\t/, $l);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
429 my $info = $a[3];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
430 my @info_arr = split(/\;/, $info);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
431 my @left_arr=();
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
432 my @right_arr=();
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
433 @left_arr = grep(/left/, @info_arr);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
434 @right_arr = grep(/right/, @info_arr);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
435
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
436 #New
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
437 my $left = join(";", @left_arr);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
438 my $right = join(";", @right_arr);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
439 $info = join(";", @info_arr);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
440
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
441 if((@left_arr) && (@right_arr)){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
442 print OUTFILE "$a[0]\t$a[1]\t$a[2]\t$left\n$a[0]\t$a[1]\t$a[2]\t$right\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
443 } else{
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
444 my $all=join("\t",@a[0..2],$info);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
445 print OUTFILE "$all\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
446 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
447 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
448
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
449 # make sure output file name is $outfile
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
450 $command=join(" ","sed -e '/ /s//\t/g'", $tmp2,"|awk 'BEGIN{OFS=\"\\t\"}(NF==4)'", "|perl -pne 's/ /\t/g'>",$outfile);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
451 system("$command");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
452 if($verbose){print "$command\n"};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
453 unlink("$tmp_name","$random_file","$tmp2","$outfile","$index","random_name","$new_blacklist") if (-z $outfile || ! -e $outfile) ;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
454 if (-z $outfile || ! -e $outfile){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
455 print "Fixed mergeBed file is empty($outfile).\nNo strucutral variants found\n\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
456 open (OUT,">$OUTNAME")||die "Can't write files here!\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
457 &print_header();
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
458 close OUT;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
459 exit;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
460 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
461
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
462 print "completed fixing mergeBed\n\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
463
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
464 ###############################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
465 # Seperate directions of soft clips
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
466 ###############################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
467 my $left_sc = join("", "left", $tmp2);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
468 my $right_sc = join("", "right", $tmp2);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
469 use FindBin qw($Bin);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
470 #my $path=$Bin;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
471
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
472 $command=join("","grep left ", $tmp2, " |sed -e '/left /s//left\;/g' | sed -e '/ /s//\t/g'|perl ".$path."/direction_filter.pl - >",$left_sc);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
473 system("$command");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
474 #print "$command\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
475 $command=join("","grep right ", $tmp2, " |sed -e '/right /s//right\;/g' | sed -e '/ /s//\t/g'|perl ".$path."/direction_filter.pl - >",$right_sc);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
476 #$command=join(" ","grep right ", $tmp2, " |sed -e '/right /s//right\;/g' | sed -e '/ /s//\t/g' >",$right_sc);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
477 system("$command");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
478 #print "$command\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
479 #die "CHECK $right_sc\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
480
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
481 ###############################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
482 # Count the number and identify directions of soft clips
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
483 ###############################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
484 print "Count the number and identify directions of soft clips\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
485 #print "looking in $outfile\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
486 $outfile=join("",$random_name,".sc.fixed.merge.bed");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
487 #system("ls -lhrt");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
488 open (INFILE,$outfile)||die "couldn't open temp file\n\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
489 my $tmp3 = join("", $random_file, "predSV");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
490 open (OUTFILE, ">$tmp3")||die "couldn't create temp file\n\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
491 while(<INFILE>){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
492 chomp;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
493 @_=split(/\t/,$_);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
494 my $count=tr/\;//;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
495 $count=$count+1;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
496 my $left=0;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
497 my $right=0;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
498
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
499 while ($_ =~ /left/g) { $left++ } # count number of right clips
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
500 while ($_ =~ /right/g) { $right++ } # count number of left clips
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
501
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
502 ###############################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
503 if ($count >= $minSoftReads){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
504 ####get longets soft-clipped read
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
505 my @clips=split(/\;|\|/,$_[3]);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
506
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
507 my ($max, $temp, $temp2, $temp3, $dir, $maxSclip) = (0) x 6;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
508 for (my $i=0; $i<$count; $i++) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
509 my $plus1=$i+1;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
510 $temp=length($clips[$i]);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
511 $temp2=$clips[$plus1];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
512 $temp3=$clips[$i];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
513
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
514 if ($temp > $max){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
515 $maxSclip=$temp3;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
516 $max =$temp;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
517 $dir=$temp2;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
518 } else {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
519 $max=$max;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
520 $dir=$dir;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
521 $maxSclip=$maxSclip;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
522 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
523 $i++;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
524 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
525 my $order2 = join("|", $left, $right);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
526 #print join ("\t",@_[0..2],$maxSclip,$max,$dir,$count,$order2) . "\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
527 print OUTFILE join ("\t",@_[0..2],$maxSclip,$max,$dir,$count,$order2) . "\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
528 } elsif($_=~/unknown/){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
529 print OUTFILE join ("\t",@_[0..2],"NA","NA","left","NA","NA|NA") . "\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
530 print OUTFILE join ("\t",@_[0..2],"NA","NA","right","NA","NA|NA") . "\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
531 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
532 ####Format is Chrom,start, end,longest Soft-clip,length of longest Soft-clip, direction of longest soft-clip,#supporting softclips,#right Sclips|#left Sclips
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
533 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
534 close INFILE;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
535 close OUTFILE;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
536
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
537 unlink("$tmp2","$tmp_name","$random_file","$tmp3","$outfile","$index","$random_name","$right_sc","$left_sc","$new_blacklist") if (-z $tmp3 || !-e $tmp3) ;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
538
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
539 if (-z $tmp3 || !-e $tmp3){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
540 print "No structural variants found while Counting the number and identify directions of soft clips.\n" ;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
541
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
542 # open (OUT,">$OUTNAME")||die "Can't write files here!\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
543 # &print_header();
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
544 # close OUT;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
545 exit;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
546 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
547
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
548 print "Done counting Softclipped reads\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
549 ###############################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
550 #### Print header information
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
551 ###############################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
552
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
553
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
554 foreach my $random_file_disc(@INPUT_BAM){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
555 print "Making the header for $random_file_disc\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
556 $SAMPLE_NAME=`samtools view -f2 -H $random_file_disc|awk '{if(\$1~/^\@RG/){sub("ID:","",\$2);name=\$2;print name}}'|head -1`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
557 $SAMPLE_NAME=~s/\n//g;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
558 if($chrom){$SAMPLE_NAME.=".".$chrom}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
559
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
560 $SAMPLE_NAME.=".vcf";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
561 open (OUT,">$SAMPLE_NAME")||die "Can't write files here!\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
562 &print_header();
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
563
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
564 # DO the bulk of the work
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
565 open (FILE,"$tmp3")|| die "Can't open file\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
566
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
567 while (<FILE>){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
568 #If left clip {+- or -- or -+ }{+- are uninformative b/c they go upstream}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
569 #If right clip {++ or -- or +-}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
570 chomp $_;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
571 my @res=();my $res;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
572 my $line = $_;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
573 my @info = split(/\t/, $_);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
574 my $i=0;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
575 my $basename=basename($random_file_disc);$i=0;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
576 if($info[5] eq "left") {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
577 $res=bulk_work("left", $line, $random_file_disc);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
578 if(!$res){$res=join("\t",".",".",".",".",".",".",".",".",".",".")};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
579 $i++;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
580 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
581 elsif ($info[5] eq "right") {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
582 $res=bulk_work("right", $line, $random_file_disc);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
583 if(!$res){$res=join("\t",".",".",".",".",".",".",".",".",".",".")};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
584 $i++;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
585 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
586 if($res){@res=split("\t",$res);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
587 print OUT join("\t",@res)."\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
588 }}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
589 close FILE;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
590 close OUT;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
591 print "Done with $random_file_disc\n\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
592 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
593
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
594
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
595
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
596 ###############################################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
597 ###############################################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
598 #### Delete temp files
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
599 my $meregedBed=join("",$random_name,".sc.direction.bed.sc.merge.bed");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
600
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
601 if(defined($temp_output)){$temp_output=$temp_output} else {$temp_output="no"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
602
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
603 if ($temp_output eq "no"){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
604 unlink("$tmp_name","$random_file","$tmp2",,"$tmp3","$outfile","$index","$random_name","$right_sc","$left_sc","$meregedBed","$random_name.bed");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
605 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
606 ####Sort VCF
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
607 #my $tmp=join(".",$random_name,"tmp");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
608 #Get header
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
609 #$cmd="grep \"#\" $OUTNAME > $tmp";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
610 #system($cmd);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
611 #sort results
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
612 #$cmd="grep -v \"#\" $OUTNAME|perl -pne 's/chr//'|sort -k1,1n -k2,2n|perl -ne 'print \"chr\".\$_' >>$tmp";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
613 #system($cmd);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
614 #$cmd="mv $tmp $OUTNAME";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
615 #system($cmd);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
616 #remove entries next to each other
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
617
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
618
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
619 print "Analysis Completed\n\nYou did it!!!\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
620 print "Finish Time : " . &spGetCurDateTime() . "\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
621 $now = time - $now;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
622 printf("\n\nTotal running time: %02d:%02d:%02d\n\n", int($now / 3600), int(($now % 3600) / 60),
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
623 int($now % 60));
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
624
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
625 exit;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
626
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
627 ###############################################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
628 sub rev_comp {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
629 my $dna = shift;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
630 my $revcomp = reverse($dna);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
631 $revcomp =~ tr/ACGTacgt/TGCAtgca/;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
632 return $revcomp;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
633 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
634
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
635
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
636 ###############################################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
637 #### to get reference base
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
638 sub getSeq{
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
639 my ($chr,$pos,$fasta)=@_;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
640 #don't require chr
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
641 #if($chr !~ /^chr/){die "$chr is not correct\n";}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
642 # die "$pos is not a number\n" if ($pos <0);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
643 my @result=();
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
644 if ($pos <0){print "$pos is not a valid position (likely caused by circular MT chromosome)\n";return;}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
645
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
646 @result = `samtools faidx $fasta $chr:$pos-$pos`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
647 if($result[1]){chomp($result[1]);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
648 return uc($result[1]);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
649 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
650 return("NA");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
651 #### after return will not be printed
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
652 ####print "RESULTS=@result\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
653 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
654
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
655 sub getBases{
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
656 my ($chr,$pos1,$pos2,$fasta)=@_;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
657 #don't require chr
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
658 #if($chr !~ /^chr/){die "$chr is not correct\n";}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
659 my @result=();
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
660 if ($pos1 <0){print "$pos1 is not a valid position (likely caused by circular MT chromosome)\n";return;};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
661
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
662 @result = `samtools faidx $fasta $chr:$pos1-$pos2`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
663 if(!$result[1]){$result[1]="NA"};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
664 chomp($result[1]);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
665 return uc($result[1]);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
666
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
667 #### after return will not be printed
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
668 ####print "RESULTS=@result\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
669 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
670 ###############################################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
671 #### to get time
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
672 sub spGetCurDateTime {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
673 my ($sec, $min, $hour, $mday, $mon, $year) = localtime();
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
674 my $curDateTime = sprintf "%4d-%02d-%02d %02d:%02d:%02d",
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
675 $year+1900, $mon+1, $mday, $hour, $min, $sec;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
676 return ($curDateTime);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
677 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
678
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
679
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
680 ###############################################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
681 #### print header
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
682 sub print_header {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
683 my $date=&spGetCurDateTime();
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
684 my $header = qq{##fileformat=VCFv4.1
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
685 ##fileDate=$date
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
686 ##source=SoftSearch.pl
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
687 ##reference=$INPUT_FASTA
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
688 ##Usage= SoftSearch.pl -l $minSoft -q $MapQ -r $minSoftReads -d $dist_To_Soft -m $minRP -u $unmated_pairs -s $num_sd -b @INPUT_BAM -f $INPUT_FASTA -o $OUTNAME
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
689 ##INFO=<ID=SVTYPE,Number=1,Type=String,Description="Type of structural variant">
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
690 ##INFO=<ID=EVENT,Number=1,Type=String,Description="ID of event associated to breakend">
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
691 ##FORMAT=<ID=lSC,Number=1,Type=Integer,Description="Length of the longest soft clips supporting the BND">
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
692 ##FORMAT=<ID=nSC,Number=1,Type=Integer,Description="Number of supporting soft-clips\">
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
693 ##FORMAT=<ID=uRP,Number=1,Type=Integer,Description="Number of unmated read pairs nearby Soft-Clips">
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
694 ##FORMAT=<ID=levD_local,Number=1,Type=Float,Description="Levenstein distance between soft-clipped bases and the area around the original soft-clipped site">
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
695 ##FORMAT=<ID=levD_distl,Number=1,Type=Float,Description="Levenstein distance between the soft-clipped bases and mate location">
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
696 ##FORMAT=<ID=CTX,Number=1,Type=Integer,Description="Number of chromosomal translocations">
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
697 ##FORMAT=<ID=DEL,Number=1,Type=Integer,Description="Number of reads supporting Large Deletions">
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
698 ##FORMAT=<ID=INS,Number=1,Type=Integer,Description="Number of reads supporting Large insertions">
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
699 ##FORMAT=<ID=NOV_INS,Number=1,Type=Integer,Description="Number of reads supporting novel sequence insertion">
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
700 ##FORMAT=<ID=INV,Number=1,Type=Integer,Description="Number of reads supporting inversions">
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
701 ##FORMAT=<ID=sDEL,Number=1,Type=Integer,Description="Number of reads supporting novel sequence insertion">
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
702 ##INFO=<ID=NO_MATE_SC,Number=1,Type=Flag,Description="When there is no softclipping of the mate read location, an appromiate position is used">
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
703 ##FORMAT=<ID=GT,Number=1,Type=String,Description="Dummy value for maintaining VCF-Spec">
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
704 #CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\t$SAMPLE_NAME\n};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
705
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
706 print OUT $header;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
707 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
708
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
709
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
710 ###############################################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
711 sub bulk_work {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
712 my ($side, $line, $file) = @_;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
713 my $local_levD = 0;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
714 my $distl_levD = 0;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
715
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
716 #my @info = split(/\t/, $line);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
717 my @plus_Reads = split(/\t/, $line);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
718 $plus_Reads[7] =~ s/\n//g;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
719
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
720 #### softclip length and softclip size.
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
721 my $lSC = $plus_Reads[4];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
722 my $nSC = $plus_Reads[6];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
723
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
724
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
725 #Get all types of compatible reads
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
726 #Get improperly paired reads (@ max distance)
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
727
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
728 #### default value for left SIDE.
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
729 #If left-clip, then look downstream for match of softclipped reads to define a deletion, but look for DRPs upstream
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
730 my $sv_type = "SVTYPE=BND";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
731 my $start_local=0; my $end_local=0;my $target_local="";my $target_drp="";my $start_drp="";my $end_drp="";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
732 if ($side =~ /left/) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
733 $start_local = $plus_Reads[1]-$dist_To_Soft;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
734 $end_local = $plus_Reads[2];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
735 $start_drp = $plus_Reads[1];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
736 $end_drp = $plus_Reads[1]+$dist_To_Soft;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
737
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
738 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
739 else{
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
740 $start_local = $plus_Reads[1];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
741 $end_local = $plus_Reads[1]+$dist_To_Soft;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
742 $start_drp = $plus_Reads[1]-$dist_To_Soft;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
743 $end_drp = $plus_Reads[1];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
744 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
745
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
746 $target_local=join("", $plus_Reads[0], ":", $start_local, "-", $end_local);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
747 $target_drp=join("", $plus_Reads[0], ":", $start_drp, "-", $end_drp);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
748 my $num_unmapped_pairs="";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
749 if ($side =~ /right/) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
750 $num_unmapped_pairs=`samtools view $new_blacklist -q $MapQ -f8 -F 1536 -c $file $target_drp`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
751 } else {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
752 $num_unmapped_pairs=`samtools view $new_blacklist -q $MapQ -f24 -F 1536 -c $file $target_drp`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
753 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
754 if($verbose){print "samtools view $new_blacklist -q $MapQ -f24 -F 1536 -c $file $target_drp\n";}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
755
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
756 $num_unmapped_pairs=~s/\n//;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
757 if($verbose){print "NUM UNMAPPED PAIRS= $num_unmapped_pairs\n";}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
758 my $REF1_base = "";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
759 my $REF2_base = "";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
760 my $INFO_1 = "";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
761 my $INFO_2 = "";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
762 my $ALT_1 = "";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
763 my $ALT_2 = "";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
764 my $isize = 0;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
765 my $QUAL = "";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
766 my $FORMAT = "GT:";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
767
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
768 #### get 8 bit rand id
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
769 my $BND1_name = join "", map { ("a".."z")[rand 26] } 1..8;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
770 my $BND2_name = join "", map { ("a".."z")[rand 26] } 1..8;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
771 $BND1_name=join "_","BND",$BND1_name;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
772 $BND2_name=join "_","BND",$BND2_name;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
773
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
774 my $counts = {CTX => 0, DEL => 0, INS => 0, INV => 0, TDUP => 0, NOV_INS => 0 };
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
775 my $event_mate_info = {CTX => "", DEL => "", INS => "", INV => "", TDUP => "", NOV_INS => "" };
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
776
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
777 #### get mate pair info and counts per event
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
778 foreach my $e (sort keys %{$counts}) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
779 my $h = get_counts_n_info($e, $side, $MapQ, $file, $dist_To_Soft, $target_drp, $upper_limit, $lower_limit);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
780
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
781 $counts->{$e} = $h->{count};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
782 $event_mate_info->{$e} = $h->{info};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
783 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
784
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
785 my $max = 0;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
786 my $type = "UNKNOWN";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
787 my $nRP = 0;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
788 my $mate_info = "NA\tNA\tNA\tNA";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
789 my $summary = "GT:";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
790
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
791 #### find max count of events and set type, nRP and info to corresponding
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
792 #### max count event.
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
793 #### also create a summary string of all counts to be added to VCF file.
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
794 foreach my $e (sort keys %{$counts}){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
795 # if ($counts->{$e} >=i $max){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
796 if ($counts->{$e} > $max){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
797 $type = $e .",". $counts->{$e};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
798 $nRP = $counts->{$e};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
799
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
800 $max = $counts->{$e};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
801
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
802 if (length($event_mate_info->{$e})) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
803 $mate_info = $event_mate_info->{$e};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
804 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
805 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
806
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
807 $summary .= $e .",". $counts->{$e} .":";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
808 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
809 #print "done with Summaryi=$summary\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
810 #### remove last colon ":" from
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
811 $summary =~ s/:$//;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
812 if (($minRP > $max)&&(!$disable_RP_only )){return};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
813
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
814 #### Run Levenstein distance on softclip in target region to find out if its a small deletion/insetion
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
815 #### passing 1: clip_seq, 2: chr, 3: start, 4: end, 5: ref file.
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
816 my $levD = new LevD;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
817 ########################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
818 ########################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
819 ########################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
820
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
821 #### redefine start and end location for LevD calc.
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
822 # $start = $plus_Reads[1]-$dist_To_Soft;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
823 # $end = $plus_Reads[2];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
824 my $num_bases_to_loc=0;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
825 my $new_start=0;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
826 my $new_end=0;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
827 my $del_seq="";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
828 my $start = $start_local;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
829 my $end = $end_local;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
830 if ($lSC=~/NA/){$lSC=0}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
831
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
832 if ($side =~ /right/) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
833 $levD->search($plus_Reads[3], $plus_Reads[0], $start, $end, $INPUT_FASTA);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
834 $local_levD = sprintf("%.2f", $levD->{relative_edit_dist});
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
835 $num_bases_to_loc=$levD->{index};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
836 $new_start = $plus_Reads[2];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
837 if ($plus_Reads[2]=~/^[0-9]/){$new_end=$plus_Reads[2]+$lSC};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
838 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
839 else{
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
840 $levD->search($plus_Reads[3], $plus_Reads[0], $start, $end, $INPUT_FASTA);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
841 $local_levD = sprintf("%.2f", $levD->{relative_edit_dist});
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
842 $num_bases_to_loc=$levD->{index};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
843 if ($plus_Reads[2]=~/^[0-9]/){$new_start=$plus_Reads[2]-$lSC};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
844 $new_end = $plus_Reads[2];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
845 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
846 return if((!$new_start)||(!$new_end));
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
847 return if ($new_start<0);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
848 $del_seq=getBases($plus_Reads[0], $new_start,$new_end,$INPUT_FASTA);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
849 ##############################################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
850 # #If there is a match, where is the start position of the match?
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
851 #
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
852 ##############################################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
853
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
854
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
855 #if $plus_Reads[3] eq "NA", then it was found without soft-clipped reads
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
856 if($plus_Reads[3] !~ /NA/){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
857 if (($local_levD < $levD_local_threshold)) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
858 return if (!$sv_only);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
859 #### add value to summary to be written to vcf file.
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
860 $summary = "GT:sDel," . $plus_Reads[6];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
861 $type = "sDEL";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
862 ###########################################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
863 ##### Printing output
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
864
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
865 #########################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
866 ##### Get DNA info
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
867 #########################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
868 #$REF1_base = getSeq($plus_Reads[0], $plus_Reads[1], $INPUT_FASTA);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
869 $REF1_base = substr($del_seq, 0, 1);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
870
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
871 #### this is alt ref. for softclip its $plus_Reads[3]
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
872 $REF2_base = $del_seq;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
873 $QUAL = 1/($local_levD + 0.001);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
874 $QUAL = sprintf("%.2f",$QUAL);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
875 $isize = length($del_seq);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
876
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
877 #### svtype = none for sDEL
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
878 #### isize = length($info[3]);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
879 #### nRP = NA
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
880 #### mate_id = NA
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
881 #### CTX,:DEL,:....sDEL,##
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
882 $INFO_1=join(";", "SVTYPE=NA", "EVENT=$type", "ISIZE=$isize");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
883
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
884 #Add Sample infomration
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
885 my $FORMAT="GT:sDEL";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
886 $FORMAT .= ":lSC:nSC:uRP:levD_local";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
887 my $SAMPLE= "0/1:";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
888 $SAMPLE .= "$plus_Reads[6]:$lSC:$nSC:$num_unmapped_pairs:$local_levD";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
889
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
890 #### remove any white spaces.
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
891 $INFO_1=~s/\s//g;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
892 $INFO_2=~s/\s//g;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
893
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
894 $BND1_name =~ s/^BND/LEVD/;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
895 # If left, then the start position is plus_Reads[1]-isize
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
896 my $start_pos=0;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
897 #Make sure Ref1 and Ref2 bases are different
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
898 if($REF2_base eq $REF1_base){$REF1_base="NA"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
899 if($side=~/left/){$start_pos=$plus_Reads[1]-$isize}else{$start_pos=$plus_Reads[1]};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
900 my $var=join("\t", $plus_Reads[0], $start_pos, $BND1_name, $REF2_base, $REF1_base, $QUAL, "PASS", $INFO_1,$FORMAT,$SAMPLE);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
901 return $var;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
902 #print OUT join ("\t", $plus_Reads[0], $start_pos, $BND1_name, $REF2_base, $REF1_base, $QUAL, "PASS", $INFO_1,$FORMAT,$SAMPLE, "\n");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
903 # return;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
904 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
905 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
906
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
907 #### Otherwise, look for DRP mate info
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
908 #if($nRP=~/NA/){print "MATE_INFO=$mate_info\tSide=$side\tline=$line\n";}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
909 my @mate_info_arr = split(/\t/, $mate_info);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
910 $nRP = $mate_info_arr[3];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
911 my $mate_chr=$mate_info_arr[0];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
912
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
913 if((! defined $nRP) || ($nRP =~ /na/i) || ($mate_chr =~ /NA/) ){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
914 #PRINT UNKNOWN
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
915 return if ($nRP =~ /na/i);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
916 #print "There is an unknown\nNRP=$nRP Mate_CHR=$mate_chr minRP=$minRP\n";die;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
917 $summary .= ":unknown," . $plus_Reads[6];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
918 $type = "unknown";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
919 $REF1_base = getSeq($plus_Reads[0], $plus_Reads[1], $INPUT_FASTA);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
920 $REF2_base = $plus_Reads[3];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
921 $BND1_name =~ s/^BND/UNKNOWN/;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
922 $QUAL = 1/($local_levD + 0.001);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
923 $QUAL = sprintf("%.2f",$QUAL);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
924 $INFO_1=join(";", "SVTYPE=unknown", "EVENT=unknown", "ISIZE=unknown");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
925 #Add Sample infomration
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
926 my $FORMAT="GT:sDEL";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
927 $FORMAT .= ":lSC:nSC:uRP:levD_local";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
928 my $SAMPLE = "0/1:";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
929 $SAMPLE .= "$plus_Reads[6]:$lSC:$nSC:$num_unmapped_pairs:$local_levD";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
930 #### remove any white spaces.
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
931 $INFO_1=~s/\s//g;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
932 #print join ("\t", $plus_Reads[0], $plus_Reads[1], $REF2_base, $REF1_base, $QUAL, "PASS", $INFO_1,$FORMAT,$SAMPLE, "\n");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
933
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
934 #print OUT join ("\t", $plus_Reads[0], $plus_Reads[1], $BND1_name, $REF1_base, $REF2_base, $QUAL, "PASS", $INFO_1,$FORMAT,$SAMPLE, "\n");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
935 my $var=join("\t", $plus_Reads[0], $plus_Reads[1], $BND1_name, $REF1_base, $REF2_base, $QUAL, "PASS", $INFO_1,$FORMAT,$SAMPLE);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
936 return $var;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
937
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
938 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
939
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
940 #### end if there is no mate info or nRP+uRP<minRP
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
941 return if (($nRP<$minRP)&&($unmated_pairs > ($num_unmapped_pairs+$nRP)));
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
942
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
943 ##################################################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
944 # Find out if mates have nearby soft-clips (to refine the breakpoints)
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
945 ##################################################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
946 #Look for evidence of soft-clipping near mate
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
947 my @mate_soft_arr = ();
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
948 my $mate_start = 0;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
949 my $mate_soft = "";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
950
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
951 @mate_info_arr = split(/\t/, $mate_info);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
952
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
953 #### mate start and end locations.
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
954 my $filename = $right_sc;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
955
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
956 $start = $mate_info_arr[1] - $dist_To_Soft;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
957 $end = $mate_info_arr[1];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
958
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
959 if ($side =~ /right/) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
960 $start = $mate_info_arr[2];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
961 $end = $mate_info_arr[2] + $dist_To_Soft;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
962
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
963 $filename = $left_sc;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
964 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
965
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
966 #### add levenstein distance to Summary
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
967 #print "Calc distal Levd\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
968 $levD->search(rev_comp($plus_Reads[3]), $mate_info_arr[0], $start, $end, $INPUT_FASTA);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
969 $distl_levD = sprintf("%.2f", $levD->{relative_edit_dist});
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
970 $distl_levD = "NA" if($plus_Reads[3] =~ /NA/);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
971 #If there is no softclips to string match, then give 0 as quality value
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
972 if ($plus_Reads[3] !~ /NA/){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
973 $QUAL=1/($distl_levD + 0.001);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
974 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
975 else {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
976 $QUAL=0;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
977 };
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
978 $QUAL=sprintf("%.2f",$QUAL);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
979 #### looking for softclips to refine break point
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
980 #### if left look in right and vice-versa.
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
981 $cmd = qq{echo -e "$mate_info_arr[0]\t$start\t$end"};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
982 $cmd .= qq{ | awk -F'\t' 'NR==3' | intersectBed -a stdin -b $filename | head -1};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
983
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
984 $mate_soft = `$cmd`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
985
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
986 $mate_soft =~ s/\n//g;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
987 @mate_soft_arr = split(/\s/, $mate_soft);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
988 my $NO_MATE_SC="";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
989 if(@mate_soft_arr){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
990 $mate_chr = $mate_soft_arr[0];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
991 $mate_start = $mate_soft_arr[1];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
992 $NO_MATE_SC="APPROXIMATE";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
993
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
994 } else{
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
995 @mate_info_arr = split(/\s/,$mate_info);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
996 $mate_chr = $mate_info_arr[0];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
997 $mate_start = $mate_info_arr[1];
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
998 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
999
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1000 #end if there is no mate info
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1001 return if ($mate_chr eq "");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1002 #end if there is no mate info and !disable_RP_only
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1003 return if (($lSC =~/NA/)&&(!$disable_RP_only));
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1004
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1005
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1006 ###########################################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1007 ##### Printing output
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1008
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1009 #########################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1010 # Get DNA info
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1011 #########################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1012 #print "PLUS_READS=$plus_Reads[0],$plus_Reads[1]\nMATE=$mate_chr,$mate_start,$INPUT_FASTA\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1013 $REF1_base = getSeq($plus_Reads[0], $plus_Reads[1], $INPUT_FASTA);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1014
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1015 ### this is alt ref. for softclip its $plus_Reads[3]
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1016 $REF2_base = getSeq($mate_chr, $mate_start, $INPUT_FASTA);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1017
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1018 #########################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1019 # print in VCF format
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1020 #########################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1021
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1022 #### abs value to account for left and right reads.
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1023 $isize = abs($plus_Reads[1]-$mate_start);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1024
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1025 my $event_type=$type;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1026 $event_type=~ s/,|[0-9]//g;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1027 $INFO_1=join(";", "$sv_type", "EVENT=$event_type", "ISIZE=$isize","MATE_ID=$BND2_name");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1028 $INFO_2=join(";", "$sv_type", "EVENT=$event_type", "ISIZE=$isize","MATE_ID=$BND1_name");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1029
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1030 #### remove any white spaces.
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1031 #### ask: did you mean to remove space from ends? eg. trim()
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1032 $INFO_1=~s/\s//g;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1033 $INFO_2=~s/\s//g;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1034
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1035 $FORMAT=$summary;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1036 $FORMAT=~ s/,|[0-9]//g;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1037 $FORMAT .= ":lSC:nSC:uRP:distl_levD";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1038 if($NO_MATE_SC){$INFO_2 .= ":NO_MATE_SC"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1039 my $SAMPLE="0/1:";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1040 $SAMPLE .=$summary;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1041 # if($NO_MATE_SC){$SAMPLE.= ":$NO_MATE_SC"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1042
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1043 $SAMPLE=~s/[A-Z|,|_]//g;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1044 my $MATE_SAMPLE=$SAMPLE;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1045 $SAMPLE .= ":$lSC:$nSC:$num_unmapped_pairs:$distl_levD";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1046 $MATE_SAMPLE .=":NA:NA:NA:NA";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1047 $SAMPLE=~s/::/:/g;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1048 $MATE_SAMPLE=~s/::/:/g;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1049
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1050 if($type !~ /INV/){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1051 $ALT_1 = join("","]",$mate_chr,":",$mate_start,"]",$REF1_base);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1052 $ALT_2 = join("",$REF2_base,"[",$plus_Reads[0],":",$plus_Reads[1],"[");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1053 # 2 321682 bnd_V T ]13:123456]T 6 PASS SVTYPE=BND
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1054 # 13 123456 bnd_U C C[2:321682[ 6 PASS SVTYPE=BND
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1055 } else {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1056 $ALT_1 = join("", "]", $plus_Reads[0], ":", $plus_Reads[1], "]", $REF2_base);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1057 $ALT_2 = join("", $REF1_base, "[", $mate_chr, ":", $mate_start, "[");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1058 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1059
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1060 if(($mate_chr) && ($plus_Reads[0])){
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1061 # print OUT join ("\t", $plus_Reads[0], $plus_Reads[1], $BND1_name, $REF1_base, $ALT_1, $QUAL,"PASS", $INFO_1, $FORMAT,$SAMPLE,"\n");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1062 # print OUT join ("\t", $mate_chr, $mate_start, $BND2_name, $REF2_base, $ALT_2, $QUAL, "PASS", $INFO_2, $FORMAT,$MATE_SAMPLE,"\n");
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1063 my $var=join ("\t", $plus_Reads[0], $plus_Reads[1], $BND1_name, $REF1_base, $ALT_1, $QUAL,"PASS", $INFO_1, $FORMAT,$SAMPLE);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1064 return $var;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1065 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1066 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1067
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1068 ###############################################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1069 ###############################################################################
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1070 sub get_counts_n_info {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1071 my ($event, $side, $mapQ, $file, $dist, $target, $upL, $lwL) = @_;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1072
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1073 my $mate_info = "";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1074 my $cmd = "";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1075
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1076 if ($event =~ /^CTX$/i) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1077 #print "CTX side $side\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1078 if ($side =~ /right/i) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1079 $cmd = qq{samtools view $new_blacklist -q $mapQ -f 32 -F 1536 $file $target};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1080 $cmd .= qq{ | perl -ane 'if(\$F[6] ne "="){\$end=\$F[7]+1; print join ("\\t",\$F[6],\$F[7],\$end,"\\n")}'};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1081 $cmd .= qq{ | sortBed | mergeBed -d $dist -n | sort -k4nr | head -1};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1082 if($verbose){print "$cmd\n"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1083 $mate_info=`$cmd`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1084 } else {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1085 $cmd = qq{ samtools view $new_blacklist -q $mapQ -f 16 -F 1536 $file $target};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1086 $cmd .= qq{ | perl -ane 'if(\$F[6] ne "="){\$end=\$F[7]+1; print join ("\\t",\$F[6],\$F[7],\$end,"\\n")}'};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1087 $cmd .= qq{ | sortBed | mergeBed -d $dist -n | sort -k4nr | head -1};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1088 if($verbose){print "$cmd\n"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1089 $mate_info=`$cmd`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1090 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1091 } elsif ($event =~ /^DEL$/i) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1092 #print "DEL side $side\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1093 if ($side =~ /right/i) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1094 $cmd = qq{samtools view $new_blacklist -q $mapQ -f 32 -F 1552 $file $target};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1095 $cmd .= qq{ | awk '{OFS="\\t"}{if((\$7 ~ /=/)&&(\$9>$upL)){end=\$8+1;print \$3,\$8,end}}'};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1096 $cmd .= qq{ | sortBed | mergeBed -d $dist -n | sort -k4nr | head -1};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1097 if($verbose){print "$cmd\n"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1098 $mate_info=`$cmd`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1099 } else {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1100 $cmd = qq{samtools view $new_blacklist -q $mapQ -F 1568 -f 16 $file $target};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1101 $cmd .= qq{ | awk '{OFS="\\t"} {if((\$7 ~ /=/)&&(\$9<-$upL)){end=\$8+1;print \$3,\$8,end}}'};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1102 $cmd .= qq{ | sortBed | mergeBed -d $dist_To_Soft -n | sort -k4nr | head -1};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1103 if($verbose){print "$cmd\n"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1104
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1105 $mate_info=`$cmd`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1106 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1107 } elsif ($event =~ /^INS$/i) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1108 #print "INS side $side\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1109 if ($side =~ /right/i) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1110 $cmd = qq{samtools view $new_blacklist -q $mapQ -f 32 -F 1552 $file $target};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1111 $cmd .= qq{ | awk '{OFS="\\t"}{if((\$7 ~ /=/)&&(\$9<$lwL && \$9 > 0 )){end=\$8+1;print \$3,\$8,end}}'};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1112 $cmd .= qq{ | sortBed | mergeBed -d $dist -n | sort -k4nr | head -1};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1113 if($verbose){print "$cmd\n"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1114 $mate_info = `$cmd`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1115 } else {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1116 $cmd = qq {samtools view $new_blacklist -q $mapQ -f 16 -F 1568 $file $target};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1117 $cmd .= qq{ | awk '{OFS="\\t"}{if((\$7 ~ /=/)&&(\$9>-$lwL && \$9 < 0 )){end=\$8+1;print \$3,\$8,end}}'};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1118 $cmd .= qq{ | sortBed | mergeBed -d $dist -n | sort -k4nr | head -1};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1119 if($verbose){print "$cmd\n"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1120 $mate_info = `$cmd`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1121 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1122 } elsif ($event =~ /^INV$/i) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1123 #print "INV side $side\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1124 if ($side =~ /right/i) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1125 $cmd = qq{samtools view $new_blacklist -q $mapQ -F 1596 $file $target};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1126 $cmd .= qq{ | awk '{OFS="\\t"}{if((\$7 ~ /=/)){end=\$8+1;print \$3,\$8,end}}'};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1127 $cmd .= qq{ | sortBed | mergeBed -d $dist_To_Soft -n | sort -k4nr | head -1};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1128 if($verbose){print "$cmd\n"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1129 $mate_info = `$cmd`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1130 } else {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1131 $cmd = qq{samtools view $new_blacklist -q $mapQ -f 48 -F 1548 $file $target};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1132 $cmd .= qq{ | awk '{OFS="\\t"}{if((\$7 ~ /=/)){end=\$8+1;print \$3,\$8,end}}'};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1133 $cmd .= qq{ | sortBed | mergeBed -d $dist -n | sort -k4nr | head -1};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1134 if($verbose){print "$cmd\n"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1135 $mate_info = `$cmd`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1136 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1137 } elsif ($event =~ /^TDUP$/i) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1138 #print "TDUP side $side\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1139 if ($side =~ /right/i) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1140 $cmd = qq{samtools view $new_blacklist -q $mapQ -f 32 -F 1552 $file $target};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1141 # $cmd .= qq{ | awk '{OFS="\\t"}{if((\$7 ~ /=/)&&(\$9>$upL)){end=\$8+1;print \$3,\$8,end}}'};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1142 $cmd .= qq{ | awk '{OFS="\\t"}{if((\$7 ~ /=/)&&(\$4>\$8)&&(\$9<0)){end=\$8+1;print \$3,\$8,end}}'};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1143 $cmd .= qq{ | sortBed | mergeBed -d $dist -n | sort -k4nr | head -1};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1144 if($verbose){print "$cmd\n"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1145 $mate_info = `$cmd`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1146 } else {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1147 $cmd = qq{samtools view $new_blacklist -q $mapQ -f 16 -F 1568 $file $target};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1148 # $cmd .= qq{ | awk '{OFS="\\t"}{if((\$7 ~ /=/)&&(\$9<-$upL )){end=\$8+1;print \$3,\$8,end}}'};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1149 $cmd .= qq{ | awk '{OFS="\\t"}{if((\$7 ~ /=/)&&(\$4<\$8)&&(\$9>0)){end=\$8+1;print \$3,\$8,end}}'};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1150 $cmd .= qq{ | sortBed | mergeBed -d $dist -n | sort -k4nr | head -1};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1151 if($verbose){print "$cmd\n"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1152 $mate_info = `$cmd`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1153 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1154 } elsif ($event =~ /^NOV_INS$/i) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1155 #print "NOV_INS side $side\n";
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1156 if ($side =~ /right/i) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1157 $cmd = qq{samtools view $new_blacklist -q $mapQ -f 8 -F 1552 $file $target};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1158 $cmd .= qq{ | awk '{OFS="\\t"}{end=\$8+1;print \$3,\$8,end}'};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1159 $cmd .= qq{ | sortBed | mergeBed -d $dist -n | sort -k4nr | head -1};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1160 if($verbose){print "$cmd\n"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1161 $mate_info = `$cmd`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1162 } else {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1163 $cmd = qq{samtools view $new_blacklist -q $mapQ -f 24 -F 1536 $file $target};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1164 $cmd .= qq{ | awk '{OFS="\\t"}{end=\$8+1;print \$3,\$8,end}'};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1165 $cmd .= qq{ | sortBed | mergeBed -d $dist -n | sort -k4nr | head -1};
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1166 if($verbose){print "$cmd\n"}
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1167 $mate_info = `$cmd`;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1168 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1169 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1170
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1171 $mate_info=~s/\n//g;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1172 my @tmp=split(/\t/, $mate_info);
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1173
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1174 my $counts = 0;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1175
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1176 if (defined $tmp[3]) {
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1177 $tmp[3] =~ s/\n//g;
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1178
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1179 $counts = $tmp[3] if (length($tmp[3]));
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1180 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1181 return ({count=>$counts, info=>$mate_info});
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1182 }
e3609c8714fb Uploaded
plus91-technologies-pvt-ltd
parents:
diff changeset
1183