1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
|
my($gtf,$submit,$hold_jid,$fasta_ref,$jobid,$prefix,@prefix_list);
#options on the command line
GetOptions( "gtf:s" => \$gtf,
"submit:s" => \$submit,
"fasta_ref:s" => \$fasta_ref);
#management options
$submit ||= "no";
if (!defined $gtf) {
die "Error : missing gtf file\n";
}
if (!defined $fasta_ref) {
die "Error : missing fasta_ref file\n";
}
my (@gtf_files);
if (!defined $gtf) {
die "Error : define input gtf file(s) with -gtf file1.gtf[,file2.gtf]\n";
} else {
@gtf_files = split(",",$gtf);
}
foreach my $gtf (@gtf_files) {
#create a directory for each sample
my @nn=split(/\_/,$gtf);
$prefix=$nn[0];
my $p=$nn[3];
my @nn2=split(/\./,$p);
$prefix=$nn2[0];
my $command = "mkdir $prefix";
qx{$command};
#generate qsub file and launch it (extraction of transcripts)
my @jobs_id1;
my $qsuboutfile1 = "step1_$prefix\_extraction_transcripts.qsub";
open (QSUB,">$prefix/$qsuboutfile1")or die "\nError : $qsuboutfile1 file could not be created: $!\n";
print QSUB qq(#!/bin/bash
#\$ -S /bin/bash
#\$ -o $qsuboutfile1.out
#\$ -e $qsuboutfile1.err
#\$ -cwd
#\$ -q workq
# command(s):
hostname >&2
~/bin/perl/annotation/cuff.pl -gtf $gtf > $prefix\_new_transcripts.gtf
gffread -w $prefix\_new_transcrits.fa -g $fasta_ref $prefix\_new_transcripts.gtf
mkdir fasta_eclate
SplitMultifasta.pl -i $prefix\_new_transcrits.fa -n 50 -o ./fasta_eclate
cd fasta_eclate
for i in *; do echo "blastall -p blastx -a 10 -d nr -i \$i -e 0.00001 -o `basename \$i .fa`_nr.blast"; done > ../blast_nr.sh
for i in *; do echo "blastall -p blastn -a 10 -d /bank/blastdb/rfam -i \$i -e 0.00001 -o `basename \$i .fa`_rfam.blast"; done > ../blast_rfam.sh
for i in *; do echo "blastall -p blastx -a 10 -d /bank/blastdb/uniprot_sprot -i \$i -e 0.00001 -o `basename \$i .fa`_uniprot.blast"; done > ../blast_uniprot.sh
cd ..
cat blast_nr.sh blast_rfam.sh blast_uniprot.sh > blast.sh
mv blast.sh ./fasta_eclate
);
close QSUB;
#launch qsub job if asked
if ($submit eq "yes") {
my $dependencies = "";
if (defined $hold_jid) {
$dependencies = "-hold_jid $hold_jid";
}
my $instruction = "cd $prefix/; qsub $dependencies $qsuboutfile1 | awk '{print \$3}' ";
$jobid=qx{$instruction};
chomp $jobid;
push (@jobs_id1,$jobid);
print STDERR ".";
} else {
push (@jobs_id1,"");
}
if ($submit eq "yes") {
print STDERR "\n";
}
if ($submit eq "yes") {
print STDERR ".";
}
# create R file to built histogram of descriptive data
my $stat = "stats_$prefix\.R";
open (ROUT,">$prefix/$prefix\_stats.R") or die "\nError : $stat file could not be created: $!\n";
print ROUT qq(##read input files
.........
);
close ROUT;
#generate qsub file and launch it (blastall on data base)
my @jobs_id2;
my $qsuboutfile2 = "step2_$prefix\_blastall.qsub";
open (QSUB,">$prefix/$qsuboutfile2")or die "\nError : $qsuboutfile2 file could not be created: $!\n";
print QSUB qq(#!/bin/bash
#\$ -S /bin/bash
#\$ -o $qsuboutfile2.out
#\$ -e $qsuboutfile2.err
#\$ -cwd
#\$ -q workq
# command(s):
#hostname >&2
cd fasta_eclate/
qarray -N blast_jobs -l mem=16G -l h_vmem=32G -pe parallel_smp 4 blast.sh
cd ..
perl ~/bin/perl/annotation/stat_1.pl -gtf_file $gtf -| perl ~/bin/perl/annotation/stat_2.pl > $prefix\_stats_1.txt
perl ~/bin/perl/annotation/stat_3.pl -fasta $prefix\_new_transcrits.fa > $prefix\_N_number.txt
);
close QSUB;
#launch qsub job if asked
if ($submit eq "yes") {
my $dependencies = "-hold_jid ".join(",",@jobs_id1);
my $instruction = "cd $prefix/;qsub $dependencies ./$qsuboutfile2 | awk '{print \$3};'";
$jobid=qx{$instruction};
chomp $jobid;
push (@jobs_id2,$jobid);
print STDERR ".";
} else {
push (@jobs_id2,"");
}
#generate qsub file and launch it (blastall on 3 databases)
my @jobs_id2_1;
my $qsuboutfile2_1 = "blast.sh";
#launch qsub job if asked
if ($submit eq "yes") {
my $dependencies = "-hold_jid ".join(",",@jobs_id1);
my $instruction = "cd $prefix/fasta_eclate/; qarray -N blast_jobs $dependencies -l mem=16G -l h_vmem=32G -pe parallel_smp 4 $qsuboutfile2_1 | awk '{print \$3}' ";
$jobid=qx{$instruction};
chomp $jobid;
push (@jobs_id2_1,$jobid);
print STDERR ".";
} else {
push (@jobs_id2_1,"");
}
if ($submit eq "yes") {
print STDERR "\n";
}
if ($submit eq "yes") {
print STDERR ".";
} |
