|
11
|
1 #!/usr/bin/perl -w
|
|
|
2 my $version=1.00;
|
|
|
3 use strict;
|
|
|
4 use warnings;
|
|
|
5 use Getopt::Long;
|
|
|
6 use Getopt::Std;
|
|
|
7 use threads;
|
|
|
8 use threads::shared;
|
|
|
9 use Parallel::ForkManager;
|
|
|
10 use lib '/leofs/biotrans/chentt/perl_module/';
|
|
|
11 #perl ../siRNA.pl -i config -g /leofs/biotrans/projects/rice/smallRNA/sRNA_package/bin/test/ref/genome.fa -f /share_bio/hs4/disk3-4/Reference/Plants/Rice_TIGR/Reference/TIGR/version_6.1/all.dir/all.gff3 -path /leofs/biotrans/projects/rice/smallRNA/sRNA_package/bin/ -o /leofs/biotrans/projects/rice/smallRNA/sRNA_package/bin/test -t 3 -rfam /leofs/biotrans/projects/rice/smallRNA/sRNA_package/bin/test/ref/Rfam.fasta -idx /leofs/biotrans/projects/rice/smallRNA/sRNA_package/bin/test/ref/genome -idx2 /leofs/biotrans/projects/rice/smallRNA/sRNA_package/bin/test/ref/rfam -deg deg -n 25 -nat class/nat_1 -repeat class/repeat_1 -cen centromere_TIGR.txt -format fastq
|
|
|
12 print "
|
|
|
13 #####################################
|
|
|
14 # #
|
|
|
15 # sRNA cluster #
|
|
|
16 # #
|
|
|
17 #####################################
|
|
|
18 ";
|
|
|
19 ###########################################################################################
|
|
|
20 my $usage="$0
|
|
|
21 Options:
|
|
|
22 -i input file# raw data file
|
|
|
23 -tag string #raw data sample name
|
|
|
24 -g genome file
|
|
|
25 -f gff file
|
|
|
26
|
|
|
27 -o workdir file
|
|
|
28 -path script path
|
|
|
29 -t int, number of threads [1]
|
|
|
30 -format fastq, fq, fasta or fa
|
|
|
31 -idx string, genome file index, file-prefix #(must be indexed by bowtie-build) The parameter
|
|
|
32 string must be the prefix of the bowtie index. For instance, if
|
|
|
33 the first indexed file is called 'h_sapiens_37_asm.1.ebwt' then
|
|
|
34 the prefix is 'h_sapiens_37_asm'.##can be null
|
|
|
35 -mis int number of allowed mismatches when mapping reads to genome, default 0
|
|
|
36 -rfam string, input file# rfam database file.
|
|
|
37 -idx2 string, rfam file index, file-prefix #(must be indexed by bowtie-build) The parameter
|
|
|
38 string must be the prefix of the bowtie index. For instance, if
|
|
|
39 the first indexed file is called 'h_sapiens_37_asm.1.ebwt' then
|
|
|
40 the prefix is 'h_sapiens_37_asm'.##can be null
|
|
|
41
|
|
|
42 -v int report end-to-end hits w/ <=v mismatches; ignore qualities,default 0; used in rfam alignment
|
|
|
43
|
|
|
44 -a string, ADAPTER string. default is ATCTCGTATG.
|
|
|
45 -n int max hits number,default 25; used in genome alignment
|
|
|
46 -d int distance of tag to merged a cluster; default 100
|
|
|
47 -p cluster method F :conventional default is F
|
|
|
48 T :NIBLES
|
|
|
49 -l int the length of the upstream and downstream,default 1000;used in position annotate
|
|
|
50
|
|
|
51 -nat natural antisense transcripts file
|
|
|
52 -repeat repeat information file out of Repeatmasker
|
|
|
53 -deg file config of de sample
|
|
|
54 -cen centromere file input
|
|
|
55 -span plot span, default 50000
|
|
|
56 ";
|
|
|
57
|
|
|
58 my %options;
|
|
|
59 GetOptions(\%options,"i:s@","tag:s@","g=s","f=s","o=s","a:s","path:s","p=s","format=s","nat:s","repeat:s","deg:s","n:i","mis:i","rfam:s","t:i","v:i","d:i","l:i","idx:s","idx2:s","cen:s","span:s","h");
|
|
|
60 #print help if that option is used
|
|
|
61 if($options{h}){die $usage;}
|
|
|
62
|
|
|
63 my @filein=@{$options{'i'}};
|
|
|
64 my @mark=@{$options{'tag'}};
|
|
|
65
|
|
|
66 #my $config=$options{'i'};
|
|
|
67 my $genome_fa=$options{'g'};
|
|
|
68 my $gff=$options{'f'};
|
|
|
69
|
|
|
70
|
|
|
71 ##########################################################################################
|
|
|
72 my $predir=`pwd`;
|
|
|
73 chomp $predir;
|
|
|
74 my $workdir=defined($options{'o'}) ? $options{'o'}:$predir;
|
|
|
75
|
|
|
76 my $path=$options{'path'};
|
|
|
77
|
|
|
78 my $t=defined($options{'t'})? $options{'t'}:1; #threads number
|
|
|
79
|
|
|
80 my $mis=defined $options{'mis'} ? $options{'mis'}:0;
|
|
|
81
|
|
|
82 my $mis_rfam=defined $options{'v'} ? $options{'v'}:0;
|
|
|
83
|
|
|
84 my $hit=defined $options{'n'}?$options{'n'}:25;
|
|
|
85
|
|
|
86 my $distance_of_merged_tag=defined $options{'d'} ? $options{'d'}:100;
|
|
|
87
|
|
|
88 my $up_down_dis=defined $options{'l'} ?$options{'l'}:1000;
|
|
|
89
|
|
|
90 my $cluster_mothod=defined $options{'p'}?$options{'p'}:"F";
|
|
|
91
|
|
|
92 my $format=$options{'format'};
|
|
|
93 #if ($format ne "fastq" && $format ne "fq" && $format ne "fasta" && $format ne "fa") {
|
|
|
94 # die "Parameter \"-format\" is error! Parameter is fastq, fq, fasta or fa\n";
|
|
|
95 #}
|
|
|
96
|
|
|
97 my $adpter="ATCTCGTATG"; #adapter
|
|
|
98 if (defined $options{'a'}) {$adpter=$options{'a'};}
|
|
|
99
|
|
|
100
|
|
|
101 my $phred_qv=64;
|
|
|
102 my $sample_number;
|
|
|
103 my ($dir,$dir_tmp);
|
|
|
104 ################################ MAIN ##################################################
|
|
|
105 print "\ncluster program start:";
|
|
|
106 my $time=Time();
|
|
|
107 make_dir_tmp();
|
|
|
108
|
|
|
109 my @clip;
|
|
|
110 my $mark;
|
|
|
111 my $sample_mark;
|
|
|
112
|
|
|
113 my $config=$dir."/input_config";
|
|
|
114 open CONFIG,">$config";
|
|
|
115 for (my $i=0;$i<@filein;$i++) {
|
|
|
116 print CONFIG $filein[$i],"\t",$mark[$i],"\n";
|
|
|
117 }
|
|
|
118 close CONFIG;
|
|
|
119 if (@filein != @mark) {
|
|
|
120 die "Maybe config file have some wrong!!!\n";
|
|
|
121 }
|
|
|
122 $sample_number=@mark;
|
|
|
123 $mark=join "\t",@mark;
|
|
|
124 $sample_mark=join "\#",@mark;
|
|
|
125
|
|
|
126
|
|
|
127 #read_config();
|
|
|
128
|
|
|
129 trim_adapter_and_filter();
|
|
|
130
|
|
|
131 my $filter_out=$dir."preProcess\/"."collapse_reads_out.fa";## raw clean data
|
|
|
132 my $data2=$filter_out; ### mirbase not mapped reads
|
|
|
133 my $data3=$dir."\/rfam_match\/rfam_not_mapped\.fa"; ### rfam not mapped reads
|
|
|
134 my $bed=$dir."cluster\/"."sample\.bed";
|
|
|
135 my $read=$dir."cluster\/"."sample_reads\.cluster";
|
|
|
136 my $read_txt=$dir."cluster\/"."cluster\.txt";
|
|
|
137 my $rpkm=$dir."cluster\/"."sample_rpkm\.cluster";
|
|
|
138 my $preprocess;
|
|
|
139 my $cluster_file;
|
|
|
140 my $annotate_dir;
|
|
|
141 my $deg_dir;
|
|
|
142 my %id;
|
|
|
143 for (my $i=0;$i<@mark ;$i++) {
|
|
|
144 $id{$mark[$i]}=$i+4;
|
|
|
145 }
|
|
18
|
146
|
|
|
147
|
|
11
|
148 group_and_filter(); #collapse reads to tags
|
|
|
149
|
|
|
150 rfam();
|
|
|
151
|
|
|
152 my @map_read;
|
|
|
153 my $map_tag=0;
|
|
|
154 genome();
|
|
|
155
|
|
|
156 bwt2bed();
|
|
|
157
|
|
|
158 cluster();
|
|
|
159
|
|
|
160 quantify();
|
|
|
161
|
|
|
162 phase();
|
|
|
163
|
|
|
164 if (defined($options{'nat'})&&defined($options{'repeat'})) {
|
|
|
165 class();
|
|
|
166 }
|
|
|
167 else{
|
|
|
168 get_genelist();
|
|
|
169 }
|
|
|
170
|
|
|
171 annotate();
|
|
|
172
|
|
|
173 genome_length();
|
|
|
174
|
|
|
175 plot();
|
|
|
176
|
|
|
177 my @pairdir;
|
|
|
178 if (defined($options{'deg'})) {
|
|
|
179 dec();
|
|
|
180 infor_merge();
|
|
|
181 }
|
|
|
182 else{infor_merge_no_dec()}
|
|
|
183 html();
|
|
|
184 print "\ncluster program end:";
|
|
|
185 Time();
|
|
|
186 ############################sub program###################################################
|
|
|
187 sub make_dir_tmp{
|
|
|
188
|
|
|
189 #make temporary directory
|
|
|
190 if(not -d "$workdir\/cluster_runs"){
|
|
|
191 mkdir("$workdir\/cluster_runs");
|
|
|
192 mkdir("$workdir\/cluster_runs\/ref\/");
|
|
|
193 }
|
|
|
194
|
|
|
195 $dir="$workdir\/cluster_runs\/";
|
|
|
196 #print STDERR "mkdir $dir\n\n";
|
|
|
197 return;
|
|
|
198 }
|
|
|
199
|
|
|
200 #sub read_config{
|
|
|
201 # open IN,"<$config";
|
|
|
202 # while (my $aline=<IN>) {
|
|
|
203 # chomp $aline;
|
|
|
204 # my @tmp=split/\t/,$aline;
|
|
|
205 # push @filein,$tmp[0];
|
|
|
206 # push @mark,$tmp[1];
|
|
|
207 # }
|
|
|
208 # close IN;
|
|
|
209 # if (@filein != @mark) {
|
|
|
210 # die "Maybe config file have some wrong!!!\n";
|
|
|
211 # }
|
|
|
212 # $sample_number=@mark;
|
|
|
213 # $mark=join "\t",@mark;
|
|
|
214 # $sample_mark=join "\#",@mark;
|
|
|
215 #}
|
|
|
216
|
|
|
217
|
|
|
218 sub trim_adapter_and_filter{
|
|
|
219 my $time=time();
|
|
|
220 $preprocess=$dir."preProcess/";
|
|
|
221 mkdir $preprocess;
|
|
|
222 my $can_use_threads = eval 'use threads; 1';
|
|
|
223 if ($can_use_threads) {
|
|
|
224 # Do processing using threads
|
|
|
225 my @filein1=@filein; my @mark1=@mark;
|
|
|
226 while (@filein1>0) {
|
|
|
227 my @thrs; my @res;
|
|
|
228 for (my $i=0;$i<$t ;$i++) {
|
|
|
229 last if(@filein1==0);
|
|
|
230 my $in=shift @filein1;
|
|
|
231 my $out=shift @mark1;
|
|
|
232 push @clip,$dir."preProcess\/$out\_clip\.fq";
|
|
|
233 $thrs[$i]=threads->create(\&clips,$in,$out);
|
|
|
234 }
|
|
|
235 for (my $i=0;$i<@thrs;$i++) {
|
|
|
236 $res[$i]=$thrs[$i]->join();
|
|
|
237 }
|
|
|
238 }
|
|
|
239 }
|
|
|
240 else {
|
|
|
241 # Do not processing using threads
|
|
|
242 for (my $i=0;$i<@filein ;$i++) {
|
|
|
243 my $in=$filein[$i];
|
|
|
244 my $out=$mark[$i];
|
|
|
245 push @clip,$dir."preProcess\/$out\_clip\.fq";
|
|
|
246 &clips($in,$out);
|
|
|
247 }
|
|
|
248 }
|
|
|
249 }
|
|
|
250
|
|
|
251 sub clips{
|
|
|
252 my ($filein,$fileout)=@_;
|
|
|
253 my $adapter=$dir."preProcess\/$fileout\_clip\.fq";
|
|
|
254 if($format eq "fq" || $format eq "fastq"){
|
|
|
255 my $clip=`fastx_clipper -a $adpter -M 6 -Q $phred_qv -i $filein -o $adapter`;
|
|
|
256 }
|
|
|
257 if($format eq "fa" || $format eq "fasta"){
|
|
|
258 my $clip=`fastx_clipper -a $adpter -M 6 -i $filein -o $adapter`;
|
|
|
259 }
|
|
|
260 #my $clean=$dir."preProcess\/$fileout\_clean.fq";
|
|
|
261 #my $filter=`filterReadsByLength.pl -i $adapter -o $clean -min 18 -max 40 `;
|
|
|
262 return $fileout;
|
|
|
263 }
|
|
|
264
|
|
|
265 sub group_and_filter{
|
|
|
266 #my ($ins,$data)=@_;
|
|
|
267 my @ins=@clip;
|
|
|
268 my $str="";
|
|
|
269 my $group_out_file=$dir."preProcess\/"."collapse_reads.fa";
|
|
|
270 #print "$$ins[0]\t$$ins[0]\n";
|
|
|
271 for (my $i=0;$i<@clip;$i++) {
|
|
|
272 $str .="-i $clip[$i] ";
|
|
|
273 #print "$$ins[$i]\n";
|
|
|
274 }
|
|
|
275 my $group=`perl $path\/collapseReads2Tags.pl $str -mark seq -o $group_out_file -format $format`;
|
|
|
276 print "perl $path\/collapseReads2Tags.pl $str -mark seq -o $group_out_file -format $format\n\n";
|
|
|
277
|
|
|
278 my $l_out=$dir."preProcess\/"."collapse_reads_18-40.fa";
|
|
|
279 my $length_f=`perl $path\/filterReadsByLength_1.pl -i $group_out_file -o $l_out -min 18 -max 40 -mark $sample_mark`;
|
|
|
280 print "perl $path\/filterReadsByLength_1.pl -i $group_out_file -o $l_out -min 18 -max 40 -mark $sample_mark\n\n";
|
|
|
281 my $cout_f=`perl $path\/filterReadsByCount.pl -i $l_out -o $filter_out -mark $sample_mark`;
|
|
|
282 print "perl $path\/filterReadsByCount.pl -i $l_out -o $filter_out -mark $sample_mark\n\n";
|
|
|
283 return 0;
|
|
|
284 }
|
|
|
285
|
|
|
286 sub rfam{
|
|
|
287 if (defined $options{'idx2'}) {
|
|
|
288 system("perl $path\/rfam.pl -i $data2 -ref $options{rfam} -v $mis_rfam -p $t -o $dir -index $options{idx2}");
|
|
|
289 }else{
|
|
|
290 system("perl $path\/rfam.pl -i $data2 -ref $options{rfam} -v $mis_rfam -p $t -o $dir");
|
|
|
291 }
|
|
|
292 my $tag=join "\\;" ,@mark;
|
|
|
293 my $rfam_count=`perl $path\/count_rfam_express.pl -i $dir\/rfam_match\/rfam_mapped.bwt -tag $tag -o $dir\/rfam_match\/rfam_non-miRNA_annotation.txt`;
|
|
|
294 return 0;
|
|
|
295 }
|
|
|
296 sub genome{
|
|
|
297 if(defined $options{'idx'}){
|
|
|
298 system("perl $path\/matching.pl -i $data3 -g $genome_fa -v $mis -p $t -r $hit -o $dir -index $options{idx}") ;
|
|
|
299 }else{
|
|
|
300 system("perl $path\/matching.pl -i $data3 -g $genome_fa -v $mis -p $t -r $hit -o $dir ") ;
|
|
|
301 }
|
|
|
302 #=================== mapping sta ===================================================
|
|
|
303 my $map_file=$dir."genome_match\/genome_mapped\.fa";
|
|
|
304 open (MAP,"<$map_file")||die"$!";
|
|
|
305 print "\n#each sample mapping reads sta:\n\n";
|
|
|
306 print "#$mark\ttotal\n";
|
|
|
307 while (my $ID=<MAP>) {
|
|
|
308 chomp $ID;
|
|
|
309 my @tmp=split/\:/,$ID;
|
|
|
310 my @exp=split/\_/,$tmp[1];
|
|
|
311 $exp[-1] =~ s/^x//;
|
|
|
312 for (my $i=0;$i<@exp ;$i++) {
|
|
|
313 $map_read[$i]+=$exp[$i];
|
|
|
314 }
|
|
|
315 $map_tag++;
|
|
|
316 my $seq=<MAP>;
|
|
|
317 }
|
|
|
318 my $map_read=join"\t",@map_read;
|
|
|
319 print "$map_read\n\n";
|
|
|
320 print "#total mapped tags:$map_read\n\n";
|
|
|
321 close MAP;
|
|
|
322 return 0;
|
|
|
323 }
|
|
|
324
|
|
|
325 sub bwt2bed{
|
|
|
326 $cluster_file=$dir."cluster\/";
|
|
|
327 mkdir ("$cluster_file");
|
|
|
328 print "sam file changed to bed file\n";
|
|
|
329 my ($file) = $dir."genome_match\/genome_mapped\.bwt";
|
|
|
330
|
|
|
331 my $sam2bed=`perl $path\/sam2Bed_bowtie.pl -i $file -mark $sample_mark -o $bed `;
|
|
|
332 print "perl $path\/sam2Bed_bowtie.pl -i $file -mark $sample_mark -o $bed\n\n";
|
|
|
333 return 0;
|
|
|
334 }
|
|
|
335
|
|
|
336 sub cluster{
|
|
|
337 print "tags is ready to merged clusters\n\n";
|
|
|
338 my ($file) =$bed;
|
|
|
339 if ($cluster_mothod eq "F") {
|
|
|
340 my $cluster=`perl $path\/conventional.pl -i $file -d $distance_of_merged_tag -n $sample_number -mark $sample_mark -o $read -t $read_txt`;
|
|
|
341 print "Using converntional method\n perl $path\/conventional.pl -i $file -d $distance_of_merged_tag -n $sample_number -mark $sample_mark -o $read -t $read_txt\n\n";
|
|
|
342 }
|
|
|
343 elsif($cluster_mothod eq "T"){
|
|
|
344 my $cluster=`perl $path\/nibls.pl -f $file -m $distance_of_merged_tag -o $read -t $read_txt -mark $sample_mark`;
|
|
|
345 print "Using nibls method\n perl $path\/nibls.pl -f $file -m $distance_of_merged_tag -o $read -t $dir\/cluster.txt -mark $sample_mark\n\n";
|
|
|
346 }
|
|
|
347 else{print "\-p is wrong!\n\n";}
|
|
|
348 return 0;
|
|
|
349 }
|
|
|
350
|
|
|
351
|
|
|
352 sub quantify{
|
|
|
353 print "clusters is ready to quantified\n\n";
|
|
|
354 my @depth=@map_read;
|
|
|
355 pop @depth;
|
|
|
356 my $depth=join ",",@depth;
|
|
|
357 my $quantify=`perl $path\/quantify.pl -i $read -d $depth -o $rpkm`;
|
|
|
358 print "perl $path\/quantify.pl -i $read -d $depth -o $rpkm\n\n\n";
|
|
|
359 return 0;
|
|
|
360 }
|
|
|
361
|
|
|
362 sub phase{
|
|
|
363 $annotate_dir=$dir."annotate\/";
|
|
|
364 mkdir ("$annotate_dir");
|
|
|
365 print "clusters is to predict phase siRNA\n";
|
|
|
366 my $phase=`perl $path\/phased_siRNA.pl -i $read_txt -o $annotate_dir\/phase.out`;
|
|
|
367 print "perl $path\/phased_siRNA.pl -i $read_txt -o $annotate_dir\/phase.out\n\n\n";
|
|
|
368 return 0;
|
|
|
369 }
|
|
|
370
|
|
|
371 sub class{
|
|
|
372 print "clusters is ready to annotate by sources\n\n";
|
|
|
373 my $nat=$options{'nat'};
|
|
|
374 my $repeat=$options{'repeat'};
|
|
|
375 my $class=`perl $path\/ClassAnnotate.pl -i $rpkm -g $gff -n $nat -r $repeat -p $annotate_dir\/phase.out -o $annotate_dir\/sample_class.anno -t $annotate_dir\/nat.out -l $dir\/ref\/genelist.txt`;
|
|
|
376 print "perl $path\/ClassAnnotate.pl -i $rpkm -g $gff -n $nat -r $repeat -p $annotate_dir\/phase.out -o $annotate_dir\/sample_class.anno -t $annotate_dir\/nat.out -l $dir\/ref\/genelist.txt\n\n";
|
|
|
377 }
|
|
|
378
|
|
|
379 sub annotate{
|
|
|
380 print "clusters is ready to annotate by gff file\n\n";
|
|
|
381 my $file;
|
|
|
382 if (defined($options{'nat'})&&defined($options{'repeat'})) {
|
|
|
383 $file="$annotate_dir\/sample_class.anno";
|
|
|
384 }
|
|
|
385 else{
|
|
|
386 $file=$rpkm;
|
|
|
387 }
|
|
|
388 my $annotate=`perl $path\/Annotate.pl -i $file -g $dir\/ref\/genelist.txt -d $up_down_dis -o $annotate_dir\/sample_c_p.anno`;
|
|
|
389 print "perl $path\/Annotate.pl -i $file -g $dir\/ref\/genelist.txt -d $up_down_dis -o $annotate_dir\/sample_c_p.anno\n\n";
|
|
|
390 return 0;
|
|
|
391 }
|
|
|
392 sub get_genelist{
|
|
|
393
|
|
|
394 my $get_genelist=`perl $path\/get_genelist.pl -i $gff -o $dir\/ref\/genelist.txt`;
|
|
|
395 print "perl $path\/get_genelist.pl -i $gff -o $dir\/ref\/genelist.txt";
|
|
|
396 }
|
|
|
397
|
|
|
398 sub dec{
|
|
|
399 print "deg reading\n\n";
|
|
|
400 my $deg_file=$options{'deg'};
|
|
|
401 open IN,"<$deg_file";
|
|
|
402 my @deg;
|
|
|
403 my $s=0;
|
|
|
404 while (my $aline=<IN>) {
|
|
|
405 chomp $aline;
|
|
|
406 next if($aline=~/^\#/);
|
|
|
407 $deg[$s]=$aline;
|
|
|
408 my @ea=split/\s+/,$aline;
|
|
|
409 push @pairdir,"$ea[0]_VS_$ea[1]\/";
|
|
|
410 #print "$deg[$s]\n";
|
|
|
411 $s++;
|
|
|
412 }
|
|
|
413 close IN;
|
|
|
414 $deg_dir=$dir."deg\/";
|
|
|
415 mkdir ("$deg_dir");
|
|
|
416 my $max_process = 10;
|
|
|
417 my $pm = new Parallel::ForkManager( $max_process );
|
|
|
418 my $number=@deg-1;
|
|
|
419 foreach(0..$number){
|
|
|
420 $pm->start and next;
|
|
|
421 &dec_pel($deg[$_]);
|
|
|
422 $pm->finish;
|
|
|
423 }
|
|
|
424 $pm->wait_all_children;
|
|
|
425 }
|
|
|
426
|
|
|
427 sub dec_pel{
|
|
18
|
428 print "\n******************\nstart:\n";
|
|
11
|
429 Time();
|
|
|
430 my $sample=shift(@_);
|
|
|
431 my @each=split/\s+/,$sample;
|
|
|
432 print "$each[0]\t$each[1]\n";
|
|
|
433 my $deg_sample_dir=$deg_dir."$each[0]_VS_$each[1]\/";
|
|
|
434 mkdir ("$deg_sample_dir");
|
|
18
|
435 print "read: $read\n";
|
|
|
436 print "deg_sample_dir: $deg_sample_dir\n";
|
|
|
437 print "$id{$each[0]}\t$each[0]\n";
|
|
|
438 print "$id{$each[1]}\t$each[1]\n";
|
|
11
|
439 my $deg=`perl $path\/DEGseq_2.pl -i $read -outdir $deg_sample_dir -column1 $id{$each[0]} -mark1 $each[0] -column2 $id{$each[1]} -mark2 $each[1]`; #-depth1 -depth2
|
|
|
440 my $time2=time();
|
|
18
|
441 print "end:\n*************************\n";
|
|
11
|
442 Time();
|
|
|
443 sleep 1;
|
|
|
444 }
|
|
|
445
|
|
|
446 sub infor_merge{
|
|
|
447 my ($input,$mark);
|
|
|
448 foreach (@pairdir) {
|
|
|
449 print "@pairdir\n";
|
|
|
450 $mark.=" -mark $_ ";
|
|
|
451 $input.=" -i $dir/deg\/$_\/output_score\.txt ";
|
|
|
452 print "$input\n$mark\n";
|
|
|
453 }
|
|
|
454 my $infor_merge=`perl $path\/SampleDEGseqMerge.pl $input $mark -f $annotate_dir\/sample_c_p.anno -n $sample_number -o $dir\/total.result `;
|
|
|
455 print "perl $path\/SampleDEGseqMerge.pl $input $mark -f $annotate_dir\/sample_c_p.anno -n $sample_number -o $dir\/total.result\n\n";
|
|
|
456 }
|
|
|
457
|
|
|
458 sub infor_merge_no_dec{
|
|
|
459 my $infor_merge_no_dec=`cp $annotate_dir\/sample_c_p.anno $dir\/total.result`;
|
|
|
460 }
|
|
|
461
|
|
|
462 sub genome_length{
|
|
|
463 my $length=`perl $path\/count_ref_length.pl -i $genome_fa -o $dir\/ref\/genome\.length`;
|
|
|
464 print "perl $path\/count_ref_length.pl -i $genome_fa -o $dir\/ref\/genome\.length\n\n"
|
|
|
465
|
|
|
466 }
|
|
|
467
|
|
|
468 sub plot{
|
|
|
469 my $plot_file="$dir\/plot\/";
|
|
|
470 mkdir ("$plot_file");
|
|
|
471 my $genome_plot="$dir\/plot\/genome\/";
|
|
|
472 mkdir ("$genome_plot");
|
|
|
473 #genome cluster
|
|
|
474 my $span=defined($options{span})?$options{span}:50000;
|
|
|
475 foreach (1..$sample_number) {
|
|
|
476 my $mark=$mark[$_-1];
|
|
|
477 my $cen="";
|
|
|
478 if (defined $options{cen}) {
|
|
|
479 $cen="-cen $options{cen}";
|
|
|
480 }
|
|
|
481 my $plot=`perl $path\/sRNA_rpkm_distribution_along_genome.pl -c $rpkm -n $_ -mark $mark -span $span -l $dir\/ref\/genome\.length $cen -o $genome_plot\/$mark\.html -out $genome_plot\/$mark\.txt`;
|
|
|
482 print "perl $path\/sRNA_rpkm_distribution_along_genome.pl -c $rpkm -n $_ -mark $mark -span $span -l $dir\/ref\/genome\.length $cen -o $genome_plot\/$mark\.html -out $genome_plot\/$mark\.txt\n\n";
|
|
|
483 }
|
|
|
484
|
|
|
485 my $chr_plot_dir="$dir\/plot\/chr\/";
|
|
|
486 mkdir("$chr_plot_dir");
|
|
|
487 my %chr;
|
|
|
488 open LEN,"<$dir\/ref\/genome\.length";
|
|
|
489 while (my $aline=<LEN>) {
|
|
|
490 next if($aline=~/^\#/);
|
|
|
491 chomp $aline;
|
|
|
492 my @temp=split/\t/,$aline;
|
|
|
493 $chr{$temp[0]}=$temp[1];
|
|
|
494 }
|
|
|
495 close LEN;
|
|
|
496 foreach my $chr (sort keys %chr) {
|
|
|
497 my $cen="";
|
|
|
498 if (defined $options{cen}) {
|
|
|
499 $cen="-cen $options{cen}";
|
|
|
500 }
|
|
|
501 my $chr_plot=`perl $path\/chr_plot.pl -l $chr{$chr} -chro $chr -g $dir\/ref\/genelist.txt -span $span -c $rpkm -mark $sample_mark -o $chr_plot_dir\/$chr\.html`;
|
|
|
502 print "perl $path\/chr_plot.pl -l $chr{$chr} -chro $chr -g $dir\/ref\/genelist.txt -span $span -c $rpkm -mark $sample_mark -o $chr_plot_dir\/$chr\.html\n";
|
|
|
503 }
|
|
|
504 }
|
|
|
505
|
|
|
506 sub html{
|
|
|
507 my $pathfile="$dir/path.txt";
|
|
|
508 open PA,">$pathfile";
|
|
|
509 print PA "$config\n";
|
|
|
510 print PA "$preprocess\n";
|
|
|
511 print PA "$dir"."rfam_match\n";
|
|
|
512 print PA "$dir"."genome_match\n";
|
|
|
513 print PA "$cluster_file\n";
|
|
|
514 print PA "$annotate_dir\n";
|
|
|
515 if (defined($deg_dir)) {
|
|
|
516 print PA "$deg_dir\n";
|
|
|
517 }
|
|
|
518 close PA;
|
|
|
519 my $html=`perl $path\/html.pl -i $pathfile -format $format -o $dir/result.html`;
|
|
|
520 }
|
|
|
521
|
|
|
522 sub Time{
|
|
|
523 my $time=time();
|
|
|
524 my ($sec,$min,$hour,$day,$month,$year) = (localtime($time))[0,1,2,3,4,5,6];
|
|
|
525 $month++;
|
|
|
526 $year+=1900;
|
|
|
527 if (length($sec) == 1) {$sec = "0"."$sec";}
|
|
|
528 if (length($min) == 1) {$min = "0"."$min";}
|
|
|
529 if (length($hour) == 1) {$hour = "0"."$hour";}
|
|
|
530 if (length($day) == 1) {$day = "0"."$day";}
|
|
|
531 if (length($month) == 1) {$month = "0"."$month";}
|
|
|
532 print "$year-$month-$day $hour:$min:$sec\n";
|
|
|
533 return("$year-$month-$day-$hour-$min-$sec");
|
|
|
534 }
|
|
|
535 #################################################################################
|
