Mercurial > repos > chawhwa > neuma

comparison NEUMA-1.2.1/auto_NEUMA_PE.pl @ 0:c44c43d185ef draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
NEUMA-1.2.1 Uploaded
author chawhwa
date Thu, 08 Aug 2013 00:46:13 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:c44c43d185ef
1 #!/usr/bin/perl
2
3 if(@ARGV<8) { print "usage: $0 [options] -L=<read_length> -D=<maxdist> -1=<input_file1(mate1)> -2=<input_file2(mate2)> -U=<Utable_prefix(fullpath, before .gU.table or .iU.table)> --g2m=<gene2NM_file> --g2s=<gene2symbol_file> -b=<bowtie_dir(eg.bin/bowtie-0.12.7)> --bi=<bowtieindex> -o=<outputdir> -s=<sample_name>\n\nOrder of arguments can be interchangeable.\n\n"; exit; }
4
5
6 ## options ##
7 ($fastoption)= grep /^-f=/, @ARGV;
8 if(!defined($fastoption)){$fastoption='q'; }
9 else { $fastoption =~s/^-f=//; if($fastoption !~/^[fq]$/) { die "ERROR: wrong file type (-f).\n"; } }
10
11 ($coding_option)= grep /^-c=/, @ARGV;
12 if(!defined($coding_option)){$coding_option='n'; }
13 else { $coding_option =~s/^-c=//; if($coding_option !~/^[nc]$/) { die "ERROR: wrong coding option (-c).\n"; } }
14
15 ($datatype)= grep /^-d=/, @ARGV;
16 if(!defined($datatype)){$datatype='R'; }
17 else { $datatype =~s/^-d=//; if($datatype !~/^[RE]$/) { die "ERROR: wrong data type (-d).\n"; } }
18
19 ($parallel)= grep /^-p=/, @ARGV;
20 if(!defined($parallel)){$parallel=1; }
21 else { $parallel =~s/^-p=//; if($parallel !~/^[\d]+$/) { die "ERROR: wrong parallel (multi-thread) option (-p).\n"; } }
22
23 ($EUMAcut)= grep /^-t=/, @ARGV;
24 if(!defined($EUMAcut)){$EUMAcut=50; }
25 else { $EUMAcut =~s/^-t=//; if($EUMAcut !~/^[\d\.]+$/) { die "ERROR: wrong EUMAcut (-t).\n"; } }
26
27 ($StrandSpecificity)= grep /^--str=/, @ARGV;
28 if(!defined($StrandSpecificity)){ $StrandSpecificity='N'; }
29 else { $StrandSpecificity =~s/^--str=//; if($StrandSpecificity !~/^[SN]+$/) { die "ERROR: wrong strand specificity (--str).\n"; } }
30
31 ($skip_init)= grep /^--skip_init/, @ARGV;
32 if(!defined($skip_init)){ $skip_init=0; }
33 else { $skip_init=1; }
34
35 ($only_init)= grep /^--only_init/, @ARGV;
36 if(!defined($only_init)){ $only_init=0; }
37 else { $only_init=1; }
38
39 if($only_init==1 && $skip_init==1) { die "ERROR: --only_init and --skip_init cannot be used simultaneously.\n"; }
40
41 ($noNIR)= grep /^--noNIR/, @ARGV;
42 if(!defined($noNIR)){ $noNIR=0; }
43 else { $noNIR=1; }
44
45 ($noNEUMA)= grep /^--noNEUMA/, @ARGV;
46 if(!defined($noNEUMA)){ $noNEUMA=0; }
47 else { $noNEUMA=1; }
48
49 if($noNIR==1 && $noNEUMA==0) { die "ERROR: --noNIR must be used together with --noNEUMA.\n"; }
50
51 ($get_gReadcount)= grep /^--gReadcount/, @ARGV;
52 if(!defined($get_gReadcount)){ $get_gReadcount=0; }
53 else { $get_gReadcount=1; }
54
55 ($mm)= grep /^--mm=/, @ARGV;
56 if(!defined($mm)){ $mm=0; }
57 else { $mm =~s/^--mm=//; if($mm !~/^[\d\.]+$/) { die "ERROR: wrong number of mismatches (--mm).\n"; } }
58
59
60
61 ## required arguments ##
62
63 ($READ_LENGTH)= grep /^-L=/, @ARGV;
64 if(!defined($READ_LENGTH)){ die "ERROR: READ_LENGTH must be specified (-L).\n"; }
65 else { $READ_LENGTH =~s/^-L=//; if($READ_LENGTH !~/^[\d]+$/) { die "ERROR: wrong READ_LENGTH (-L).\n"; } }
66
67 ($MAXDIST)= grep /^-D=/, @ARGV;
68 if(!defined($MAXDIST)){ die "ERROR: MAXDIST must be specified (-D).\n"; }
69 else { $MAXDIST =~s/^-D=//; if($MAXDIST !~/^[\d]+$/) { die "ERROR: wrong MAXDIST (-D).\n"; } }
70
71 ($input_file1)= grep /^-1=/, @ARGV;
72 if(!defined($input_file1)){ die "ERROR: input_file1 must be specified (-1).\n"; }
73 else { $input_file1 =~s/^-1=//; if(!-e $input_file1) { die "ERROR: wrong input_file1 (-1).\n"; } }
74
75 ($input_file2)= grep /^-2=/, @ARGV;
76 if(!defined($input_file2)){ die "ERROR: input_file2 must be specified (-2).\n"; }
77 else { $input_file2 =~s/^-2=//; if(!-e $input_file2) { die "ERROR: wrong input_file2 (-2).\n"; } }
78
79 ($Utable_prefix)= grep /^-U=/, @ARGV;
80 if($only_init==0 && !defined($Utable_prefix)){ die "ERROR: Utable_prefix must be specified (-U).\n"; }
81 else { $Utable_prefix =~s/^-U=//; if(!-e $Utable_prefix.".gU.table" || !-e $Utable_prefix.".iU.table") { die "ERROR: wrong Utable_prefix (-U).\n"; } }
82
83 ($gene2NM_file)= grep /^--g2m=/, @ARGV;
84 if($only_init==0 && !defined($gene2NM_file)){ die "ERROR: gene2NM_file must be specified (--g2m).\n"; }
85 else { $gene2NM_file =~s/^--g2m=//; if(!-e $gene2NM_file) { die "ERROR: wrong gene2NM_file (--g2m).\n"; } }
86
87 ($gene2symbol_file)= grep /^--g2s=/, @ARGV;
88 if($only_init==0 && !defined($gene2symbol_file)){ die "ERROR: gene2symbol_file must be specified (--g2s).\n"; }
89 else { $gene2symbol_file =~s/^--g2s=//; if(!-e $gene2symbol_file) { die "ERROR: wrong gene2symbol_file (--g2s).\n"; } }
90
91 ($bowtie_dir)= grep /^-b=/, @ARGV;
92 if(!defined($bowtie_dir)){ die "ERROR: bowtie_dir must be specified (-b).\n"; }
93 else { $bowtie_dir =~s/^-b=//; if(!-d $bowtie_dir || !-e "$bowtie_dir/bowtie") { die "ERROR: wrong bowtie_dir(please avoid using '~') (-b).\n"; } }
94
95 ($bowtieindex)= grep /^--bi=/, @ARGV;
96 if(!defined($bowtieindex)){ die "ERROR: bowtieindex must be specified (--bi).\n"; }
97 else { $bowtieindex =~s/^--bi=//; if(!-e "$bowtieindex.1.ebwt") { die "ERROR: wrong bowtieindex (--bi).\n"; } }
98
99 ($basedir)= grep /^-o=/, @ARGV;
100 if(!defined($basedir)){ die "ERROR: output dir must be specified (-o).\n"; }
101 else { $basedir =~s/^-o=//; if($basedir=~/~/) { die "ERROR: wrong base_dir(please avoid using '~') (-o).\n"; } }
102
103 ($sample)= grep /^-s=/, @ARGV;
104 if(!defined($sample)){ die "ERROR: sample name must be specified (-s).\n"; }
105 else { $sample =~s/^-s=//; }
106
107
108 ##############
109
110
111 $gUtable_file = $Utable_prefix.".gU.table";
112 $iUtable_file = $Utable_prefix.".iU.table";
113
114 my $coding='';
115 if($coding_option eq 'n') {$coding='';}
116 elsif($coding_option eq 'c') {$coding='-C';}
117 else { die "ERROR: Error: invalid coding option.\n"; }
118
119 if($datatype eq 'R') { $mapping_stat_column = 4; }
120 elsif($datatype eq 'E') { $mapping_stat_column = 6; }
121 else { die "ERROR: Error: wrong datatype.\n"; }
122
123 if($StrandSpecificity eq 'S') { $bowtie_strand_option = "--norc"; }
124 else { $bowtie_strand_option = ""; }
125
126 ($scriptdir) = $0 =~ /(.+)\/[^\/]+$/;
127
128 $bowtieoutdir = "$basedir/bowtieout";
129 $readcount_dir = "$basedir/readcount";
130 $lendis_dir = "$basedir/insertlendis";
131 $EUMAdir = "$basedir/EUMA";
132 $FVKM_dir = "$basedir/FVKM";
133 $LVKM_dir = "$basedir/LVKM";
134
135 $mapping_stat_file = "$basedir/mapping_stat.$sample";
136
137 ($midrefix) = $bowtieindex=~/([^\/]+)$/;
138
139
140 if(!-d $bowtieoutdir) { `mkdir -p $bowtieoutdir`; }
141 if(!-d $readcount_dir) { `mkdir -p $readcount_dir`; }
142 if(!-d $lendis_dir) { `mkdir -p $lendis_dir`; }
143 if(!-d $EUMAdir) { `mkdir -p $EUMAdir`; }
144 if(!-d $FVKM_dir) { `mkdir -p $FVKM_dir`; }
145 if(!-d $LVKM_dir) { `mkdir -p $LVKM_dir`; }
146
147
148 #`dos2unix $gene2NM_file`;
149 #`dos2unix $gene2symbol_file`;
150
151
152
153 #$perlcommand = "perl -I $scriptdir";
154 $perlcommand = "perl";
155 $bestbowtieout_script = "$scriptdir/filter.best.from.bowtieout.3.pl";
156 $bowtie2insertlendis_script = "$scriptdir/bowtieout2insertlendis.PE.pl";
157 $bowtieout2mappingstat_script = "$scriptdir/bowtieout2mappingstat.3.pl";
158 $bowtieout2readcount_script = "$scriptdir/bowtie2genecount.11.pl";
159
160 $calc_gEUMA_script = "$scriptdir/calculate_gEUMA.2.pl";
161 $calc_iEUMA_script = "$scriptdir/calculate_iEUMA.2.pl";
162 $NIR2LVKM_script = "$scriptdir/NIR2LVKM.pl";
163
164
165 system("date");
166
167
168 if($skip_init==0){
169 print STDERR "Mapping reads using bowtie...\n";
170 system("$bowtie_dir/bowtie -$fastoption $coding --minins 0 --maxins $MAXDIST -v $mm -a --suppress 5,6,7 -p $parallel $bowtie_strand_option $bowtieindex -1 $input_file1 -2 $input_file2 > $bowtieoutdir/$sample.$midrefix.maxins$MAXDIST.mm$mm.bowtieout");
171
172 print STDERR "Filtering best-matching alignments...\n";
173 if($mm > 0) {
174 system("$perlcommand $bestbowtieout_script -d $datatype --rm $bowtieoutdir/$sample.$midrefix.maxins$MAXDIST.mm$mm.bowtieout 1");
175 }
176 else {
177 system("mv $bowtieoutdir/$sample.$midrefix.maxins$MAXDIST.mm$mm.bowtieout $bowtieoutdir/$sample.$midrefix.maxins$MAXDIST.mm$mm.best.bowtieout");
178 }
179
180 print STDERR "Computing fragment length distribution...\n";
181 system("$perlcommand $bowtie2insertlendis_script $bowtieoutdir/$sample.$midrefix.maxins$MAXDIST.mm$mm.best.bowtieout $READ_LENGTH > $lendis_dir/$sample.$midrefix.maxins$MAXDIST.mm0.i.insertlendis");
182 }
183
184 if($only_init==0) {
185
186 print STDERR "Computing the mapping stat...\n";
187 system("$perlcommand $bowtieout2mappingstat_script -d $datatype -m $MAXDIST -l $READ_LENGTH $bowtieoutdir/$sample.$midrefix.maxins$MAXDIST.mm$mm.best.bowtieout $sample 1 > $mapping_stat_file");
188 print STDERR "Mapping stat :\n";
189 system("cat $mapping_stat_file");
190
191 if($noNIR==0){
192 print STDERR "Computing NIRs...\n";
193 system("$perlcommand $bowtieout2readcount_script -d $datatype -m $MAXDIST -l $READ_LENGTH --gNIR -f -s $sample -v $gene2NM_file $bowtieoutdir/$sample.$midrefix.maxins$MAXDIST.mm$mm.best.bowtieout 1 > $readcount_dir/$sample.$midrefix.maxins$MAXDIST.mm$mm.gNIR");
194 system("$perlcommand $bowtieout2readcount_script -d $datatype -m $MAXDIST -l $READ_LENGTH --iNIR -f -s $sample -v $gene2NM_file $bowtieoutdir/$sample.$midrefix.maxins$MAXDIST.mm$mm.best.bowtieout 1 > $readcount_dir/$sample.$midrefix.maxins$MAXDIST.mm$mm.iNIR");
195 }
196
197 if($get_gReadcount==1) {
198 print STDERR "Computing total gene read counts (gReadcount)...\n";
199 print STDERR "$perlcommand $bowtieout2readcount_script -d $datatype -m $MAXDIST -l $READ_LENGTH -f -s $sample -v $gene2NM_file $bowtieoutdir/$sample.$midrefix.maxins$MAXDIST.mm$mm.best.bowtieout 1 > $readcount_dir/$sample.$midrefix.maxins$MAXDIST.mm$mm.gReadcount\n";
200 system("$perlcommand $bowtieout2readcount_script -d $datatype -m $MAXDIST -l $READ_LENGTH -f -s $sample -v $gene2NM_file $bowtieoutdir/$sample.$midrefix.maxins$MAXDIST.mm$mm.best.bowtieout 1 > $readcount_dir/$sample.$midrefix.maxins$MAXDIST.mm$mm.gReadcount");
201 }
202
203 if($noNEUMA==0){
204 print STDERR "Computing gEUMA and iEUMA...\n";
205 system("$perlcommand $calc_gEUMA_script $gUtable_file $lendis_dir/$sample.$midrefix.maxins$MAXDIST.mm0.i.insertlendis $READ_LENGTH > $EUMAdir/$sample.$midrefix.maxins$MAXDIST.mm$mm.gEUMA");
206 system("$perlcommand $calc_iEUMA_script $iUtable_file $lendis_dir/$sample.$midrefix.maxins$MAXDIST.mm0.i.insertlendis $READ_LENGTH > $EUMAdir/$sample.$midrefix.maxins$MAXDIST.mm$mm.iEUMA");
207
208 print STDERR "Computing FVKM and LVKM...\n";
209 system("$perlcommand $NIR2LVKM_script $sample $readcount_dir $midrefix.maxins$MAXDIST.mm$mm $EUMAdir $FVKM_dir $LVKM_dir $EUMAcut $mapping_stat_column $mapping_stat_file 1 2 $gene2NM_file $gene2symbol_file $datatype $scriptdir");
210
211 }
212 }
213
214 print STDERR "Done.\n";
215
216 system("date");
217
218