perl/bsim/testLVR.pl

   1 #!/usr/bin/env perl
   2 use strict;
   3 use warnings;
   4 use Data::Dump qw(ddx);
   5
   6 use FindBin 1.51 qw($RealBin);
   7 use lib "$RealBin";
   8 use simVirusInserts;
   9
  10 sub InsertParts2RealPos($$$$$) {
  11         my ($r1fr,$PartsRef,$MappedChr,$simLMR,$VirSLR) = @_;
  12         my @Poses;
  13         for my $Part (@$PartsRef) {
  14                 $Part =~ /^(\d+)(\w+)$/ or die "[$Part]";
  15                 my ($slen,$type)=($1,$2);
  16                 my $tmp;
  17                 if ($type eq 'V') {
  18                         if ( ($VirSLR->[0] eq '+' and $r1fr eq 'f') or ($VirSLR->[0] eq '-' and $r1fr eq 'r') ) {
  19                                 $tmp = $VirSLR->[1]; # LeftPos => not to: + $slen;
  20                         } else {
  21                                 $tmp = 1+ $VirSLR->[2] - $slen;
  22                         }
  23                 } else {
  24                         $tmp = 1+ $simLMR->[1] - $slen if $type eq 'L'; # `1+ $simLMR->[1] - $slen` can also be `1+ $simLMR->[0] + $innerPos` for 'L' && Read1。
  25                         $tmp = 1+ $simLMR->[1] if $type eq 'R'; # LeftPos => not to add $slen
  26                 }
  27                 push @Poses,$tmp.$type."+$slen";
  28         }
  29         return @Poses;
  30 }
  31 sub getRealPos($$$$$$$$$$) {
  32         my ($r1fr,$innerPos, $InsertSize,$ReadLen,$VirSLR, $MappedChr,$rSMS,$r12,$strand,$simLMR)=@_;
  33         my @ret;        # 暂且忽略 $strand
  34         my $VirFrag = $VirSLR->[2] - $VirSLR->[1] +1;
  35         my ($FiveT,$ThreeT) = getInsertPos($r1fr,$innerPos,$InsertSize,$ReadLen,$r12);
  36         my @Parts = InsertPos2InsertParts($InsertSize,$ReadLen,$VirFrag, $FiveT,$ThreeT);
  37         my @Poses = InsertParts2RealPos($r1fr,\@Parts,$MappedChr,$simLMR,$VirSLR);
  38         ddx [$FiveT,$ThreeT,@Parts];
  39         return @Poses;
  40 }
  41
  42 my $str = 'sf176_Ref_43898202_43898351_43898501_Vir_+_1878_1922_R_150_90';
  43 $str =~ /^s([fr])(\d+)_Ref_(\d+)_(\d+)_(\d+)_Vir_([+-])_(\d+)_(\d+)_R_(\d+)_(\d+)$/ or die "$str";      # _([\d\|\-LVR]+)
  44 my ($r1fr,$innerPos,$RefLeft,$RefMiddle,$RefRight,$VirStrand,$VirLeft,$VirRight,$InsertSize,$ReadLen) = ($1,$2,$3,$4,$5,$6,$7,$8,$9,$10);
  45
  46 my ($r12R1,$r12R2,$r1SMS,$r2SMS)=(1,2);
  47 my ($strandR1,$strandR2)=('+','+');
  48 my ($refR1,$refR2)=('NA','NA');
  49
  50 my @Poses1 = getRealPos($r1fr,$innerPos,$InsertSize,$ReadLen, [$VirStrand,$VirLeft,$VirRight], $refR1,$r1SMS,$r12R1,$strandR1, [$RefLeft,$RefMiddle,$RefRight]);
  51 my @Poses2 = getRealPos($r1fr,$innerPos,$InsertSize,$ReadLen, [$VirStrand,$VirLeft,$VirRight], $refR2,$r2SMS,$r12R2,$strandR2, [$RefLeft,$RefMiddle,$RefRight]);