Mercurial > repos > portiahollyoak > temp

comparison scripts/pickSoftClipping.over.pl @ 21:9672fe07a232 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/portiahollyoak/Tools commit 0fea84d05f8976b8360a8b4943ecb01b87e3ade0-dirty
author mvdbeek
date Mon, 05 Dec 2016 09:58:47 -0500
parents 28d1a6f8143f
children
comparison
equal deleted inserted replaced
20:6e02b9179a24 21:9672fe07a232
36 my $upper=$a[4]+100; 36 my $upper=$a[4]+100;
37 my $chr_num=$a[2]; 37 my $chr_num=$a[2];
38 $chr_num =~ s/chr//; 38 $chr_num =~ s/chr//;
39 if (($chrs{$a[2]} == 1) && (! defined $chrs{$chr_num})) {$chr_num=$a[2];} 39 if (($chrs{$a[2]} == 1) && (! defined $chrs{$chr_num})) {$chr_num=$a[2];}
40 system("samtools view -bu $title $chr_num\:$lower\-$upper > temp.bam"); 40 system("samtools view -bu $title $chr_num\:$lower\-$upper > temp.bam");
41 system("samtools view -Xf 0x2 temp.bam > temp.sam"); 41 system("samtools view -f 0x2 temp.bam > temp.sam");
42 42
43 my $leftseq=""; 43 my $leftseq="";
44 my $rightseq=""; 44 my $rightseq="";
45 45
46 my $ll=$a[3]-150; 46 my $ll=$a[3]-150;
77 while(<in>) 77 while(<in>)
78 { 78 {
79 chomp; 79 chomp;
80 my @f=split/\t/,$_,12; 80 my @f=split/\t/,$_,12;
81 ## read number 1 or 2 81 ## read number 1 or 2
82 my ($rnum)=$f[1]=~/(\d)$/; 82 #my ($rnum)=$f[1]=~/(\d)$/;
83 my $rnum=1;
84 if (($f[1] & 128) == 128) {$rnum=2;}
83 85
84 ## XT:A:* 86 ## XT:A:*
85 my ($xt)=$f[11]=~/XT:A:(.)/; 87 my ($xt)=$f[11]=~/XT:A:(.)/;
86 88
87 my $CIGAR=$f[5]; 89 my $CIGAR=$f[5];
92 my $coor=-10; 94 my $coor=-10;
93 my $overcoor=-10; 95 my $overcoor=-10;
94 my $strand=""; 96 my $strand="";
95 my @z=split(/M/, $f[5]); 97 my @z=split(/M/, $f[5]);
96 98
97 if (($f[5]=~/S$/)&&($f[1]=~/r/)) 99 if (($f[5]=~/S$/)&&(($f[1] & 16) == 16))
98 { 100 {
99 my (@cigar_m)=$f[5]=~/(\d+)M/g; 101 my (@cigar_m)=$f[5]=~/(\d+)M/g;
100 my (@cigar_d)=$f[5]=~/(\d+)D/g; 102 my (@cigar_d)=$f[5]=~/(\d+)D/g;
101 my (@cigar_s)=$f[5]=~/(\d+)S/g; 103 my (@cigar_s)=$f[5]=~/(\d+)S/g;
102 my (@cigar_i)=$f[5]=~/(\d+)I/g; 104 my (@cigar_i)=$f[5]=~/(\d+)I/g;
114 # print "$clipseq\t$rightseq\t$overcoor\t"; 116 # print "$clipseq\t$rightseq\t$overcoor\t";
115 if ($overcoor > -1) {$overcoor += ($a[4] - 149);} 117 if ($overcoor > -1) {$overcoor += ($a[4] - 149);}
116 } 118 }
117 # print "\n"; 119 # print "\n";
118 } 120 }
119 elsif (($f[1]=~/R/)&&($z[0]=~/S/)) { 121 elsif ((($f[1] & 32) == 32)&&($z[0]=~/S/)) {
120 $coor=$f[3]; $strand="+"; 122 $coor=$f[3]; $strand="+";
121 my (@clipped)=$z[0]=~/(\d+)S/g; 123 my (@clipped)=$z[0]=~/(\d+)S/g;
122 my $cliplen=sum(@clipped); 124 my $cliplen=sum(@clipped);
123 # print "$f[0]\n"; 125 # print "$f[0]\n";
124 # print "$cliplen\t"; 126 # print "$cliplen\t";