Mercurial > repos > jjohnson > crest
diff StructVar.pm @ 0:acc8d8bfeb9a
Uploaded
| author | jjohnson |
|---|---|
| date | Wed, 08 Feb 2012 16:59:24 -0500 |
| parents | |
| children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/StructVar.pm Wed Feb 08 16:59:24 2012 -0500 @@ -0,0 +1,1240 @@ +package StructVar; + +# $Id: StructVar.pm, v 1.00 2010/09/24 jianmin.wang Exp $ + +use strict; +use Bio::DB::Sam; +use Carp; +use Data::Dumper; +use English; +use Bio::SearchIO; +use GeneModel; +use Cwd; +use SVExtTools; +use List::MoreUtils qw/ uniq /; +use SVUtil qw( prepare_reads read_fa_file get_sclip_reads rev_comp); +use SCValidator qw(LEFT_CLIP RIGHT_CLIP); + +use constant UNK => 0; +use constant CTX => 1; +use constant INV => 2; +use constant INS => 3; +use constant DEL => 4; +use constant ITX => 5; + +# type of SVs and thier string +my %type2str = ( + 0 => 'UNKNOWN', # unknow type, not single event + 1 => 'CTX', # cross chromosome translocation + 2 => 'INV', # inversion, requires both breakpoints available + 3 => 'INS', # (tandem) insertion + 4 => 'DEL', # deletion + 5 => 'ITX', # within chromosome translocation +); +# class static variables +my $sam_d; +my $sam_g; +my $gm; +my $assembler; +my $aligner; +my $mapper; +my @filters = (\&_germline_sclip_filter, \&_type_distance_filter, \&_germline_indel_filter); +my $RNASeq; +my $read_len; +my $genome; +my $tr_max_indel_size; #tandem repeat mediated indel size +my $tr_min_size; +my $tr_max_size; #tandem repeat max_size of single repeat +my $tr_min_num; #tandem repeat minimum number of occurence + +my $max_bp_dist = 15; +my $germline_seq_width = 100; +my $germline_search_width = 50; +my $min_percent_cons_of_read = 0.75; + +my %default_filters = ( + 'germline_sclip' => \&_germline_sclip_filter, + 'type_distance' => \&_type_distance_filter, + 'germline_indel' => \&_germline_indel_filter, +# 'mapping_artifact' => \&_mapping_artifact_filter, + 'mapping_quality' => \&_mapping_quality_filter, + 'tandem_repeat' => \&_tandem_repeat_filter, +); + +sub new { + my $class = shift; + my %param = @_; + my $self = {}; + $self->{first_bp} = {}; + $self->{second_bp} = {}; + $self->{first_bp} = $param{-FIRST_BP} if($param{-FIRST_BP}); + $self->{second_bp} = $param{-SECOND_BP} if($param{-SECOND_BP}); + $self->{consensus} = $param{-CONSENSUS} ? $param{-CONSENSUS} : ""; + $self->{type} = ""; + bless $self, ref($class) || $class; + return $self; +} + +sub add_filter { + my $self = shift; + my $name = shift; + my $f = shift; + croak "you must specify a filter" if(!$f); + croak "the filter must be a subroutine" if(ref($f) ne 'CODE'); + $default_filters{$name} = $f; +} + +sub add_RNASeq_filter { + $default_filters{'RNASeq_strand'} = \&_RNASeq_strand_filter; + $default_filters{'RNASeq_INS'} = \&_RNASeq_INS_filter; +} + +sub remove_filter { + my $self = shift; + my $name = shift; + if($default_filters{$name}){ + delete $default_filters{$name}; + } + else { + print STDERR "Filter $name is not in filter list"; + } +} + +sub tr_max_indel_size { #tandem repeat mediated indel size + my $self = shift; + my $value = shift; + $tr_max_indel_size = $value if($value); + return $tr_max_indel_size; +} + +sub min_percent_cons_of_read { + my $self= shift; + my $value = shift; + $min_percent_cons_of_read = $value if($value); + return $min_percent_cons_of_read; +} + +sub tr_min_size { + my $self = shift; + my $value = shift; + $tr_min_size = $value if($value); + return $tr_min_size; +} +sub germline_seq_width { + my $self = shift; + my $value = shift; + $germline_seq_width = $value if($value); + return $germline_seq_width; +} + +sub germline_search_width { + my $self = shift; + my $value = shift; + $germline_search_width = $value if($value); + return $germline_search_width; +} + +sub max_bp_dist { + my $self = shift; + my $value = shift; + $max_bp_dist = $value if($value); + return $max_bp_dist; +} + +sub tr_max_size { #tandem repeat max_size of single repeat + my $self = shift; + my $value = shift; + $tr_max_size = $value if($value); + return $tr_max_size; +} + +sub tr_min_num { #tandem repeat minimum number of occurence + my $self = shift; + my $value = shift; + $tr_min_num = $value if($value); + return $tr_min_num; +} +sub genome { + my $self = shift; + my $value = shift; + $genome = $value if($value); + return $genome; +} +sub sam_d { + my $self = shift; + my $value = shift; + $sam_d = $value if($value); + return $sam_d; +} + +sub sam_g { + my $self = shift; + my $value = shift; + $sam_g = $value if($value); + return $sam_g; +} + +sub assembler { + my $self = shift; + my $value = shift; + $assembler = $value if($value); + return $assembler; +} + +sub aligner { + my $self = shift; + my $value = shift; + $aligner = $value if($value); + return $aligner; +} + +sub RNASeq { + my $self = shift; + my $value = shift; + $RNASeq = 1 if($value); + return $RNASeq; +} + +sub mapper { + my $self = shift; + my $value = shift; + $mapper = $value if($value); + return $mapper; +} + +sub gene_model { + my $self = shift; + my $value = shift; + $gm = $value if($value); + return $gm; +} + +sub read_len { + my $self = shift; + my $value = shift; + $read_len = $value if($value); + return $read_len; + +} +sub first_bp { + my $self = shift; + my %param = @_; + if(%param) { + $self->{first_bp}{left_chr} = $param{-LEFT_CHR}; + $self->{first_bp}{left_pos} = $param{-LEFT_POS}; + $self->{first_bp}{left_ort} = $param{-LEFT_ORT}; + $self->{first_bp}{right_chr} = $param{-RIGHT_CHR}; + $self->{first_bp}{right_pos} = $param{-RIGHT_POS}; + $self->{first_bp}{right_ort} = $param{-RIGHT_ORT}; + $self->{first_bp}{sc_reads} = $param{-SC_READS}; + $self->{first_bp}{pos} = $param{-POS}; + $self->{first_bp}{cover} = $param{-COVER}; + $self->{first_bp}{sc_seq} = $param{-SC_SEQ}; + $self->{first_bp}{chr} = $param{-CHR}; + } + return $self->{first_bp}; +} + +sub second_bp { + my $self = shift; + my %param = @_; + if(%param) { + $self->{second_bp}{left_chr} = $param{-LEFT_CHR}; + $self->{second_bp}{left_pos} = $param{-LEFT_POS}; + $self->{second_bp}{left_ort} = $param{-LEFT_ORT}; + $self->{second_bp}{right_chr} = $param{-RIGHT_CHR}; + $self->{second_bp}{right_pos} = $param{-RIGHT_POS}; + $self->{second_bp}{right_ort} = $param{-RIGHT_ORT}; + $self->{second_bp}{sc_reads} = $param{-SC_READS}; + $self->{second_bp}{pos} = $param{-POS}; + $self->{second_bp}{cover} = $param{-COVER}; + $self->{second_bp}{sc_seq} = $param{-SC_SEQ}; + $self->{second_bp}{chr} = $param{-CHR}; + } + return $self->{second_bp}; +} + +sub type { + my $self = shift; + my ($bp1, $bp2) = ($self->{first_bp}, $self->{second_bp}); + if($bp1->{left_chr} ne $bp1->{right_chr}) { + return CTX if($bp1->{left_ort} eq $bp2->{left_ort} && + $bp1->{right_ort} eq $bp2->{right_ort}); + + return UNK; + } + if( $bp1->{left_ort} eq $bp1->{right_ort} && + $bp2->{left_ort} eq $bp2->{right_ort} ) { # Insertion or deletion + return DEL if( $bp1->{left_pos} <= $bp1->{right_pos} && + $bp2->{left_pos} <= $bp2->{right_pos}); + return INS if( $bp1->{left_pos} >= $bp1->{right_pos} && + $bp2->{left_pos} >= $bp2->{right_pos}); + return UNK; + } + if( $bp1->{left_ort} ne $bp1->{right_ort} && + $bp2->{left_ort} ne $bp2->{right_ort} ) { + return INV if( ($bp1->{left_ort} eq '+' && $bp2->{left_ort} eq '-') + || ($bp2->{left_ort} eq '+' && $bp1->{left_ort} eq '-')); + return ITX if( ($bp1->{left_ort} eq '+' && $bp2->{left_ort} eq '+') + || ($bp1->{left_ort} eq '-' && $bp2->{left_ort} eq '-') ); + return UNK; + } + return UNK; +} + +sub to_string { + my $self = shift; + my $type = $self->type; + my ($bp1, $bp2) = ($self->{first_bp}, $self->{second_bp}); + my ($left_len, $right_len) = (0, 0); + if(exists $bp1->{left_len}) { + ($left_len, $right_len ) = ($bp1->{left_len}, $bp1->{right_len}); + } + else { + $left_len = length($bp1->{sc_seq}) if($bp1->{sc_seq}); + $right_len = length($bp2->{sc_seq}) if($bp2->{sc_seq}); + } + if($type == INV && $bp1->{pos} > $bp2->{pos}) { + ($bp1, $bp2) = ($bp2, $bp1); + ($left_len, $right_len) = ($right_len, $left_len); + } + my ($cover1, $cover2) = ($bp1->{cover}, $bp2->{cover}); + + my ($chrA, $chrB, $ortA, $ortB) = ($bp1->{left_chr}, $bp1->{right_chr}, + $bp1->{left_ort}, $bp1->{right_ort}); + my ($posA, $posB, $countA, $countB) = ($bp1->{pos}, $bp2->{pos}, $bp1->{sc_reads}, $bp2->{sc_reads}); + if($bp1->{chr} eq $bp1->{right_chr} && $bp1->{pos} == $bp1->{right_pos}) { + $posA = $bp2->{pos}; + $posB = $bp1->{pos}; + $countA = $bp2->{sc_reads}; + $countB = $bp1->{sc_reads}; + ($cover1, $cover2) = ($cover2, $cover1); + ($left_len, $right_len) = ($right_len, $left_len); + } + + my $rtn = join("\t", ($chrA, $posA, $ortA, $countA, $chrB, $posB, $ortB, $countB, + $type2str{$type}, $cover1, $cover2, $left_len, $right_len, $self->get_statistics(100), + $self->{c_start}, $self->{start_chr}, $self->{t_start}, + $self->{c_end}, $self->{end_chr}, $self->{t_end}, $self->{consensus})); + if($self->{type} == INV) { + $rtn = join("\t", ($rtn, $self->{c_start2}, $self->{start_chr2}, $self->{t_start2}, + $self->{c_end2}, $self->{end_chr2}, $self->{t_end}, $self->{consensus2})); + } + return $rtn; +} + +sub set_consensus { + my $self = shift; + my ($validator, $paired ) = @_; + if(!$assembler) { + print STDERR "No assembler set, no consensus generated\n"; + return; + } + my %all_reads; + my ($bp1, $bp2) = ($self->{first_bp}, $self->{second_bp}); + my $N = 0; + + my @s_reads1 = $self->_get_reads($sam_d, $bp1); + my @s_reads2 = $self->_get_reads($sam_d, $bp2); + + if($self->{type} == INV) { + $self->{consensus} = _generate_consensus(\@s_reads1); + $self->{consensus2} = _generate_consensus(\@s_reads2); + } + else { + push @s_reads1, @s_reads2; + $self->{consensus} = _generate_consensus(\@s_reads1); + } +} + +sub _generate_consensus { + my $s_reads = shift; + my $N = scalar(@{$s_reads}); + return $s_reads->[0]{seq} if($N == 1); + my $fa_name = prepare_reads($s_reads); + my ($contig_file, $sclip_count) = $assembler->run($fa_name); + if(!$contig_file || -s $contig_file == 0) { + system ("rm $fa_name"); system("rm $fa_name*"); + return; + } + my $n = 0; + my $contig_name; + foreach my $c (keys %{$sclip_count}) { + if($sclip_count->{$c} > $n) { + $n = $sclip_count->{$c}; + $contig_name = $c; + } + } + return if($N * 8 > $n * 10); + my $contig_seqs = read_fa_file($contig_file); + return $contig_seqs->{$contig_name}; +} + +sub _get_reads{ + my $self = shift; + my ($sam, $bp) = @_; + + my ($chr, $pos) = ($bp->{chr}, $bp->{pos}); + my ($start, $end) = ($pos, $pos); + if($bp->{all_pos}) { + my @tmp = @{$bp->{all_pos}}; + @tmp = sort {$a <=> $b} @tmp; + ($start, $end) = ($tmp[0], $tmp[$#tmp]); + } + my %reads; + return if(scalar @{$bp->{reads}} == 0); + foreach my $r (@{$bp->{reads}}){ + $reads{$r} = 1; + } + + #my ($sam, $chr, $start, $end, $reads) = + # ($args{-SAM}, $args{-CHR}, $args{-START}, $args{-END}, $args{-READS}) ; + + my @rtn; + my $range = $chr; + $range = $chr . ":" . $start . "-" . $end; + my($s, $e) = ($start, $end); + $sam->fetch($range, sub { + my $a = shift; + return unless (exists $reads{$a->qname}); + return unless ($a->start && $a->end); + return unless (($a->start >= $start && $a->start <= $end) + || ($a->end >= $start && $a->end <= $end)); + $s = $a->start if($s > $a->start); + $e = $a->end if($e < $a->end); + my $qscore = $a->qscore; + delete $reads{$a->qname}; + push @rtn, { + name => $a->qname, + seq => $a->query->dna, + qual => $qscore, + }; + } ); +# $bp->{range} = [$s, $e]; + if($self->{type} == INV) { + if($start == $bp->{left_pos} || $end == $bp->{left_pos}) { + $bp->{left_range} = [$s, $e]; + } + else { + $bp->{right_range} = [$s, $e]; + } + return @rtn; + } + if($chr eq $self->{first_bp}{left_chr} && ($start == $self->{first_bp}{left_pos} || $end == $self->{first_bp}{left_pos})){ + $self->{first_bp}{left_range} = [$s, $e]; + } + if($chr eq $self->{first_bp}{right_chr} && ($start == $self->{first_bp}{right_pos} || $end == $self->{first_bp}{right_pos})){ + $self->{first_bp}{right_range} = [$s, $e]; + } + + return @rtn; +} + +my $start_dir; +BEGIN { + $start_dir = getcwd; +} + +sub _cat_genes { + my ($ort, @genes) = @_; + my $rtn; + my @names; + foreach my $g (@genes) { + push @names, $g->val->name; + } + @names = uniq @names; + $rtn = join(";", @names); + $rtn = $rtn . "($ort)" if($rtn); + return $rtn; +} + +sub get_genes { + my $self = shift; + return ['NA', 'NA', 'NA'] unless $gm; + + my ($gene1, $gene2); + my $dist; + my ($bp1, $bp2) = ($self->{first_bp}, $self->{second_bp}); + my ($ort1, $ort2) = ($bp1->{left_ort}, $bp1->{right_ort}); + my ($chr1, $chr2) = ($bp1->{left_chr}, $bp1->{right_chr}); + my ($pos1, $pos2) = ($bp1->{pos}, $bp2->{pos}); + ($pos1, $pos2) = ($pos2, $pos1) if($bp1->{right_pos} == $pos1); + my $range1 = $ort1 eq "+" ? [$pos1 - 5, $pos1] : [$pos1, $pos1 + 5]; + my $range2 = $ort2 eq "+" ? [$pos2, $pos2 + 5] : [$pos2 - 5, $pos2]; + + my ($tmp1, $tmp2) = ($chr1 =~ m/chr/ ? $chr1 : "chr" . $chr1, + $chr2 =~ m/chr/ ? $chr2 : "chr" . $chr2); + my ($r_tree1, $f_tree1) = ($gm->sub_model($tmp1, "-"), $gm->sub_model($tmp1, "+")); + my ($r_tree2, $f_tree2) = ($gm->sub_model($tmp2, "-"), $gm->sub_model($tmp2, "+")); + + my ($g_ort1, $g_ort2); + my @genes; + @genes = $r_tree1->intersect($range1) if($r_tree1); + $gene1 = _cat_genes("-", @genes) if(scalar @genes > 0 ); + undef @genes; + @genes = $f_tree1->intersect($range1) if($f_tree1); + if(scalar @genes > 0) { + my $tmp = _cat_genes("+", @genes); + if($gene1) { + $gene1 .= ";" . $tmp; + } + else { + $gene1 = $tmp; + } + } + undef @genes; + @genes = $r_tree2->intersect($range2) if($r_tree2); + $gene2 = _cat_genes("-", @genes) if(scalar @genes > 0 ); + undef @genes; + @genes = $f_tree2->intersect($range2) if($f_tree2);; + if(scalar @genes > 0) { + my $tmp = _cat_genes("+", @genes); + if($gene2) { + $gene2 .= ";" . $tmp; + } + else { + $gene2 = $tmp; + } + } + + $gene1 = 'NA' unless $gene1; + $gene2 = 'NA' unless $gene2; + my $type = $self->{type}; + if( $type == INS or $type == DEL ) { + ($pos1, $pos2) = ($pos2, $pos1) if($pos1 > $pos2); + @genes = $r_tree1->intersect([$pos1 - 5, $pos2 + 5]); + if(scalar @genes >= 2) { + $dist = scalar(@genes) - 2; + $dist .= "(-)"; + } + @genes = $f_tree1->intersect([$pos1 - 5, $pos2 + 5]); + if(scalar @genes >= 2) { + if($dist) { + $dist .= ":" . (scalar (@genes) - 2); + } + else { + $dist = scalar (@genes) - 2; + } + $dist .= "(+)"; + } + } + $dist = 'NA' unless $dist; + return [$gene1, $gene2, $dist]; +} + +sub to_full_string { + my $self = shift; + my $type = $self->{type}; + my ($bp1, $bp2) = ($self->{first_bp}, $self->{second_bp}); + return join("\t", ( $bp1->{left_chr}, $bp1->{left_pos}, $bp1->{left_ort}, + $bp1->{right_chr}, $bp1->{right_pos}, $bp1->{right_ort}, $bp1->{sc_reads}, + $bp2->{left_chr}, $bp2->{left_pos}, $bp2->{left_ort}, + $bp2->{right_chr}, $bp2->{right_pos}, $bp2->{right_ort}, $bp2->{sc_reads}, + $type2str{$type})); +} + +sub filter { + my $self = shift; + print STDERR "SV filter starting....\n"; + $self->{type} = $self->type; + $self->set_consensus; + + if(!$self->{consensus} || length $self->{consensus} < $read_len * $min_percent_cons_of_read) { + $self->{error} = "No Consensus or consensus too short"; + print STDERR "FAILED\n"; + return; + } + if($self->{type} == INV && (!$self->{consensus2} || length $self->{consensus2} < $read_len * $min_percent_cons_of_read)) { + $self->{error} = "No Consensus or consensus too short"; + print STDERR "FAILED\n"; + return; + } + + foreach my $f (values(%default_filters) ){ + if(! $f->($self)){ + print STDERR "FAILED\n"; + return; + } + } + print STDERR "PASSED\n"; + return 1; +} + +sub _mapping_quality_filter { + my $self = shift; + if($self->{type} == INV) { + my $tmp1 = $self->_mapping_quality($self->{first_bp}, $self->{consensus}, 1); + my $tmp2 = $self->_mapping_quality($self->{second_bp}, $self->{consensus2}, 2); + return $tmp1 && $tmp2; + } + return $self->_mapping_quality($self->{first_bp}, $self->{consensus}, 1); +} + +sub _mapping_quality { + my $self = shift; + my $bp = shift; + my $con = shift; + my $which = shift; + $which = "" unless $which == 2; + my $l = length($con); + print STDERR "Mapping quality filter ... "; + + open my $QUERY, ">query.fa" or croak "can't open query.fa : $OS_ERROR"; + print $QUERY ">query\n", $con, "\n"; + close $QUERY; + + my ($range1, $range2) = $self->_find_bp_ref_range($bp, $con); + my ($s1, $e1, $t_s1, $t_e1) = $self->_map_consensus($range1); + return unless $e1; + my ($chr, $s, $e) = split /[:|-]/, $range2; + my $offset = 0; + if($s < $t_e1 && $e > $t_s1) { #overlap + my $tmp = substr $con, 0, $s; + if($s1 < $l - $e1) { + $offset = $e1; + $tmp = substr $con, $e1 + 1; + } + open my $QUERY, ">query.fa" or croak "can't open query.fa : $OS_ERROR"; + print $QUERY ">query\n", $tmp, "\n"; + close $QUERY; + } + my ($s2, $e2, $t_s2, $t_e2) = $self->_map_consensus($range2); + return unless $e2; + $s2 += $offset; + $e2 += $offset; + if($s1 < $s2) { + if($e1 < $s2) { + $self->{"insert" . $which} = substr($con, $e1, $s2 - $e1); + } + $self->{"c_start" . $which} = $s1; + $self->{"t_start" . $which} = $t_s1; + $self->{"start_chr" . $which} = $bp->{left_chr}; + $self->{"c_end" . $which} = $e2; + $self->{"t_end" . $which} = $t_e2; + $self->{"end_chr". $which} = $bp->{right_chr}; + $s1 = 1; + ($e1, $s2) = ($s2, $e1); + $e2 = $l; + } + else { + if($e2 < $s1) { + $self->{"insert" . $which} = substr($con, $e2, $s1 - $e2); + } + $self->{"c_start" . $which} = $s2; + $self->{"t_start" . $which} = $t_s2; + $self->{"start_chr" . $which} = $bp->{right_chr}; + $self->{"c_end" . $which} = $e1; + $self->{"t_end" . $which} = $t_e1; + $self->{"end_chr" . $which } = $bp->{left_chr}; + $s2 = 1; + ($e2, $s1) = ($s1, $e2); + $e1 = $l; + } + + my $n = 1; + foreach my $r( ($range1, $range2) ) { + my $r_a = $r eq $range1 ? $range2 : $range1; + my($chr, $s, $e) = split /[:|-]/, $r; + my($chr_a, $s_a, $e_a) = split /[:|-]/, $r_a; + my($s_c, $e_c) = $n == 1 ? ($s1, $e1) : ($s2, $e2); + $n++; + open my $QUERY, ">query.fa" or croak "can't open query.fa : $OS_ERROR"; + print $QUERY ">query\n", substr($self->{consensus}, $s_c, $e_c - $s_c + 1), "\n"; + close $QUERY; + my $hit = $aligner->run(-TARGET => $genome . ":" . $r, -QUERY => "query.fa"); + my $h = $hit->{query}; + return unless $h; + my $hsp = $h->next_hsp; + my @matches = $hsp->matches; + $hit = $mapper->run(-QUERY => "query.fa"); + foreach $h (@{$hit->{query}}) { + next if($h->{tchr} eq $chr && $h->{tstart} >= $s - $l && $h->{tend} <= $e + $l); # the same one + return if($h->{matches} >= $matches[0]); #find a better or as good as match + return if($h->{tchr} eq $chr_a && $h->{matches} - $matches[0] >= -5 && $self->{type} == CTX); + } + } + print STDERR "PASSED\n"; + return 1; +} + +sub _map_consensus { + my $self = shift; + my $range = shift; + my ($s_c, $e_c, $s_t, $e_t); + + my $hits = $aligner->run(-TARGET => $genome . ":" . $range, -QUERY => "query.fa"); + my ($chr, $s, $e) = split /[:|-]/, $range; + my $h = $hits->{query}; + return unless $h; # can't find a good enough match + my $hsp = $h->next_hsp; + my @matches = $hsp->matches; + ($s_c, $e_c, $s_t, $e_t) = ($hsp->start('query'), $hsp->end('query'), + $hsp->start('hit'), $hsp->end('hit')); + return if( ($e_c - $s_c + 1) * 0.97 > $matches[0]); # the alignment is not good + ($s_t, $e_t) = $hsp->strand == 1 ? ($s_t + $s, $e_t + $s) : ($e_t + $s, $s_t + $s); + return ($s_c, $e_c, $s_t, $e_t); +} + +sub _find_bp_ref_range { + my $self = shift; + my $bp = shift; + my $con = shift; + + my $l = int(length($con)/2); + my $ext = $l * 2; + $ext = 100000 if($RNASeq); + if($self->{type} == DEL && abs($bp->{left_pos} - $bp->{right_pos}) < $l ){ + $l = abs($bp->{left_pos} - $bp->{right_pos}) - 1; + } + + my($range1, $range2); + my ($chr, $p) = ($bp->{left_chr}, $bp->{left_pos}); + my $ort = $bp->{left_ort}; + my ($s, $e); + if($bp->{left_range} ) { + $s = $bp->{left_range}->[0] - $l; + $e = $bp->{left_range}->[1] + $l; + } + else { + $s = $e = $p; + $s = $ort eq "+" ? $p - $ext : $p - $l; + $e = $ort eq "+" ? $p + $l : $p + $ext; + } + $s = 1 if($s < 1); + $e = $bp->{left_pos} + $l if($self->{type} == DEL); + $range1 = $chr . ":" . $s . "-" . $e; + + ($chr, $p) = ($bp->{right_chr}, $bp->{right_pos}); + $ort = $bp->{left_ort}; + if($bp->{right_range} ) { + $s = $bp->{right_range}->[0] - $l; + $e = $bp->{right_range}->[1] + $l; + } + else { + $s = $e = $p; + $s = $ort eq "-" ? $p - $ext : $p - $l; + $e = $ort eq "-" ? $p + $l : $p + $ext; + } + $s = $bp->{right_pos} - $l if($self->{type} == DEL); + $s = 1 if($s < 1); + $range2 = $chr . ":" . $s . "-" . $e; + + return ($range1, $range2); +} + +sub _find_ref_range { + my $self = shift; + my $l = int(length($self->{consensus})/2); + my ($range1, $range2); + my ($bp1, $bp2) = ($self->{first_bp}, $self->{second_bp}); + + my ($chr, $p) = ($bp1->{chr}, $bp1->{pos}); + my ($s, $e); + my $ext = $l*2 ; + $ext = 100000 if($RNASeq); + if($self->{type} == DEL && abs($bp1->{pos} - $bp2->{pos}) < $l ){ + $l = abs($bp1->{pos} - $bp2->{pos}) - 1; + } + + my $ort = $bp1->{pos} == $bp1->{left_pos} ? + ($bp1->{left_ort} eq "+" ? "-" : "+") : $bp1->{right_ort}; + if($bp1->{range}) { + $s = $bp1->{range}->[0] - $l; + $s = 1 if($s < 1); + $e = $bp1->{range}->[1] + $l; + } + else { + $s = $e = $p; + $s = $ort eq "+" ? $p - $l : $p - $ext; + $e = $ort eq "+" ? $p + $ext : $p + $l; + } + $s = 1 if($s < 1); + $e = $bp1->{pos} + $l if($self->{type} == DEL); + + $range1 = $chr . ":" . $s . "-" . $e; + + ($chr, $p) = ($bp2->{chr}, $bp2->{pos}); + $ort = $bp2->{pos} == $bp2->{left_pos} ? + ($bp2->{left_ort} eq '+' ? '-' : '+') : $bp2->{right_ort}; + if($self->{type} == INV) { + $ort = $bp1->{pos} == $bp1->{left_pos} ? + $bp1->{right_ort} : ($bp1->{left_ort} eq '+' ? '-' : '+'); + } + + if($bp2->{range} && $self->{type} != INV) { + $s = $bp2->{range}->[0] - $l; + $e = $bp2->{range}->[1] + $l; + } + else { + $s = $e = $p; + $s = $ort eq "+" ? $p - $l : $p - $ext; + $e = $ort eq "+" ? $p + $ext : $p + $l; + } + $s = $bp2->{pos} - $l if($self->{type} == DEL); + $s = 1 if($s < 1); + + $range2 = $chr . ":" . $s . "-" . $e; + return($range1, $range2); +} + +my $GAP = -5; +my $MATCH = 2; +my $MISMATCH = -5; +sub _refine_bp { + my $self = shift; + my $h = shift; + my $ref_seq = shift; + my $con = $self->{consensus}; + my $hsp = $h->next_hsp; + my ($i, $j) = ($hsp->start("hit"), $hsp->start("query")); + if($hsp->strand == -1) { + $con = rev_comp($con); + $j = length($con) - $hsp->end("query"); + } + + my ($score, $g, $max_score, $s) = (0, 0, 0); + my ($current_s_r, $current_s_c) = ($i, $j); + my ($s_r, $s_c, $e_r, $e_c); + my $p; + my ($n_matches, $n_max_matches) = (0, 0); + my @bl_r = @{$hsp->gap_blocks('hit')}; + my @bl_c = @{$hsp->gap_blocks('query')}; + for( my $n = 0; $n < scalar @bl_r; $n++) { + my $m = $bl_r[$n]->[1]; + if($p) { #gaps + if (!$RNASeq || $bl_r[$n]->[0] - $p < 25) { + $score += $GAP * ($bl_r[$n]->[0] - $p); + $score = 0 if($score < 0 ) + } + $i += $m; + $j += $bl_c[$n]->[1]; + } + + while($m) { #matches / mismatches + ($current_s_r, $current_s_c, $n_matches) = ($i, $j, 0) + if($score == 0); #reset + + if(substr($ref_seq, $i, 1) eq substr($con, $j, 1)) { + $score += $MATCH; + $n_matches++; + } + else { + $score += $MISMATCH; + } + if($score > $max_score) { + $n_max_matches = $n_matches; + $max_score = $score; + ($s_r, $s_c) = ($current_s_r, $current_s_c); + ($e_r, $e_c) = ($i, $j); + } + $n_matches = $score = 0 if($score < 0 ); + $m--; $i++; $j++; + } + $p = $bl_r[$n]->[0] + $bl_r[$n]->[1]; + } + return($n_max_matches, $s_r, $e_r, $s_c, $e_c); +} + +sub _mapping_artifact_filter { + my $self = shift; + if($self->{type} == INV) { + return ($self->_mapping_artifact($self->{consensus}) + || $self->_mapping_artifact($self->{consensus2}) ); + } + return $self->_mapping_artifact($self->{consensus}); +} + +sub _mapping_artifact { + my $self = shift; + my $con = shift; + my $l = length $con; + print STDERR "Mapping artifact filter ... "; + $self->{error} = "Mapping Artifact"; + $l = $l - 5; + my $options = " minIdentity=97 -out=psl -nohead"; + $options = $options . " -maxIntron=5 " unless $RNASeq; + my $exp = $l; + if($self->{type} == DEL && $RNASeq) { + $exp = $self->{second_bp}{pos} - $self->{first_bp}{pos}; + $options = $options . " -maxIntron=" . ($exp + 100) if($exp > 75000); + } + + open my $QUERY, ">query.fa" or croak "can't open query.fa : $OS_ERROR"; + print $QUERY ">query\n", $con, "\n"; + close $QUERY; + my $tmp_mapping = $mapper->run(-QUERY => "query.fa", -OPTIONS => $options); + system("rm query.fa"); system("rm query.fa*"); + + $l = $l + 5; + if($RNASeq) { + #can't find a decent mapping for DEL event + return if($self->{type} == DEL && !$tmp_mapping->{query}); + foreach my $h(@{$tmp_mapping->{query}}) { + #next if($h->{tend} - $h->{tstart} < $exp - 10); + return if ($self->{type} != DEL && $h->{qend} - $h->{qstart} > $l - 10); + my $chr = $h->{tchr}; + $chr = "chr" . $chr unless $chr =~ m/^chr/; + my ($r_tree, $f_tree) = ($gm->sub_model($chr, "-"), $gm->sub_model($chr, "+")); + return unless ($f_tree && $r_tree); + my @f_genes = $f_tree->intersect([$h->{tstart}, $h->{tend}]); + my @r_genes = $r_tree->intersect([$h->{tstart}, $h->{tend}]); + return if(scalar @f_genes <= 1 && scalar @r_genes <= 1); + my @f_tmp_genes; + my @r_tmp_genes; + + foreach my $g (@f_genes) { + foreach my $block (@{$h->{blocks}}) { + if($g->val->overlap($block)){ + push @f_tmp_genes, $g; + last; + } + } + } + foreach my $g (@r_genes) { + foreach my $block (@{$h->{blocks}}) { + if($g->val->overlap($block)){ + push @r_tmp_genes, $g; + last; + } + } + } + return if(scalar @f_tmp_genes < 1 && scalar @r_tmp_genes < 1); + } + } + else { + foreach my $h (@{$tmp_mapping->{query}}){ + return if($h->{tend} - $h->{tstart} < $l + 10 && $h->{qend} - $h->{qstart} > $l - 5); + } + } + $self->{error} = undef; + return 1; +} + +sub _germline_sclip_filter { + my $self = shift; + print STDERR "Germline sclip filter\n"; + return 1 if(!$sam_g); + $self->{error} = "Germline Event"; + my $rtn = _compare_seq($self->{first_bp}) && + _compare_seq($self->{second_bp}); + $self->{error} = undef if($rtn); + return $rtn; + +} + +sub _compare_seq { + my $bp = shift; + my $seq; + if($bp->{pos} == $bp->{left_pos}) { + my $seg = $sam_d->segment($bp->{right_chr}, $bp->{right_pos} - + $germline_seq_width, $bp->{right_pos} + $germline_seq_width); + $seq = $seg->dna; + } + else { + my $seg = $sam_d->segment($bp->{left_chr}, $bp->{left_pos} - $germline_seq_width, + $bp->{left_pos} + $germline_seq_width); + $seq = $seg->dna; + } + my $sclip_seq = _get_sclip_seqs($sam_g, $bp->{chr}, $bp->{pos} - $germline_search_width, + $bp->{pos} + $germline_search_width); + + + open my $TARGET, ">target.fa" or croak "can't open target.fa : $OS_ERROR"; + print $TARGET ">target\n", $seq, "\n"; + close $TARGET; + + if(keys(%{$sclip_seq})) { + open my $QUERY, ">query.fa" or croak "can't open query.fa : $OS_ERROR"; + foreach my $s (keys(%{$sclip_seq})) { + print $QUERY ">", $s, "\n", $sclip_seq->{$s}, "\n"; + } + close($QUERY); + my $hits = $aligner->run(-TARGET => "target.fa", -QUERY => 'query.fa'); + return if(keys(%{$hits})); + } + return 1; +} + +sub _type_distance_filter { + my $self = shift; + my $type = $self->{type}; + $self->{error} = "Type Distance"; + print STDERR "Type distance filter\n"; + if($type == UNK) { #those events are not sure, always ignore them + print STDERR $self->to_full_string, "\n"; + return; + } + my ($bp1, $bp2) = ($self->{first_bp}, $self->{second_bp}); + return 1 if($bp2->{sc_reads} == 0); + my ($d1, $d2); + if($type == INS || $type == DEL) { + ($d1, $d2) = ($bp1->{left_pos} - $bp2->{left_pos}, + $bp1->{right_pos} - $bp2->{right_pos}); + } + if( $type == INV ){ + ($d1, $d2) = ($bp1->{left_pos} - $bp2->{left_pos}, + $bp2->{right_pos} - $bp1->{right_pos}); + } + if( $type == ITX ) { + ($d1, $d2) = ($bp1->{left_pos} - $bp2->{right_pos}, + $bp2->{left_pos} - $bp1->{right_pos}); + } + if($type == CTX ) { + if($bp1->{left_ort} eq $bp1->{right_ort} || + $bp1->{sc_reads} == 0 || $bp2->{sc_reads} == 0) { + ($d1, $d2) = ($bp1->{left_pos} - $bp2->{left_pos}, $bp1->{right_pos} - $bp2->{right_pos}) ; + } + else { + ($d1, $d2) = ($bp1->{left_pos} - $bp2->{right_pos}, $bp2->{left_pos} - $bp1->{right_pos}); + } + } + print STDERR $self->to_full_string, "\t", $d1, "\t", $d2, "\n"; + if($d1 <= 1 && $d2 <= 1 || abs($d1 - $d2) <= 5) { + $self->{error} = undef; + return 1; + } + if(abs($d1) > $max_bp_dist || abs($d2) > $max_bp_dist){ + print STDERR "Distance fail\n"; + return; + } + $self->{error} = undef; + return 1; # the distance is small enough +} + +sub _germline_indel_filter { + my $self = shift; + my $type = $self->{type}; + print STDERR "Germline INDEL FILTER test\n"; + return if(abs($self->{first_bp}{pos} - $self->{second_bp}{pos}) < 40); + return 1 if(!$sam_g); + return 1 if( $type != INS && $type != DEL); + return 1 if( abs($self->{first_bp}{pos} - $self->{second_bp}{pos}) > 50 ); + $self->{error} = "Germline INDEL"; + my ($start, $end ) = ($self->{first_bp}{left_pos}, $self->{first_bp}{right_pos}); + my $indel_len = abs($end - $start); + ($start, $end) = ($end, $start) if($start > $end); + + my $itr = $sam_g->features(-iterator => 1, -seq_id => $self->{first_bp}{chr}, + -start => $start, -end => $end); + while( my $a = $itr->next_seq ) { + next if($type == DEL && $a->cigar_str !~ m/D/); + next if($type == INS && $a->cigar_str !~ m/I/); + my $cigar_array = $a->cigar_array; + my $pos = $a->pos; + foreach my $ca (@{$cigar_array}) { + if($ca->[0] eq 'M') { + $pos += $ca->[1]; + next; + } + if($ca->[0] eq 'I' && $type == INS ){ + return if(abs($ca->[1] - $indel_len) <= 20 && abs($pos - $start) <= 20); + next; + } + if($ca->[0] eq 'D' && $type == DEL) { + return if(abs($ca->[1] - $indel_len) <= 20 && abs($pos - $start) <= 20); + $pos += $ca->[1]; + next; + } + } + } + $self->{error} = undef; + return 1; +} + +sub _tandem_repeat_filter { + my $self = shift; + my $con = $self->{consensus}; + print STDERR "low compexity filter\n"; + open my $QUERY, ">query.fa" or croak "can't open query.fa : $OS_ERROR"; + print $QUERY ">query\n", $self->{consensus}, "\n"; + close $QUERY; + + system("ptrfinder -seq query.fa -repsize $tr_min_size,$tr_max_size -minrep $tr_min_num > query.fa.rep"); + if(-e "query.fa.rep") { + open my $REP, "<query.fa.rep" or croak "can't open query.fa.rep:$OS_ERROR"; + while( my $line = <$REP>) { + chomp $line; + my($pattern, $len, $times, $s, $e) = + $line =~ m/PATTERN\s(\w+)\sLENGTH\s(\d+)\sTIMES\s(\d+)\sSTART\s(\d+)\sSTOP\s(\d+)\sID/; + if($self->{type} == DEL || $self->{type} == INS){ + if(abs($self->{first_bp}{left_pos} - $self->{first_bp}{right_pos}) < $tr_max_indel_size) { + print STDERR "Tandem repeat mediated INDEL!\n"; + close $REP; + return; + } + } + else{ + if(($s < 5 || length($self->{consensus}) - $e < 5) && $len * $times > 30 ) { + print STDERR "Tandem repeat mediated events!\n"; + close $REP; + return; + } + } + } + } + return 1; +} + +sub _get_sclip_seqs { + my ($sam, $chr, $start, $end) = @_; + my %rtn; + my $range = $chr . ":" . $start . "-" . $end; + + $sam->fetch($range, sub { + my $a = shift; + my $cigar_str = $a->cigar_str; + return if($cigar_str !~ /S/); + my ($sclip_len, $pos, $seq, $qual); + my @cigar_array = @{$a->cigar_array}; + if($cigar_array[0]->[0] eq 'S' ) { + $sclip_len = $cigar_array[0]->[1]; $pos = $a->start; + return if($pos < $start or $pos > $end); # the softclipping position is not in range + $seq = substr($a->query->dna, 0, $sclip_len ); + $rtn{$a->qname} = $seq if($sclip_len >= 10); + } + #if($cigar_str =~ m/S(\d+)$/) { + if($cigar_array[$#cigar_array]->[0] eq 'S') { + $sclip_len = $cigar_array[$#cigar_array]->[1]; + $pos = $a->end; + return if($pos < $start or $pos > $end); # the softclipping position is not in range + $seq = substr($a->qseq, $a->l_qseq - $sclip_len ); + $rtn{$a->qname} = $seq if($sclip_len >= 10); + } + } + ); + return \%rtn; +} + +sub _RNASeq_strand_filter { + my $self = shift; + my $type = $self->type; + print STDERR "RNASeq strand filter\n"; + return 1 unless $gm; + my ($bp1, $bp2) = ($self->{first_bp}, $self->{second_bp}); + my ($ort1, $ort2) = ($bp1->{left_ort}, $bp1->{right_ort}); + my ($chr1, $chr2) = ($bp1->{left_chr}, $bp1->{right_chr}); + my ($pos1, $pos2) = ($bp1->{pos}, $bp2->{pos}); + ($pos1, $pos2) = ($pos2, $pos1) if($bp1->{right_pos} == $pos1); + my ($tmp1, $tmp2) = ($chr1 =~ m/chr/ ? $chr1 : "chr" . $chr1, + $chr2 =~ m/chr/ ? $chr2 : "chr" . $chr2); + $tmp1 = "chrM" if($tmp1 eq "chrMT"); + $tmp2 = "chrM" if($tmp2 eq "chrMT"); + my ($r_tree1, $f_tree1) = ($gm->sub_model($tmp1, "-"), $gm->sub_model($tmp1, "+")); + my ($r_tree2, $f_tree2) = ($gm->sub_model($tmp2, "-"), $gm->sub_model($tmp2, "+")); + my ($g_ort1, $g_ort2); + return 1 unless($r_tree1 && $r_tree2 && $f_tree1 && $f_tree2); + $g_ort1 = "-" if(scalar($r_tree1->intersect([$pos1 - 5, $pos1])) > 0); + $g_ort1 = "+" if(scalar($f_tree1->intersect([$pos1 - 5, $pos1])) > 0); + $g_ort2 = "-" if(scalar($r_tree2->intersect([$pos2, $pos2 + 5])) > 0); + $g_ort2 = "+" if(scalar($f_tree2->intersect([$pos2, $pos2 + 5])) > 0); + return 1 unless($g_ort1 && $g_ort2); + return 1 if($g_ort1 eq $ort1 && $g_ort2 eq $ort2); + return 1 if($g_ort1 ne $ort1 && $g_ort2 ne $ort2); + return; +} + +# INS filter check to make sure that the INS part overlap with any gene +# if the INS part only in intron will return as a false positive +sub _RNASeq_INS_filter { + my $self = shift; + my $type = $self->{type}; + return 1 if($type != INS); + print STDERR "RNAseq INS filter\n"; + my ($bp1, $bp2) = ($self->{first_bp}, $self->{second_bp}); + my ($pos1, $pos2) = ($bp1->{pos}, $bp2->{pos}); + ($pos1, $pos2) = ($pos2, $pos1) if($pos2 < $pos1); + my $chr = $bp1->{chr}; + $chr = "chr" . $chr unless $chr =~ m/chr/; + my ($r_tree, $f_tree) = ($gm->sub_model($chr, "-"), $gm->sub_model($chr, "+")); + foreach my $tree( ($r_tree, $f_tree) ) { + my @tmp; + @tmp = $f_tree->intersect([$pos1, $pos2]); + foreach my $g (@tmp) { + return 1 if($g->val->overlap([$pos1, $pos2])); + } + } + return; +} + +sub _RNASeq_DEL_filter { + my $self = shift; + my $type = $self->type; + return 1 if($type != DEL); + print STDERR "RNAseq DEL filter\n"; + + my ($gene1, $gene2, $dist) = @{$self->get_genes}; + return if($gene1 eq $gene2); + my ($d_plus, $d_minus) = (0, 0); + $d_plus = $1 if($dist =~ m/(\d+)\(\+\)/); + $d_minus = $1 if($dist =~ m/(\d+)\(\-\)/); + return if($d_plus == 0 && $d_minus == 0); + my ($bp1, $bp2) = ($self->{first_bp}, $self->{second_bp}); + return if($bp1->{cover} == 0 || $bp2->{cover} == 0); + my($left_len, $right_len); + if(exists $bp1->{left_len}) { + ($left_len, $right_len ) = ($bp1->{left_len}, $bp1->{right_len}); + } + else { + $left_len = length($bp1->{sc_seq}) if($bp1->{sc_seq}); + $right_len = length($bp2->{sc_seq}) if($bp2->{sc_seq}); + } + return if($left_len < 30 || $right_len < 30); + return if($bp1->{sc_reads} < 10 || $bp2->{sc_reads} < 10); + return unless $bp1->{left_gene}->overlap([$bp1->{pos} - $left_len, $bp1->{pos}]); + return unless $bp1->{right_gene}->overlap([$bp2->{pos} + $right_len, $bp2->{pos}]); + return 1; +} + +sub get_statistics { + my $self = shift; + my $half_width = shift; + my $sam = $self->sam_d; + my @rtn; + + foreach my $bp ( ($self->{first_bp}, $self->{second_bp})) { + my ($chr, $pos ) = ($bp->{chr}, $bp->{pos}); + my $range = $chr . ":" . ($pos - $half_width) . "-" . ($pos + $half_width); + my ($n_seq, $n_rep, $total_pid) = (0, 0, 0); + + $sam->fetch($range, sub { + my $a = shift; + return unless ($a->start && $a->end); + return unless ($a->start >= $pos - $half_width && $a->end <= $pos + $half_width); + $n_seq++; + $total_pid += _cal_pid($a); + if($a->has_tag("XT")) { + $n_rep++ if($a->aux_get("XT") ne "U"); + } + } + ); + if($n_seq == 0) { + push @rtn, (0, 0); + } + else { + push @rtn, ($total_pid/$n_seq, $n_rep/$n_seq); + } + } + return @rtn; +} + +sub _cal_pid { + my $a = shift; + my ($ref, $matches, $query) = $a->padded_alignment; + my ($n_match, $n) = (0, 0); + for( my $i = 0; $i < length($matches); $i++) { + my $c = substr $matches, $i, 1; + $n_match++ if($c eq "|"); + $n++; + } + return 0 if($n == 0); + return $n_match/$n; +} + + +1; +__END__; + +=pod + +=head1 NAME +