Mercurial > repos > jjohnson > crest
view SVUtil.pm @ 1:4f6952e0af48 default tip
CREST - add crest.loc.sample
author | Jim Johnson <jj@umn.edu> |
---|---|
date | Wed, 08 Feb 2012 16:08:01 -0600 |
parents | acc8d8bfeb9a |
children |
line wrap: on
line source
package SVUtil; use strict; use Carp; use English; use Bio::SeqIO; use Bio::DB::Sam; use File::Spec; use File::Path qw(remove_tree); use File::Temp qw/ tempfile tempdir /; use Cwd; use SCValidator qw(LEFT_CLIP RIGHT_CLIP); require Exporter; our @ISA = ('Exporter'); our @EXPORT = qw(is_new_pair clip_fa_file prepare_reads parse_range get_input_bam get_work_dir is_PCR_dup read_fa_file find_smallest_cover get_sclip_reads rev_comp); our @EXPORT_OK = qw($use_scratch $spacer_seq $spcer_seq_len $work_dir); our $use_scratch = 1; our $work_dir; sub rev_comp { my $str = shift; $str = reverse $str; $str =~ tr/ACTG/TGAC/; return $str; } # something about binormial distribution sub choose { my ($n,$k) = @_; my ($result,$j) = (1,1); return 0 if $k>$n||$k<0; $k = ($n - $k) if ($n - $k) <$k; while ($j <= $k ) { $result *= $n--; $result /= $j++; } return $result; } sub dbinorm { #desity funtion my($n, $k, $p) = @_; return unless($n > 0 && $k >= 0 && $n >= $k && $p >= 0 && $p <=1); my $prob = log(choose($n, $k)) + ($n-$k)*log(1-$p) + $k*log($p); return exp($prob); } sub find_smallest_cover { my ($c, $p, $crit) = @_; my @s_cover; for( my $i = 1; $i < $c; $i++) { my $n = 1; while(1) { my $tmp = 0; for( my $j = 0; $j <= $i && $j < $n; $j++) { $tmp += dbinorm($n, $j, $p); } if( $tmp < $crit) { $s_cover[$i] = $n - 1; last; } $n++; } } return \@s_cover; } sub read_fa_file { my $file = shift; my $in = Bio::SeqIO->new(-file => $file, -format => 'Fasta'); my %seqs; while( my $seq=$in->next_seq()) { $seqs{$seq->display_id} = $seq->seq; } return \%seqs; } # rmdup: remove PCR dumplication related code # you can replace it using your own criteria sub is_new_pair { my ($a, $sclip_len, $pairs) = @_; foreach my $p (@{$pairs}) { return undef if($a->start == $p->[0] && $a->end == $p->[1] && $a->mate_start == $p->[2] && $a->mate_end == $p->[3] && $sclip_len == $p->[4]); } return 1; } # #clipping sequnce related code # # the spacer sequence, you can replace it with your own one my $spacer_seq = 'ACGTACGT' . 'CGTA' x 6 . 'ATGCATGC'; my $spacer_seq_len = length($spacer_seq); # clip the reads that with spacer_seq at one end sub clip_fa_file { my ($file, $ort) = @_; my $out = $file . ".clip.fa"; my $in = Bio::SeqIO->new( -file => $file, -format => 'Fasta'); open my $OUT, ">$out" or croak "can't open $file.fa:$OS_ERROR"; while( my $seq=$in->next_seq()) { print $OUT ">", $seq->display_id, "\n"; if($ort eq RIGHT_CLIP) { print $OUT $seq->subseq($spacer_seq_len + 1, $seq->length), "\n"; } else { print $OUT $seq->subseq(1, $seq->length - $spacer_seq_len), "\n"; } } close($OUT); return $out; } # prepare the reads into a fasta file # due to most assembler can't handle short read well, we add # the spacer_seq to the start/end of each sequence arrording to # the softclip postion sub prepare_reads { my($r_reads, $clip) = @_; my ($fh, $fa_name) = tempfile( DIR => $work_dir, SUFFIX => '.fa'); my $qual_name = $fa_name . ".qual"; my $spacer_qual = ' 50 ' x $spacer_seq_len; open my $qual_fh, ">$qual_name" or croak "can't open $qual_name:$OS_ERROR"; foreach my $r (@{$r_reads}) { print $fh ">", $r->{name}, "\n"; my $tmp = $r->{seq}; $tmp = defined $clip ? ($clip == RIGHT_CLIP ? $spacer_seq . $tmp : $tmp . $spacer_seq) : $tmp; print $fh $tmp, "\n"; next unless $r->{qual}; print $qual_fh ">", $r->{name}, "\n"; $tmp = join(" ", @{$r->{qual}}); $tmp = defined $clip ? ($clip == RIGHT_CLIP ? $spacer_qual . $tmp : $tmp . $spacer_qual) : $tmp; print $qual_fh $tmp, "\n"; } close $qual_fh; return $fa_name; } # check PCR duplicate sub is_PCR_dup { my ($a, $pairs, $sclip_len) = @_; my ($s, $e, $ms, $me ) = ($a->start, $a->end, $a->mate_start ? $a->mate_start : 0, $a->mate_end ? $a->mate_end : 0); foreach my $p (@{$pairs}) { return 1 if($s == $p->[0] && $e == $p->[1] && $ms == $p->[2] && $me == $p->[3] && $sclip_len == $p->[4]); } return; } # parse the range of input, format is: chr:start-end # start and end is optional sub parse_range { my $range = shift; my ($chr, $start, $end); my @field = split /:/, $range; $chr = $field[0]; # $chr = substr($chr, 3) if($chr =~ /^chr/); if(@field > 1) { @field = split /-/, $field[1]; $start = $field[0]; $end = $field[1] if($field[1]); if($start !~ m/^\d+$/ or $end !~ m/^\d+$/) { croak "wrong range format, need to be: chr:start-end"; } } return ($chr, $start, $end); } # St. Jude only # grab the bam file from the bam dir my $raw_bam_dir = "/nfs_exports/genomes/1/PCGP/BucketRaw"; sub get_input_bam { my $sample = shift; my ($group) = $sample =~ /(.*?)\d+/; my $tmp_bam_dir = File::Spec->catdir($raw_bam_dir, $group); opendir(my $dh, $tmp_bam_dir) or croak "can't open directory: $raw_bam_dir: $OS_ERROR"; my @files = grep { /^$sample-.*bam$/ } readdir($dh); close $dh; return File::Spec->catfile($raw_bam_dir, $group, $files[0]); } #scratch related code sub get_work_dir { my %options = @_; my $scratch_dir; if($options{-SCRATCH}){ $scratch_dir = (exists $options{-DIR} ? $options{-DIR} : (-e '/scratch_space' ? '/scratch_space' : undef) ); } my $dir; if($scratch_dir) { $scratch_dir = File::Spec->rel2abs($scratch_dir); if($ENV{'JOB_ID'}) { $dir = File::Spec->catdir($scratch_dir, $ENV{'JOB_ID'}) } else { $dir = tempdir(DIR => $scratch_dir); } } else { $dir = tempdir(); } mkdir($dir) if( ! -e $dir); $work_dir = $dir; return $dir; } sub get_sclip_reads { my %args = @_; my ($sam, $chr, $start, $end, $clip, $min_len, $validator, $paired, $rmdup, $full_seq, $extra_base) = ($args{-SAM}, $args{-CHR}, $args{-START}, $args{-END}, $args{-CLIP}, $args{-MINCLIPLEN}, $args{-VALIDATOR}, $args{-PAIRED}, $args{-RMDUP}, $args{-FULLSEQ}, $args{-EXTRA}); $extra_base = 0 unless $extra_base; my @rtn; $min_len = 0 unless $min_len; my $max_sclip_len = 0; my $range = $chr; $range = $chr . ":" . $start . "-" . $end if($start && $end); my @pairs; my $strand_validator = $paired ? 1 : 0; $sam->fetch($range, sub { my $a = shift; my $cigar_str = $a->cigar_str; my @cigar_array = @{$a->cigar_array}; return if($cigar_str !~ m/S/); return unless($a->start && $a->end); #return if(!$a->proper_pair && $paired); if($paired && !$strand_validator) { $validator->add_validator("strand_validator"); } if($paired && !$a->proper_pair) { $validator->remove_validator("strand_validator"); $strand_validator = 0; } my ($sclip_len, $pos, $seq, $qual); my @tmp = $a->qscore; if($cigar_array[0]->[0] eq 'S' && $clip == LEFT_CLIP ) { $pos = $a->start; return if($pos < $start or $pos > $end); # the softclipping position is not in range $sclip_len = $cigar_array[0]->[1]; #print $cigar_str, "\t", scalar @pairs, "\n"; return unless($validator && $validator->validate($a, LEFT_CLIP)); $max_sclip_len = $sclip_len if($max_sclip_len < $sclip_len); return if(@pairs > 0 && $rmdup && is_PCR_dup($a, \@pairs, $sclip_len)); #it's a PCR dumplicate $seq = $full_seq ? $a->query->dna : substr($a->query->dna, 0, $sclip_len + $extra_base); @tmp = $full_seq ? @tmp : @tmp[0 .. ($sclip_len -1 + $extra_base)]; push @pairs, [$a->start, $a->end, $a->mate_start ? $a->mate_start : 0, $a->mate_end ? $a->mate_end : 0, $sclip_len] if($rmdup); my $qscore = $a->qscore; push @rtn, { name => $a->qname, seq => $seq, qual => \@tmp, full_seq => $a->query->dna, full_qual => $qscore, pos => $pos, }; } if($cigar_array[$#cigar_array]->[0] eq 'S' && $clip == RIGHT_CLIP ) { $pos = $a->end; return if($pos < $start or $pos > $end); # the softclipping position is not in range # print $cigar_str, "\t", scalar @pairs, "\n"; return unless $validator->validate($a, RIGHT_CLIP); $sclip_len = $cigar_array[$#cigar_array]->[1]; $max_sclip_len = $sclip_len if($max_sclip_len < $sclip_len); return if(@pairs > 0 && $rmdup && is_PCR_dup($a, \@pairs, $sclip_len)); #it's a PCR dumplicate $seq = $full_seq ? $a->qseq : substr($a->qseq, $a->l_qseq - $sclip_len - $extra_base ); @tmp = $full_seq ? @tmp : @tmp[($a->l_qseq - $sclip_len - $extra_base) .. ($a->l_qseq - 1)]; push @pairs, [$a->start, $a->end, $a->mate_start ? $a->mate_start : 0, $a->mate_end ? $a->mate_end : 0, $sclip_len] if($rmdup); my $qscore = $a->qscore; push @rtn, { name => $a->qname, seq => $seq, qual => \@tmp, full_seq => $a->query->dna, full_qual => $qscore, pos => $pos, }; } } ); if($max_sclip_len >= $min_len) { return @rtn; } else{ undef @rtn; return @rtn; } } my $start_dir; BEGIN { $start_dir = getcwd; } END { chdir($start_dir); remove_tree($work_dir) if($use_scratch && $work_dir && $start_dir ne $work_dir); } 1;