Homec4science
Diffusion eSCT pipeline interoperability (master) Browse

R2915/example_pipelines/MelArray/container/break_point_mgroup_2017_6.pl
f5726a6975d6master

Last Change
Raw File

break_point_mgroup_2017_6.pl
View Options

#!/usr/bin/env perl
#
# 2017/6/14
# Yong Kong
#
use warnings;
use strict;
use Getopt::Long;
use File::Temp qw/ :POSIX /;
use Module::Load;
load Bio::Seq;
#push (@INC,"~/.cpan/build/");
my $bam;
my $chr;
my $start;
my $stop;
my $min_chr_num = 2;
my $max_dist = 1000; # group of starting pos
my $max_pair_dist = 1000; # pair mate dist
my $max_dist_to_bkp = 1000;
my $min_match = 8; # min of match in cigar
my $genome = "/home/yk336/genomes/human/b37/bwa-0.7.12/b37.fasta";
my $flank_seq_left = 200;
my $flank_seq_right = 200;
my $group_cutoff = 3;
my $test = 0;
GetOptions (
'bam=s' => \$bam,
'chr=s' => \$chr,
'start=i' => \$start,
'stop=i' => \$stop,
'min_chr_num=i' => \$min_chr_num,
'max_dist=i' => \$max_dist,
'max_pair_dist=i' => \$max_pair_dist,
'max_dist_to_bkp=i' => \$max_dist_to_bkp,
'min_match=i' => \$min_match,
'genome=s' => \$genome,
'flank_seq_left=i' => \$flank_seq_left,
'flank_seq_right=i' => \$flank_seq_right,
'group_cutoff=i' => \$group_cutoff,
'test' => \$test,
);
die "need bam chr start stop"
unless ($bam && $chr && $start && $stop);
my @modes = (
"A+ => B+",
"B- => A-",
"B+ => A+",
"A- => B-",
"A+ => B-",
"B+ => A-",
"A- => B+",
"B- => A+",
);
my %modes_cnts = ();
for my $m (@modes) {
$modes_cnts{$m} = 0;
}
my %modes_reads = ();
# A B
# AlBr =====>------ -----=====>
# <=====------ -----<=====
#
# ArBl -----======> ====>------
# -----<====== <====------
#
# AlBl =====>------ ====>------
# <=====------ <====------
#
# ArBr -----======> -----=====>
# -----<====== -----<=====
#
## inversion is Albl
my %mode_map = (
"AlBr" => "A+ => B+ OR B- => A-",
"ArBl" => "B+ => A+ OR A- => B-",
"AlBl" => "A+ => B- OR B+ => A-",
"ArBr" => "A- => B+ OR B- => A+",
);
my $tmp = tmpnam();
print "$tmp\n" if ($test);
# 0x400=1024 is PCR or optical duplicate;
my $cmd = "samtools view -F 0x400 $bam ${chr}:${start}-${stop} | grep SA: > $tmp";
print "$cmd\n" if ($test);
system($cmd);
&parse_bam($tmp);
=pod
From SAM specs:
SA Z Other canonical alignments in a chimeric alignment, formatted as a semicolon-delimited
list: (rname,pos,strand,CIGAR,mapQ,NM ;)+. Each element in the list represents a
part of the chimeric alignment. Conventionally, at a supplementary line, the first element
points to the primary line.
=cut
sub parse_bam {
my ($fn, ) = @_;
my $f = read_in_tab_delim_file($fn);
my @s1 = (); # seg1
my @s2 = (); # seg2
my @mate = (); ###### mate should go with seg2
my @ids = ();
my @seqs = ();
#
# chr1 chr2
# seq1 >>> seq2 >>>
# <<< mate
#
#
# chr1 chr2
# mate >>>
# <<< seq2 <<< seq1
#
my %reads = (); ### same chr, reads duplicates: sivT -chr 7 -start 140494267 -stop 140500161
my %uniquereads = (); ### picard does not label chimeric reads as pcr duplicates; see 5B/fusion_BRAF_7_140494267_140500161.txt
for my $e (@{$f}) {
### 2016/7/22 how about the pairs???
my $id = $e->[0];
#next if ($reads{$id});
#$reads{$id}++;
my $seq = $e->[9];
next if ($uniquereads{$seq});
$uniquereads{$seq}++;
my $flag = $e->[1];
my $rst = ($flag & 0x10) ? "-" : "+"; # read strand
my $pst = ($flag & 0x20) ? "-" : "+"; # mate strand
my $rchr = $e->[2];
my $rs = $e->[3];
my $rcigar = $e->[5];
my $pchr = $e->[6];
my $ps = $e->[7];
my $plen = $e->[8];
push(@s1, [$rchr, $rs, $rst, $rcigar]);
push(@mate, [$pchr, $ps, $pst]);
push(@ids, $id);
push(@seqs, $seq);
if ($test) {
print join("\t",
$flag,
$rchr,
$rs,
$rcigar,
$rst,
$pchr,
$ps,
$pst,
$plen,
$seq,
), "\n";
}
for (my $i=11; ; $i++) {
my $t = $e->[$i];
if ($t =~ /^SA:Z:(.*);/) { # use only the 1st of SA
my $sa = $1;
my @tmp = split(",", $sa);
print "\t$sa\n" if ($test);
my $chr = $tmp[0];
my $pos = $tmp[1];
my $st = $tmp[2];
my $cigar = $tmp[3];
push(@s2, [$chr, $pos, $st, $cigar]);
last;
}
}
}
### s1 and s2 should be the same size, as they are from the same sam line
if ($test) {
print "size of s1: ", scalar @s1, "\n";
print "size of s2: ", scalar @s2, "\n";
}
if (@s1) {
&analyze(\@s1, \@s2, \@mate, \@ids, \@seqs,);
}
}
sub analyze {
my ($s1, $s2, $m, $ids, $seqs, ) = @_;
my ($chr1, $chr1num, $gsidx1) = group($s1);
my ($chr2, $chr2num, $gsidx2) = group($s2);
if ($gsidx1 && $gsidx2) {
for my $gs1 (@{$gsidx1}) {
for my $gs2 (@{$gsidx2}) {
print "group $gs1 vs group $gs2...\n" if ($test);
for my $m (@modes) {
$modes_cnts{$m} = 0;
}
%modes_reads = ();
analyze_1_group($s1, $s2, $m, $ids, $seqs, $chr1, $chr2, $gs1, $gs2);
}
}
}
}
sub print_out {
my ($bpnum, $pair_support, $chr1, $a1, $v1, $chr2, $a2, $v2, ) = @_;
if ($a1 && $a1 > 0) {
my $mate_sup = &mates($chr2, $a1, $a2,);
print join("\t",
"chrA", "bkpA", "varA",
"chrB", "bkpB", "varB",
"#read_bkp",
"#read_bkp_mate",
"#read_pairs_over_bkp",
), "\n";
my ($pos1, $pos2) = (int($a1), int($a2));
print join("\t",
$chr1, $pos1, $v1, $chr2, $pos2, $v2,
$bpnum, $pair_support, $mate_sup,
), "\n";
print "\nNumber of reads (example: \"B+ => A-\" means chrB forward strand fuses to chrA negative strand):\n";
my %mymode = ();
for my $i (0..$#modes) {
#if ($modes_cnts{$modes[$i]} > 0) {
print join("\t", $modes[$i], $modes_cnts{$modes[$i]}, ), "\n";
#}
if ($modes[$i] eq "A+ => B+" || $modes[$i] eq "B- => A-") {
$mymode{"AlBr"} += $modes_cnts{$modes[$i]} if ($modes_cnts{$modes[$i]} > 0);
} elsif ($modes[$i] eq "B+ => A+" || $modes[$i] eq "A- => B-") {
$mymode{"ArBl"} += $modes_cnts{$modes[$i]} if ($modes_cnts{$modes[$i]} > 0);
} elsif ($modes[$i] eq "A+ => B-" || $modes[$i] eq "B+ => A-") {
$mymode{"AlBl"} += $modes_cnts{$modes[$i]} if ($modes_cnts{$modes[$i]} > 0);
} elsif ($modes[$i] eq "A- => B+" || $modes[$i] eq "B- => A+") {
$mymode{"ArBr"} += $modes_cnts{$modes[$i]} if ($modes_cnts{$modes[$i]} > 0);
}
}
my @cnt_sorted = sort{ $b->[1] <=> $a->[1]} map([$_, $mymode{$_}], keys %mymode);
print "\nThe fusion mode(s) are (see sequences below for details):\n";
for my $e (@cnt_sorted) {
print join("\t", @{$e}, $mode_map{$e->[0]}), "\n";
}
if (@cnt_sorted > 1) {
print "Since there are more than one fusion mode, the reads are not consistent.\n";
}
my ($seqA, $seqB, $seqa_det, $seqb_det);
my $mainmode = $cnt_sorted[0][0];
if ($mainmode eq "AlBr") {
$seqA = subseq($genome, $chr1, ($pos1 - $flank_seq_left) > 0 ? ($pos1 - $flank_seq_left) : 1, $pos1);
$seqB = subseq($genome, $chr2, $pos2, $pos2 + $flank_seq_right);
$seqa_det = join("_", $chr1, ($pos1 - $flank_seq_left) > 0 ? ($pos1 - $flank_seq_left) : 1, $pos1);
$seqb_det = join("_", $chr2, $pos2, $pos2 + $flank_seq_right);
} elsif ($mainmode eq "ArBl") {
$seqA = subseq($genome, $chr1, $pos1, $pos1 + $flank_seq_right);
$seqB = subseq($genome, $chr2, ($pos2 - $flank_seq_left) > 0 ? ($pos2 - $flank_seq_left) : 1, $pos2);
$seqa_det = join("_", $chr1, $pos1, $pos1 + $flank_seq_right);
$seqb_det = join("_", $chr2, ($pos2 - $flank_seq_left) > 0 ? ($pos2 - $flank_seq_left) : 1, $pos2);
} elsif ($mainmode eq "AlBl") {
$seqA = subseq($genome, $chr1, ($pos1 - $flank_seq_left) > 0 ? ($pos1 - $flank_seq_left) : 1, $pos1);
$seqB = subseq($genome, $chr2, ($pos2 - $flank_seq_left) > 0 ? ($pos2 - $flank_seq_left) : 1, $pos2);
$seqa_det = join("_", $chr1, ($pos1 - $flank_seq_left) > 0 ? ($pos1 - $flank_seq_left) : 1, $pos1);
$seqb_det = join("_", $chr2, ($pos2 - $flank_seq_left) > 0 ? ($pos2 - $flank_seq_left) : 1, $pos2);
} elsif ($mainmode eq "ArBr") {
$seqA = subseq($genome, $chr1, $pos1, $pos1 + $flank_seq_right);
$seqB = subseq($genome, $chr2, $pos2, $pos2 + $flank_seq_right);
$seqa_det = join("_", $chr1, $pos1, $pos1 + $flank_seq_right);
$seqb_det = join("_", $chr2, $pos2, $pos2 + $flank_seq_right);
}
my $rseqA = revComp($seqA);
my $rseqB = revComp($seqB);
if ($mainmode eq "AlBr") {
$seqA = $seqA . ">>>";
$rseqA = ">>>" . $rseqA;
$seqB = ">>>" . $seqB;
$rseqB = $rseqB . ">>>";
} elsif ($mainmode eq "ArBl") {
$seqA = ">>>" . $seqA ;
$rseqA = $rseqA . ">>>";
$seqB = $seqB . ">>>";
$rseqB = ">>>" . $rseqB;
} elsif ($mainmode eq "AlBl") {
$seqA = $seqA . ">>>";
$rseqA = ">>>" . $rseqA;
$seqB = $seqB . ">>>";
$rseqB = ">>>" . $rseqB;
} elsif ($mainmode eq "ArBr") {
$seqA = ">>>" . $seqA ;
$rseqA = $rseqA . ">>>";
$seqB = ">>>" . $seqB;
$rseqB = $rseqB . ">>>";
}
print "\nSeqA: $seqa_det (forward strand on top, reverse-complement on bottom; >>> indicates direction of fusion):\n";
print $seqA, "\n";
print $rseqA, "\n";
print "\nSeqB: $seqb_det (forward strand on top, reverse-complement on bottom; >>> indicates direction of fusion):\n";
print $seqB, "\n";
print $rseqB, "\n";
print "\nReads (fusion_mode, id, chr1, pos1, strand1, cigar1, chr2, pos2, strand2, cigar2, seq):\n";
for my $i (0..$#modes) {
if ($modes_cnts{$modes[$i]} > 0) {
for my $r (@{ $modes_reads{$modes[$i]} }) {
print join("\t", $modes[$i], @{$r},), "\n";
}
}
}
print "="x80;
print "\n\n";
}
}
### find mates in the regions that support this bkp
sub mates {
my ($c, $p1, $p2, ) = @_; # $c is chr2
my $mate_sup = 0;
my $tmp = tmpnam();
print "mates: $tmp\n" if ($test);
my $cmd = "samtools view -F1024 $bam ${chr}:${start}-${stop} | awk '(\$7 == \"$c\" || \$7 == \"=\")' | grep -v SA:Z> $tmp";
print "$cmd\n" if ($test);
system($cmd);
my $f = read_in_tab_delim_file($tmp);
for my $e (@{$f}) {
my $rs = $e->[3]; # read pos
my $mchr = $e->[6]; # mate chr
my $ps = $e->[7]; # mate pos
#print "\t$rs\t$ps\n" if ($test);
if ( (($mchr eq $c ) || ($c eq $chr && $mchr eq "=") ) &&
abs($ps - $p2) < $max_dist_to_bkp &&
abs($rs - $p1) < $max_dist_to_bkp
) {
$mate_sup++;
}
}
return $mate_sup;
}
### group based on cordinates
sub group {
my ($s, ) = @_; # chr, pos, st, cigar
my $total = @{$s};
my %h = ();
for my $e (@{$s}) {
my $chr = $e->[0];
$h{$chr}++;
}
# top chrs
my @ks = sort {$h{$b} <=> $h{$a}} keys %h;
my $topchr = $ks[0];
my $topnum = $h{$topchr};
print "$total\t$topchr\t$topnum\n" if ($test);
my %groups = ();
#my $maxgroup = -1;
my @groupindex = ();
# top chr must have $min_chr_num number of reads
if ($topnum >= $min_chr_num) {
my $gid = 0;
# sort index
my @itt = sort { $s->[$a][1] <=> $s->[$b][1] } 0..$#{$s};
my $prev_pos = $s->[$itt[0]][1];
#print "ITT: ", join(" ", @itt, ), "\n";
#print "FIRST POS: $prev_pos\n";
# the original order has not been changed
for my $i (0..$#{$s}) {
my $idx = $itt[$i];
my $chr = $s->[$idx][0];
my $pos = $s->[$idx][1];
if ($chr ne $topchr) {
push(@{$s->[$idx]}, -1); # group tag: not in the top chr
} else {
if (abs($pos - $prev_pos) > $max_dist) {
$prev_pos = $pos; # 2016/7/22: use the leftest pos of the group, not the pos of the previous member
$gid++;
}
push(@{$s->[$idx]}, $gid); # group tag
$groups{$gid}++;
}
if ($test) {
print join("\t",
@{$s->[$idx]},
), "\n";
}
}
#$maxgroup = (sort {$groups{$b} <=> $groups{$a}} keys %groups)[0];
#print "Max group: $maxgroup with $groups{$maxgroup} reads\n" if ($test);
for my $k (sort {$groups{$b} <=> $groups{$a}} keys %groups) {
if ($groups{$k} >= $group_cutoff) {
print "group $k...\n" if ($test);
push(@groupindex, $k);
}
}
if ($test) {
for my $k (sort {$groups{$b} <=> $groups{$a}} keys %groups) {
print "\t$k\t$groups{$k}\n";
}
}
}
###### hack for seepTm ~/SPORE/break_point.pl -bam recal_reads.bam -chr 3 -start 12641914 -stop 12645634 -genome /home/yk336/genomes/human/b37/bwa-0.7.12/b37.fasta -test
### $maxgroup = 3 if ($maxgroup == 4);
return ($topchr, $topnum, \@groupindex);
}
# filter out small indels in cigar
sub check_cigar {
my ($cigar, ) = @_;
#print ">$cigar<\n";
my @nums = split(/\D/, $cigar);
my @lets = split(/\d+/, $cigar);
shift(@lets);
my @nums2 = ();
my @lets2 = ();
for my $i (0..$#nums) {
#print ">$lets[$i]<\n";
if ($nums[$i] > $min_match) {
push(@nums2, $nums[$i]);
push(@lets2, $lets[$i]);
}
}
return(\@nums2, \@lets2);
}
sub bkp {
my ($n, $l, $p, $s, ) = @_;
my ($bk, $segment,);
#print "$l->[0]\t$l->[1]\n" if ($test);
if ($l->[0] eq "M" && $l->[1] =~ /[SH]/) {
$bk = $p + $n->[0];
$segment = ($s eq "+") ? 1 : 2;
} elsif ($l->[0] =~ /[SH]/ && $l->[1] eq "M") {
$bk = $p;
$segment = ($s eq "+") ? 2 : 1;
}
return ($bk, $segment);
}
sub avg_var {
my ($l, ) = @_;
my ($s0, $s1, $s2) = (0, 0, 0);
for my $e (@{$l}) {
$s0++;
$s1 += $e;
$s2 += $e*$e;
}
my $avg = $s0 > 0 ? $s1/$s0 : undef;
my $var = $s0 > 0 ? $s2/$s0 - $avg*$avg : undef;
return ($avg, $var);
}
sub analyze_1_group {
my ($s1, $s2, $m, $ids, $seqs, $chr1, $chr2, $gs1, $gs2, ) = @_;
my $bpnum = 0; # number of reads to support bkp
my $pair_support = 0;
my ($a1, $v1, ) = (-1, -1);
my ($a2, $v2, ) = (-1, -1);
my (@bks1, @bks2) = ((), ());
if ($test) {
print "s1...\n";
for my $i (0..$#{$s1}) {
print "$i...\t";
print join("\t", @{$s1->[$i]}), "\n";
}
print "s2...\n";
for my $i (0..$#{$s2}) {
print "$i...\t";
print join("\t", @{$s2->[$i]}), "\n";
}
}
for my $i (0..$#{$s1}) {
if ($s1->[$i][-1] == $gs1 && $s2->[$i][-1] == $gs2) {
my $s1_chr = $s1->[$i][0];
my $s1_pos = $s1->[$i][1];
my $s1_st = $s1->[$i][2];
my $s1_cigar = $s1->[$i][3];
my $s2_chr = $s2->[$i][0];
my $s2_pos = $s2->[$i][1];
my $s2_st = $s2->[$i][2];
my $s2_cigar = $s2->[$i][3];
my $m_chr = $m->[$i][0];
my $m_pos = $m->[$i][1];
my $m_st = $m->[$i][2];
$m_chr = $s1_chr if ($m_chr eq "=");
my ($s1_num, $s1_let) = check_cigar($s1_cigar);
my ($s2_num, $s2_let) = check_cigar($s2_cigar);
if ($test) {
print join("\t",
"s1",
@{$s1->[$i]},
@{$s1_num},
@{$s1_let},
$ids->[$i],
$seqs->[$i],
), "\n";
print join("\t",
"s2",
@{$s2->[$i]},
@{$s2_num},
@{$s2_let},
), "\n";
print join("\t",
"m",
@{$m->[$i]},
), "\n";
}
if (@{$s1_num} < 2 || @{$s1_let} < 2 ||
@{$s2_num} < 2 || @{$s2_let} < 2) {
next;
} else {
# now we have 47M29S or 46H30M
my ($bk1, $seg1) = bkp($s1_num, $s1_let, $s1_pos, $s1_st);
my ($bk2, $seg2) = bkp($s2_num, $s2_let, $s2_pos, $s2_st);
if ($seg1 == $seg2) {
warn "$seg1 == $seg2";
exit 0;
}
### check if the mate goes with seg2
if ( ($seg1 == 2 && $m_chr eq $s1_chr && $m_st ne $s1_st &&
abs($m_pos - $s1_pos) < $max_pair_dist)
||
($seg2 == 2 && $m_chr eq $s2_chr && $m_st ne $s2_st &&
abs($m_pos - $s2_pos) < $max_pair_dist) ) {
$pair_support++;
}
### 2016/7/22: why swap? will exchange of s1 and s2 cause problem?
### swap seg1 and seg2
if ($s1_chr eq $s2_chr) {
if ($bk1 > $bk2) {
($bk1, $bk2) = ($bk2, $bk1);
($seg1, $seg2) = ($seg2, $seg1);
($s1, $s2) = ($s2, $s1);
}
}
### 2016/7/22: Is this a bug???
#push(@{$s1}, $bk1);
#push(@{$s2}, $bk2);
push(@bks1, $bk1);
push(@bks2, $bk2);
$bpnum++;
print "BKP: $bk1 ($seg1)\t$bk2 ($seg2)\n" if ($test);
##### all possible 8 modes
##### note different meanings of notations
##### "2-1+" means seg1 is the 2nd fragment mapped to the negative strand,
##### while seq2 is the 1st fragment mapped to the forward strand
##### "B+ => A-" means chrB forward strand fuses to chrA negative strand
##### note a -> b
##### also the complement a and b exchange
my $mode;
if ($seg1 == 1 && $s1_st eq "+" && $seg2 == 2 && $s2_st eq "+") {
$mode = "A+ => B+";
# 1+2+ AlBr
# --->---a b--->---
# ------- -------
} elsif ($seg1 == 2 && $s1_st eq "-" && $seg2 == 1 && $s2_st eq "-") {
$mode = "B- => A-";
# 2-1- AlBr
# ------ -------
# ---<---b a---<---
} elsif ($seg1 == 2 && $s1_st eq "+" && $seg2 == 1 && $s2_st eq "+") {
$mode = "B+ => A+";
# 2+1+ ArBl
# b--->--- --->---a
# ------- -------
} elsif ($seg1 == 1 && $s1_st eq "-" && $seg2 == 2 && $s2_st eq "-") {
$mode = "A- => B-";
# 1-2- ArBl
# ------- -------
# a---<--- ---<---b
} elsif ($seg1 == 1 && $s1_st eq "+" && $seg2 == 2 && $s2_st eq "-") {
$mode = "A+ => B-";
# 1+2- AlBl
# --->---a -------
# ------- ---<---b
} elsif ($seg1 == 2 && $s1_st eq "-" && $seg2 == 1 && $s2_st eq "+") {
$mode = "B+ => A-";
# 2-1+ AlBl
# ------- --->---a
# ---<---b -------
} elsif ($seg1 == 1 && $s1_st eq "-" && $seg2 == 2 && $s2_st eq "+") {
$mode = "A- => B+";
# 1-2+ ArBr
# ------- b--->---
# a---<--- -------
} elsif ($seg1 == 2 && $s1_st eq "+" && $seg2 == 1 && $s2_st eq "-") {
$mode = "B- => A+";
# 2+1- ArBr
# b--->--- -------
# ------- a---<---
}
$modes_cnts{$mode}++;
push(@{$modes_reads{$mode}}, [$ids->[$i], @{$s1->[$i]}[0..3], @{$s2->[$i]}[0..3], $seqs->[$i]]);
}
}
}
print "#read: $bpnum\n" if ($test);
($a1, $v1, ) = avg_var(\@bks1);
($a2, $v2, ) = avg_var(\@bks2);
print join("\t", $a1, $v1,), "\n" if ($test);
print join("\t", $a2, $v2,), "\n" if ($test);
print_out($bpnum, $pair_support, $chr1, $a1, $v1, $chr2, $a2, $v2, );
}
sub read_in_tab_delim_file {
my ($fn, ) = @_;
my @r = ();
my $fh = &get_in_fh($fn);
while (my $li = <$fh>) {
chomp($li);
next if ($li =~ /^\#/);
my @t = split(/\t/, $li);
push(@r, [ @t ]);
}
close($fh);
return \@r;
}
# subseq from indexed fasta file
sub subseq {
my ($db, $chr, $start, $end, ) = @_;
$chr = "\"$chr\"";
my $prog = "samtools faidx";
my $cmd = "$prog $db $chr:${start}-${end}";
my $out = `$cmd`;
my @t = split("\n", $out);
return join('', @t[1..$#t])
}
sub revComp {
my ($seq, ) = @_;
my $seq_o = Bio::Seq->new(-seq => $seq, -alphabet => 'dna');
return $seq_o->revcom->seq();
}
sub get_in_fh {
my ($fn, ) = @_;
my $fh;
open($fh, "<$fn") || die "cannot open file $fn!";
return $fh;
}
c4science · Help