annotate bin/align.pm @ 22:59c728d31244 draft

Uploaded
author brasset_jensen
date Sun, 13 May 2018 15:44:09 -0400
parents 8ea13dab3435
children 36e45a7b3eff
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
1
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
1 package align;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
2
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
3 use strict;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
4 use warnings;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
5 use File::Basename;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
6 use String::Random;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
7
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
8 use FindBin;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
9 use lib $FindBin::Bin;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
10 use Rcall qw ( histogram );
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
11
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
12 use Exporter;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
13 our @ISA = qw( Exporter );
20
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
14 our @EXPORT = qw( &rpms_rpkm_te &BWA_call &to_build &get_unique &sam_sorted_bam &get_hash_alignment &sam_to_bam_bg &sam_count &sam_count_mis &rpms_rpkm &get_fastq_seq &extract_sam );
1
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
15
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
16 sub to_build
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
17 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
18 my ( $toBuildTabP, $log, $newdir ) = @_;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
19
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
20 foreach my $pairs ( @{ $toBuildTabP } )
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
21 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
22 if ( $pairs->[0] == 1 )
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
23 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
24 my $sym = $newdir.basename(${$pairs->[1]}).'_symlink.fa';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
25 symlink( ${$pairs->[1]}, $sym );
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
26 ${$pairs->[1]} = $sym;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
27 build_index ( ${$pairs->[1]}, $log );
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
28 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
29 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
30 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
31
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
32 sub build_index
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
33 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
34 my $to_index = shift;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
35 my $log = shift;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
36 my $index_log = $to_index.'_index.err';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
37 `bwa index $to_index 2> $index_log`;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
38 print $log "Creating index for $to_index\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
39 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
40
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
41 sub get_unique
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
42 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
43 my ( $sam, $s_uni, $out_prefix, $col_prefix, $details, $report ) = @_;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
44
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
45 my $fout = $col_prefix.'_all_mappers.fastq';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
46 my $funi = $col_prefix.'_unique_mappers.fastq';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
47 my $frej = $col_prefix.'_unmapped.fastq';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
48
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
49 my $repartition = $out_prefix.'distribution.txt';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
50 my $png_rep = $out_prefix.'distribution.png';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
51 my ( %duplicates, %genome_hits) ;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
52
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
53 #alignement to the first reference
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
54 my @return = sam_parse( $sam, $fout, $funi, $frej, $s_uni, \%duplicates, \%genome_hits, $report );
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
55 my $ref_fai = $return[4];
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
56 my $mappers = $return[5];
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
57 my $mappers_uni = $return[6];
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
58 my $size_mappedHashR = $return[7];
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
59
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
60 if ( $details == 1 )
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
61 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
62 #print number of duplicates and hits number
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
63 my ($pourcentage, $total) =(0,0);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
64
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
65 $total += $_ foreach values %{$size_mappedHashR};
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
66 open (my $rep, '>'.$repartition) || die "cannot create $repartition $!\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
67 print $rep "size\tnumber\tpercentage\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
68 foreach my $k (sort{$a cmp $b} keys (%{$size_mappedHashR}))
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
69 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
70 $pourcentage = 0;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
71 $pourcentage = $size_mappedHashR->{$k} / $total * 100 unless $total ==0;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
72
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
73 print $rep "$k\t$size_mappedHashR->{$k}\t";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
74 printf $rep "%.2f\n",$pourcentage;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
75 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
76
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
77 histogram($size_mappedHashR, $png_rep, $total);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
78
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
79
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
80 my $dup = $out_prefix.'dup_mapnum.txt';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
81 my $dup_u = $out_prefix .'dup_unique.txt';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
82 my $dup_r = $out_prefix .'dup_nonmapp.txt';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
83 open(my $tab,">".$dup) || die "cannot open output txt file\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
84 open(my $tab_r,">".$dup_r) || die "cannot open output txt file\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
85 open(my $tab_u,">".$dup_u) || die "cannot open output txt file\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
86 print $tab "sequence\tcount\tmapnum\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
87 print $tab_u "sequence\tcount\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
88 print $tab_r "sequence\tcount\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
89 foreach my $k (sort {$duplicates{$b} <=> $duplicates{$a}}keys %duplicates)
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
90 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
91 $duplicates{$k} = 0 unless exists($duplicates{$k});
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
92 $genome_hits{$k} = 0 unless exists($genome_hits{$k});
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
93 if ($genome_hits{$k} != 0) { print $tab $k."\t".$duplicates{$k}."\t".$genome_hits{$k}."\n"; }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
94 else {print $tab_r $k."\t".$duplicates{$k}."\n";}
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
95 if ($genome_hits{$k} == 1) { print $tab_u $k."\t".$duplicates{$k}."\n"; }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
96 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
97 close $dup; close $dup_r; close $dup_u;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
98 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
99 return ( $ref_fai, $mappers, $mappers_uni );
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
100 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
101
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
102 sub sam_parse
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
103 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
104 my ( $sam, $fastq_accepted, $fastq_accepted_unique, $fastq_rejected, $sam_unique, $duplicate_hashR, $best_hit_number_hashR, $report ) = @_ ;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
105 my ($reads, $mappers, $mappersUnique, @garbage, %size_num, %size_num_spe, %number, %numberSens, %numberReverse, %unique_number, %numberNM, %numberM, %size);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
106 $mappers = $mappersUnique = $reads = 0;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
107
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
108 open my $fic, '<', $sam || die "cannot open $sam $!\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
109 open my $accepted, '>', $fastq_accepted || die "cannot create $fastq_accepted $! \n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
110 open my $unique, '>', $fastq_accepted_unique || die "cannot create $fastq_accepted_unique $! \n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
111 open my $rejected, '>', $fastq_rejected || die "cannot create $fastq_rejected $! \n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
112 open my $sam_uni, '>', $sam_unique || die "cannot create $sam_unique $! \n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
113 my $sequence = '';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
114 while(<$fic>)
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
115 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
116 chomp $_;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
117 if ($_ =~ /^\@[A-Za-z][A-Za-z](\t[A-Za-z][A-Za-z0-9]:[ -~]+)+$/ || $_ =~ /^\@CO\t.*/ )
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
118 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
119 if ($_ =~ /\@SQ\tSN:(.*)\tLN:(\d*)/)
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
120 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
121 $size{$1} = $2;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
122 $unique_number{$1} = 0;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
123 $number{$1} = 0;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
124 $numberNM{$1} = 0;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
125 $numberM{$1} = 0;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
126 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
127 print $sam_uni $_."\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
128 next;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
129 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
130 $reads++;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
131 my @line = split (/\t/,$_);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
132 $sequence = $line[9];
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
133 if ($line[1] & 16)
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
134 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
135 $sequence =reverse($sequence);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
136 $sequence =~ tr/atgcuATGCU/tacgaTACGA/;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
137 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
138 if ($line[1] == 16 || $line[1] == 0)
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
139 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
140 my $len = length($sequence);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
141 $size_num{$len} ++;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
142 $size_num_spe{$line[2]}{$len}++;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
143 $mappers ++;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
144
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
145 ${$best_hit_number_hashR}{$sequence} = $1 if ($line[13] =~ /X0:i:(\d*)/ || $line[14] =~/X0:i:(\d*)/ );
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
146 ${$duplicate_hashR}{$sequence}++;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
147 $number{$line[2]}++;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
148 $numberSens{$line[2]}++ if $line[1] == 0 ;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
149 $numberReverse{$line[2]}++ if $line[1] == 16 ;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
150 print $accepted "\@".$line[0]."\n".$sequence."\n+\n".$line[10]."\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
151
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
152 if ($line[11] eq "XT:A:U")
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
153 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
154 $unique_number{$line[2]}++;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
155 $mappersUnique++;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
156 print $unique "\@".$line[0]."\n".$sequence."\n+\n".$line[10]."\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
157 print $sam_uni $_."\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
158 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
159 if ($_ =~ /.*XM:i:(\d+).*/)
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
160 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
161 if ($1 == 0){$numberNM{$line[2]}++;}else{$numberM{$line[2]}++;}
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
162 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
163 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
164 else
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
165 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
166 ${$best_hit_number_hashR}{$sequence} = 0;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
167 ${$duplicate_hashR}{$sequence}++;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
168 print $rejected "\@".$line[0]."\n".$sequence."\n+\n".$line[10]."\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
169 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
170 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
171 close $fic; close $accepted; close $unique; close $rejected; close $sam_uni;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
172
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
173 print $report "Parsing $sam file\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
174 print $report "\treads: $reads\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
175 print $report "\tmappers: $mappers\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
176 print $report "\tunique mappers: $mappersUnique\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
177 print $report "-----------------------------\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
178 return (\%number, \%unique_number, \%numberSens, \%numberReverse, \%size, $mappers, $mappersUnique, \%size_num, \%size_num_spe, \%numberNM, \%numberM );
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
179 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
180
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
181 sub get_hash_alignment
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
182 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
183 my ($index, $mismatches, $accept, $reject, $outA, $outR, $fastq, $number_of_cpus, $name, $sam, $report, $fai_f) = @_ ;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
184 my ($reads, $mappers, $unmapped) = (0,0,0);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
185 my $accep_unique;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
186 BWA_call ( $index, $fastq, $sam, $mismatches, $number_of_cpus, $report );
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
187
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
188 open my $fic, '<', $sam || die "cannot open $sam $!\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
189 open my $accepted, '>', $outA || die "cannot open $outA\n" if $accept == 1;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
190 open my $rejected, '>', $outR || die "cannot open $outR\n" if $reject == 1;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
191 open my $fai, '>', $fai_f || die "cannot open $fai_f\n" if $fai_f;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
192 #if ($name eq "snRNAs") {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
193 # open ( $accep_unique, ">".$1."-unique.fastq") if $outR =~ /(.*)\.fastq/;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
194 #}
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
195 my $sequence = '';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
196 while(<$fic>)
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
197 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
198 chomp $_;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
199 if( $_ =~ /^\@[A-Za-z][A-Za-z](\t[A-Za-z][A-Za-z0-9]:[ -~]+)+$/ || $_ =~ /^\@CO\t.*/ )
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
200 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
201 if ($fai_f && $_ =~ /\@SQ\tSN:(.*)\tLN:(\d*)/)
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
202 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
203 print $fai $1."\t".$2."\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
204 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
205 next;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
206 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
207 $reads++;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
208 my @line = split (/\t/,$_);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
209 $sequence = $line[9];
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
210 if ($line[1] & 16)
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
211 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
212 $sequence =reverse($sequence);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
213 $sequence =~ tr/atgcuATGCU/tacgaTACGA/;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
214 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
215 if ($line[1] & 16 || $line[1] == 0)
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
216 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
217 $mappers ++;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
218 if ($accept == 1 )
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
219 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
220 print $accepted "\@".$line[0]."\n".$sequence."\n+\n".$line[10]."\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
221 # print $accep_unique "\@".$line[0]."\n".$sequence."\n+\n".$line[10]."\n" if ($name eq "snRNAs" && $line[11] eq "XT:A:U");
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
222 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
223 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
224 else
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
225 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
226 print $rejected "\@".$line[0]."\n".$sequence."\n+\n".$line[10]."\n" if $reject == 1;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
227 $unmapped++;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
228 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
229 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
230 # close $accep_unique if ($name eq "bonafide_reads");
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
231 close $fic;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
232 close $accepted if $accept == 1;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
233 close $rejected if $reject ==1;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
234 close $fai if $fai_f;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
235 print $report "\treads: $reads\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
236 print $report "\tmappers: $mappers\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
237 print $report "\tunmapped: $unmapped\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
238 print $report "-----------------------------\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
239 return ($mappers, $unmapped);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
240 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
241
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
242 sub sam_count
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
243 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
244 my $sam = shift;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
245 my ( %number, %size );
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
246
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
247 open my $fic, '<', $sam || die "cannot open $sam file $!\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
248 while(<$fic>)
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
249 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
250 chomp $_;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
251 if ($_ =~ /^\@[A-Za-z][A-Za-z](\t[A-Za-z][A-Za-z0-9]:[ -~]+)+$/ || $_ =~ /^\@CO\t.*/ )
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
252 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
253 if ($_ =~ /\@SQ\tSN:(.*)\tLN:(\d*)/)
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
254 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
255 $size{$1} = $2;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
256 $number{$1} = 0;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
257 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
258 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
259 else
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
260 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
261 my @line = split (/\t/,$_);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
262 if ( $line[1] & 16 || $line[1] == 0 )
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
263 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
264 $number{$line[2]}++;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
265 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
266 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
267 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
268 close $fic;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
269 return ( \%number, \%size );
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
270 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
271
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
272 sub sam_count_mis
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
273 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
274
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
275 my $sam = shift;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
276 my ( %number, %numberNM, %numberM, %size);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
277
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
278 open my $fic, '<', $sam || die "cannot open $sam file $!\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
279 while(<$fic>)
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
280 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
281 chomp $_;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
282 if ($_ =~ /^\@[A-Za-z][A-Za-z](\t[A-Za-z][A-Za-z0-9]:[ -~]+)+$/ || $_ =~ /^\@CO\t.*/ )
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
283 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
284 if ($_ =~ /\@SQ\tSN:(.*)\tLN:(\d*)/)
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
285 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
286 $size{$1} = $2;
20
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
287 $number{$1} = [0,0,0,0,0];
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
288 $numberNM{$1} = [0,0,0,0,0];
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
289 $numberM{$1} = [0,0,0,0,0];
1
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
290 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
291 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
292 else
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
293 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
294 my @line = split (/\t/,$_);
20
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
295 my @seq = split //, $line[9];
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
296 if ( $line[1] == 16 || $line[1] == 0 )
1
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
297 {
20
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
298 $number{ $line[2] }->[0]++;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
299 if ($line[1] == 0)
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
300 {
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
301 $number{$line[2]}->[1]++;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
302 $number{$line[2]}->[3]++ if $seq[0] eq 'T';
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
303 }
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
304 else
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
305 {
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
306 $number{$line[2]}->[2]++ if $seq[9] eq 'A';
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
307 }
1
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
308 if ($_ =~ /.*XM:i:(\d+).*/)
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
309 {
20
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
310 if ( $1 == 0 )
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
311 {
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
312 $numberNM{$line[2]}->[0]++;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
313 if ($line[1] == 0)
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
314 {
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
315 $numberNM{$line[2]}->[1]++;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
316 $numberNM{$line[2]}->[3]++ if $seq[0] eq 'T';
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
317 }
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
318 else
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
319 {
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
320 $numberNM{$line[2]}->[2]++;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
321 $numberNM{$line[2]}->[4]++ if $seq[9] eq 'A';
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
322 }
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
323 }
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
324 else
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
325 {
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
326 $numberM{$line[2]}->[0]++;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
327 if ($line[1] == 0)
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
328 {
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
329 $numberM{$line[2]}->[1]++;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
330 $numberM{$line[2]}->[3]++ if $seq[0] eq 'T';
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
331 }
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
332 else
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
333 {
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
334 $numberM{$line[2]}->[2]++;
22
59c728d31244 Uploaded
brasset_jensen
parents: 20
diff changeset
335 $numberM{$line[2]}->[4]++ if $seq[9] eq 'A';
20
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
336 }
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
337 }
1
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
338 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
339 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
340 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
341 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
342 return (\%number, \%size, \%numberNM, \%numberM );
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
343 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
344
20
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
345 sub rpms_rpkm_te
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
346 {
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
347 my ( $counthashR, $sizehashR, $mapped, $out_file, $piRNA_number, $miRNA_number, $bonafide_number ) =@_;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
348 open(my $out, ">".$out_file) || die "cannot open normalized file $! \n";
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
349 print $out "ID\treads counts\tRPKM\tsens reads counts\treverse reads counts\t% of sens 1U\t% of reverse 10A";
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
350 print $out "\tper million of piRNAs" if ($piRNA_number != 0);
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
351 print $out "\tper million of miRNAs" if ($miRNA_number != 0);
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
352 print $out "\tper million of bonafide reads" if ($bonafide_number != 0);
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
353 print $out "\n";
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
354 foreach my $k ( sort keys %{$counthashR} )
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
355 {
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
356 my ($rpkm, $pirna, $mirna, $bonafide) = (0,0,0,0);
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
357
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
358 $rpkm = ( $counthashR->{$k}->[0] * 1000000000) / ( $sizehashR->{$k} * $mapped) if ( $sizehashR->{$k} * $mapped) != 0 ;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
359 print $out $k."\t".$counthashR->{$k}->[0]."\t"; printf $out "%.2f",$rpkm;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
360 print $out "\t".$counthashR->{$k}->[1]."\t".$counthashR->{$k}->[2] ;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
361 my $U1 = 0;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
362 $U1 = $counthashR->{$k}->[3] / $counthashR->{$k}->[1] * 100 if $counthashR->{$k}->[1] != 0;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
363 my $A10 = 0;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
364 $A10 = $counthashR->{$k}->[4] / $counthashR->{$k}->[2] * 100 if $counthashR->{$k}->[2] != 0;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
365
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
366 print $out "\t".$U1."\t".$A10;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
367
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
368 if ($piRNA_number != 0 )
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
369 {
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
370 $pirna = ( $counthashR->{$k}->[0] * 1000000) / $piRNA_number;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
371 printf $out "\t%.2f",$pirna;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
372 }
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
373 if ($miRNA_number != 0 )
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
374 {
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
375 $mirna = ( $counthashR->{$k}->[0] * 1000000) / $miRNA_number;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
376 printf $out "\t%.2f",$mirna;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
377 }
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
378 if ($bonafide_number != 0 )
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
379 {
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
380 $bonafide = ( $counthashR->{$k}->[0] * 1000000) / $bonafide_number;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
381 printf $out "\t%.2f",$bonafide;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
382 }
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
383 print $out "\n";
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
384 }
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
385 close $out;
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
386 }
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
387
8ea13dab3435 Uploaded
brasset_jensen
parents: 1
diff changeset
388
1
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
389 sub sam_to_bam_bg
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
390 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
391 my ( $sam, $scale, $number_of_cpus ) = @_;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
392 my ( $bam_sorted, $bedgraphM, $bedgraphP, $view_err, $sort_err ) = ( '', '', '', '', '' );
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
393 if ( $sam =~ /(.*?).sam$/ )
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
394 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
395 $bam_sorted = $1.'_sorted.bam';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
396 $bedgraphP= $1.'_plus.bedgraph';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
397 $bedgraphM = $1.'_minus.bedgraph';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
398 $view_err = $1.'_view.err';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
399 $sort_err = $1.'_sort.err';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
400 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
401 `samtools view -Shb --threads $number_of_cpus $sam 2> $view_err | samtools sort -O BAM --threads $number_of_cpus /dev/stdin 2> $sort_err > $bam_sorted`;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
402 `bedtools genomecov -scale $scale -strand + -bga -ibam $bam_sorted > $bedgraphP`;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
403 `bedtools genomecov -scale $scale -strand - -bga -ibam $bam_sorted > $bedgraphM`;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
404 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
405
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
406 sub sam_sorted_bam
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
407 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
408 my ( $sam, $number_of_cpus ) = @_;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
409 my ( $bam_sorted, $view_err, $sort_err ) = ( '', '', '' );
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
410 if ( $sam =~ /(.*?).sam$/ )
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
411 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
412 $bam_sorted = $1.'_sorted.bam';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
413 $view_err = $1.'_view.err';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
414 $sort_err = $1.'_sort.err';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
415
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
416 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
417 `samtools view -Shb --threads $number_of_cpus $sam 2> $view_err | samtools sort -O BAM --threads $number_of_cpus /dev/stdin 2> $sort_err > $bam_sorted`;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
418 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
419
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
420 sub BWA_call
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
421 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
422 my ( $index, $fastq, $sam, $mismatches, $number_of_cpus, $report ) = @_;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
423 my ( $aln_err, $samse_err, $seq_num ) = ( $sam.'_aln.err', $sam.'_samse.err', 0 );
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
424 print $report "-----------------------------\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
425 print $report "bwa aln -t $number_of_cpus -n $mismatches $index $fastq 2> $aln_err | bwa samse $index /dev/stdin $fastq 2> $samse_err > $sam\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
426 `bwa aln -t $number_of_cpus -n $mismatches $index $fastq 2> $aln_err | bwa samse $index /dev/stdin $fastq 2> $samse_err > $sam `;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
427 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
428
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
429 sub rpms_rpkm
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
430 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
431 my ( $counthashR, $sizehashR, $mapped, $out_file, $piRNA_number, $miRNA_number, $bonafide_number ) =@_;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
432 open(my $out, ">".$out_file) || die "cannot open normalized file $! \n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
433 print $out "ID\treads counts\tRPKM";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
434 print $out "\tper million of piRNAs" if ($piRNA_number != 0);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
435 print $out "\tper million of miRNAs" if ($miRNA_number != 0);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
436 print $out "\tper million of bonafide reads" if ($bonafide_number != 0);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
437 print $out "\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
438 foreach my $k ( sort keys %{$counthashR} )
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
439 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
440 my ($rpkm, $pirna, $mirna, $bonafide) = (0,0,0,0);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
441
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
442 $rpkm = ( $counthashR->{$k} * 1000000000) / ( $sizehashR->{$k} * $mapped) if ( $sizehashR->{$k} * $mapped) != 0 ;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
443 print $out $k."\t".$counthashR->{$k}."\t"; printf $out "%.2f",$rpkm;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
444
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
445 if ($piRNA_number != 0 )
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
446 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
447 $pirna = ( $counthashR->{$k} * 1000000) / $piRNA_number;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
448 printf $out "\t%.2f",$pirna;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
449 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
450 if ($miRNA_number != 0 )
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
451 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
452 $mirna = ( $counthashR->{$k} * 1000000) / $miRNA_number;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
453 printf $out "\t%.2f",$mirna;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
454 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
455 if ($bonafide_number != 0 )
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
456 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
457 $bonafide = ( $counthashR->{$k} * 1000000) / $bonafide_number;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
458 printf $out "\t%.2f",$bonafide;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
459 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
460 print $out "\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
461 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
462 close $out;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
463 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
464
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
465 sub extract_sam
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
466 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
467 my ( $hashRef, $sam_in, $sam_out, $sam_uni_out, $fastq_out, $fastq_uni_out ) = @_;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
468
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
469 open my $s_in, '<', $sam_in || die "cannot open $sam_in file $!\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
470
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
471 open my $f_out, '>', $fastq_out || die "cannot create $fastq_out $!\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
472 open my $f_uni_out, '>', $fastq_uni_out || die "cannot create $fastq_uni_out $!\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
473
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
474 open my $s_out, '>', $sam_out || die "cannot create $sam_out file $!\n" if defined ($hashRef);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
475 open my $s_uni_out, '>', $sam_uni_out || die "cannot create $sam_uni_out file $!\n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
476
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
477 my $sequence = '';
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
478 while(<$s_in>)
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
479 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
480 if ($_ =~ /^\@[A-Za-z][A-Za-z](\t[A-Za-z][A-Za-z0-9]:[ -~]+)+$/ || $_ =~ /^\@CO\t.*/ )
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
481 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
482 print $s_out $_ if defined ($hashRef);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
483 print $s_uni_out $_;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
484 next;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
485 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
486 my @line = split (/\t/,$_);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
487 $sequence = $line[0];
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
488 if ( (! defined ($hashRef) )|| ( exists $hashRef->{$sequence} && $hashRef->{$sequence} == 1 ) )
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
489 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
490 my $arn = $line[9];
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
491 if ($line[1] & 16)
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
492 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
493 $arn =reverse($arn);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
494 $arn =~ tr/atgcuATGCU/tacgaTACGA/;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
495 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
496
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
497 if ( ( $line[1] == 16 || $line[1] == 0 ) )
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
498 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
499 print $f_out "\@".$line[0]."\n".$arn."\n+\n".$line[10]."\n" ;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
500 print $s_out $_ if defined ($hashRef);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
501 if ( $line[11] eq "XT:A:U" )
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
502 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
503 print $f_uni_out "\@".$line[0]."\n".$arn."\n+\n".$line[10]."\n" ;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
504 print $s_uni_out $_ ;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
505 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
506 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
507 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
508 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
509 close $s_in; close $s_out if defined ($hashRef);
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
510 close $s_uni_out; close $f_out; close $f_uni_out;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
511 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
512
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
513 sub get_fastq_seq
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
514 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
515 my $fastq = shift;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
516 my %hash; my $cmp = 0;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
517
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
518 open my $fic, '<', $fastq || die "cannot open input file $! \n";
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
519 while(<$fic>)
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
520 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
521 chomp $_;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
522 $cmp++;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
523 if ($cmp % 4 == 1)
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
524 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
525 die "file do not contain a @ at line $cmp\n" unless ($_ =~ /^\@/ );
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
526 if ($_ =~ /^\@(.*)\s.*/) { $hash{$1} = 1;}
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
527 elsif ($_ =~ /^\@(.*)/) { $hash{$1} = 1;}
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
528 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
529 elsif ($cmp % 4 == 3 )
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
530 {
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
531 die "file do not contain a + at line $cmp\n" unless $_ =~ /^\+/;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
532 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
533 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
534 close $fic;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
535 return \%hash;
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
536 }
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
537
1df6aaac800e Deleted selected files
brasset_jensen
parents:
diff changeset
538 1;