Mercurial > repos > portiahollyoak > temp

comparison scripts/pickClippedFastq.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 36
37 my $lower=$a[1]-15; 37 my $lower=$a[1]-15;
38 my $upper=$a[2]+15; 38 my $upper=$a[2]+15;
39 if (($lower > 0)&&($upper > 0)) 39 if (($lower > 0)&&($upper > 0))
40 { 40 {
41 system("samtools view -hXf 0x2 $ARGV[0].sorted.bam $a[0]\:$lower\-$upper > temp.sam"); 41 system("samtools view -hf 0x2 $ARGV[0].sorted.bam $a[0]\:$lower\-$upper > temp.sam");
42 42
43 open in,"temp.sam"; 43 open in,"temp.sam";
44 my %pe1; 44 my %pe1;
45 my %pe2; 45 my %pe2;
46 while(<in>) 46 while(<in>)
47 { 47 {
48 chomp; 48 chomp;
49 my @f=split/\t/,$_,12; 49 my @f=split/\t/,$_,12;
50 ## read number 1 or 2 50 ## read number 1 or 2
51 my ($rnum)=$f[1]=~/(\d)$/; 51 #my ($rnum)=$f[1]=~/(\d)$/;
52 my $rnum=1;
53 if (($f[1] & 128) == 128) {$rnum=2;}
52 54
53 ## XT:A:* 55 ## XT:A:*
54 my ($xt)=$f[11]=~/XT:A:(.)/; 56 my ($xt)=$f[11]=~/XT:A:(.)/;
55 57
56 if ($f[5]=~/S/) { 58 if ($f[5]=~/S/) {
60 my $strand=""; 62 my $strand="";
61 my $final=""; 63 my $final="";
62 my $clipseq=""; 64 my $clipseq="";
63 my @z=split(/M/, $f[5]); 65 my @z=split(/M/, $f[5]);
64 66
65 if (($f[5]=~/S$/)&&($f[1]=~/r/)) 67 if (($f[5]=~/S$/)&&(($f[1] & 16) == 16))
66 { 68 {
67 my (@cigar_m)=$f[5]=~/(\d+)M/g; 69 my (@cigar_m)=$f[5]=~/(\d+)M/g;
68 my (@cigar_d)=$f[5]=~/(\d+)D/g; 70 my (@cigar_d)=$f[5]=~/(\d+)D/g;
69 my (@cigar_s)=$f[5]=~/(\d+)S/g; 71 my (@cigar_s)=$f[5]=~/(\d+)S/g;
70 my (@cigar_i)=$f[5]=~/(\d+)I/g; 72 my (@cigar_i)=$f[5]=~/(\d+)I/g;
77 if ($cliplen >= 15) { 79 if ($cliplen >= 15) {
78 $clipseq=substr($f[9], length($f[9])-$cliplen, $cliplen); 80 $clipseq=substr($f[9], length($f[9])-$cliplen, $cliplen);
79 } 81 }
80 } 82 }
81 83
82 elsif (($f[1]=~/R/)&&($z[0]=~/S/)) 84 elsif ((($f[1] & 32) == 32)&&($z[0]=~/S/))
83 { 85 {
84 $coor=$f[3]; $strand="+"; 86 $coor=$f[3]; $strand="+";
85 87
86 my (@clipped)=$z[0]=~/(\d+)S/g; 88 my (@clipped)=$z[0]=~/(\d+)S/g;
87 my $cliplen=sum(@clipped); 89 my $cliplen=sum(@clipped);
119 open in,"temp.sam"; 121 open in,"temp.sam";
120 while(<in>) 122 while(<in>)
121 { 123 {
122 chomp; 124 chomp;
123 my @f=split/\t/,$_,12; 125 my @f=split/\t/,$_,12;
124 my ($rnum)=$f[1]=~/(\d)$/; 126 #my ($rnum)=$f[1]=~/(\d)$/;
127 my $rnum=1;
128 if (($f[1] & 128) == 128) {$rnum=2;}
125 my ($xt)=$f[11]=~/XT:A:(.)/; 129 my ($xt)=$f[11]=~/XT:A:(.)/;
126 130
127 if ($f[5]=~/S/) { 131 if ($f[5]=~/S/) {
128 132
129 my $coor=-10; 133 my $coor=-10;
130 my $strand=""; 134 my $strand="";
131 my $final=""; 135 my $final="";
132 my $clipseq=""; 136 my $clipseq="";
133 my @z=split(/M/, $f[5]); 137 my @z=split(/M/, $f[5]);
134 138
135 if (($f[5]=~/S$/)&&($f[1]=~/r/)) 139 if (($f[5]=~/S$/)&&(($f[1] & 16) == 16))
136 { 140 {
137 my (@cigar_m)=$f[5]=~/(\d+)M/g; 141 my (@cigar_m)=$f[5]=~/(\d+)M/g;
138 my (@cigar_d)=$f[5]=~/(\d+)D/g; 142 my (@cigar_d)=$f[5]=~/(\d+)D/g;
139 my (@cigar_s)=$f[5]=~/(\d+)S/g; 143 my (@cigar_s)=$f[5]=~/(\d+)S/g;
140 my (@cigar_i)=$f[5]=~/(\d+)I/g; 144 my (@cigar_i)=$f[5]=~/(\d+)I/g;
147 if ($cliplen >= 6) { 151 if ($cliplen >= 6) {
148 $clipseq=substr($f[9], length($f[9])-$cliplen, $cliplen); 152 $clipseq=substr($f[9], length($f[9])-$cliplen, $cliplen);
149 } 153 }
150 } 154 }
151 155
152 elsif (($f[1]=~/R/)&&($z[0]=~/S/)) 156 elsif ((($f[1] & 32) == 32)&&($z[0]=~/S/))
153 { 157 {
154 $coor=$f[3]; $strand="+"; 158 $coor=$f[3]; $strand="+";
155 159
156 my (@clipped)=$z[0]=~/(\d+)S/g; 160 my (@clipped)=$z[0]=~/(\d+)S/g;
157 my $cliplen=sum(@clipped); 161 my $cliplen=sum(@clipped);