diff scripts/pickUniqPos_MEM.pl @ 0:28d1a6f8143f draft

planemo upload for repository https://github.com/portiahollyoak/Tools commit 132bb96bba8e7aed66a102ed93b7744f36d10d37-dirty

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/pickUniqPos_MEM.pl	Mon Apr 25 13:08:56 2016 -0400
@@ -0,0 +1,56 @@
+#!/share/bin/perl
+use Bio::Seq;
+use List::Util qw(sum);
+
+die "perl $0 <sam>\n" if @ARGV<1;
+open in,$ARGV[0];
+my %pe;
+while(<in>)
+{
+	chomp;
+	my @f=split/\t/,$_,12;
+	## read number 1 or 2
+	my ($rnum)=$f[1]=~/(\d)$/;
+
+	## XT:A:* 
+        my $xt="";
+        my @a=split(/\s+/, $_);
+        my $as=0;
+        my $xs=0;
+        for my $i (11..$#a) {
+            if ($a[$i] =~ /^AS:i:/) {
+                $a[$i] =~ s/AS:i://;
+                $as=$a[$i];
+            }
+            elsif ($a[$i] =~ /^XS:i:/) {
+                $a[$i] =~ s/XS:i://;
+                $xs=$a[$i];
+            }
+            if (($xs > 0) && ($as-$xs <= $ARGV[1])) {$xt="R";}
+            else {$xt="U";}
+        }
+
+	my $strand="+";
+	## revcomp
+	if($f[1]=~/r/)
+        {
+                my $seq=Bio::Seq->new(-seq=>$f[9]);
+                $f[9]=$seq->revcom->seq;
+		$strand="-";
+        }
+
+	## parse CIGAR
+	if($xt eq "U")
+        {
+                # CIGAR
+                my (@cigar_m)=$f[5]=~/(\d+)M/g;
+                my (@cigar_d)=$f[5]=~/(\d+)D/g;
+                my (@cigar_s)=$f[5]=~/(\d+)S/g;
+                my (@cigar_i)=$f[5]=~/(\d+)I/g;
+                my $aln_ln=sum(@cigar_m,@cigar_d);
+		
+		print $f[2],"\t",$f[3]-1,"\t",$f[3]-1+$aln_ln,"\t$f[0]/$rnum\t",$f[9],"\t",$strand,"\n";
+	}
+}
+close in;
+