Mercurial > repos > sarahinraauzeville > star

comparison sm_STAR2_V2.pl @ 6:b7ed4ee33fb6 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Deleted selected files
author sarahinraauzeville
date Mon, 09 Dec 2019 10:55:33 -0500
parents c4fc8ff6e280
children
comparison
equal deleted inserted replaced
5:c4fc8ff6e280 6:b7ed4ee33fb6
1 #!/usr/bin/perl -w
2
3 # usage : perl sm_STAR.pl <read1.fastq.gz> <read2.fastq.gz>
4 # 10/02/2014 - Wrapper du traitement des données RNAseq
5 # Sarah Maman
6 # Copyright (C) 2014 INRA
7 # This program is free software: you can redistribute it and/or modify
8 # it under the terms of the GNU General Public License as published by
9 # the Free Software Foundation, either version 3 of the License, or
10 # (at your option) any later version.
11 #
12 # This program is distributed in the hope that it will be useful,
13 # but WITHOUT ANY WARRANTY; without even the implied warranty of
14 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 # GNU General Public License for more details.
16 #
17 # You should have received a copy of the GNU General Public License
18 # along with this program. If not, see <http://www.gnu.org/licenses/>.
19 #
20 use strict;
21 use File::Basename;
22 use Getopt::Long;
23 use lib "$ENV{'MY_GALAXY_DIR'}";
24 use GalaxyPath;
25
26 my $cfg = GalaxyPath->new( -file => $ENV{"GALAXY_CONFIG_FILE"});
27 my $PATH = $cfg->my_path( 'workPath', 'MYWORKSPACE' );
28 my $STAR = $cfg->my_path( 'toolsPath', 'STAR_PATH' );
29
30
31
32 my $Nthreads;
33 my $genome_path;
34 my $reads_selector;
35 my $input_read;
36 my $Read1fastqgz;
37 my $Read2fastqgz;
38 my $alignIntronMin;
39 my $alignIntronMax;
40 my $outFilterMismatchNmax;
41 my $orientation;
42 my $refownfastaref;
43 my $refselector;
44 my $refowngtf;
45 my $compress;
46 my $cufflinks;
47 my $outputfile;
48 my $outputfileT;
49 my $outputlogSJ;
50 my $outputlogfinal;
51
52
53 Getopt::Long::Configure( 'no_ignorecase', 'bundling' );
54 GetOptions (
55 'runThreadN=i' => \$Nthreads,
56 'genomeDir=s' => \$genome_path,
57 'refselector=s' => \$refselector,
58 'refownfastaref=s' => \$refownfastaref,
59 'refowngtf=s' => \$refowngtf,
60 'compress=s' => \$compress,
61 'cufflinks=s' => \$cufflinks,
62 'readsselector=s'=> \$reads_selector,
63 'readFilesIn1=s' => \$Read1fastqgz,
64 'readFilesIn2=s' => \$Read2fastqgz,
65 'readsinputread=s' => \$input_read,
66 'alignIntronMin=i' => \$alignIntronMin,
67 'alignIntronMax=i' => \$alignIntronMax,
68 'outFilterMismatchNmax=i' => \$outFilterMismatchNmax,
69 'orientation=s' => \$orientation,
70 'outputfile=s' => \$outputfile,
71 'outputfileT=s' => \$outputfileT,
72 'outputlogfinal=s' => \$outputlogfinal,
73 'outputlogSJ=s' => \$outputlogSJ
74 ) or die "Usage: Error in command line arguments\n";
75
76 my $cmd1 = ''; my $cmd2 ='';
77 my $cmd3 = ''; my $cmd4 ='';
78
79 #STAR --runThreadN 4 --runMode genomeGenerate --genomeDir /work/smaman/TP_RNAseq/INDEX/ --genomeFastaFiles ITAG2.3_genomic_Ch6.fasta --sjdbGTFfile ITAG_pre2.3_gene_models_Ch6.gtf --sjdbOverhang 100
80
81 #smaman@node001 /work/smaman/TP_RNAseq $ ls -ltrah INDEX
82 #-rw-r--r-- 1 smaman BIOINFO 331 17 juil. 11:55 genomeParameters.txt
83 #-rw-r--r-- 1 smaman BIOINFO 387K 17 juil. 11:55 exonGeTrInfo.tab
84 #-rw-r--r-- 1 smaman BIOINFO 53K 17 juil. 11:55 geneInfo.tab
85 #-rw-r--r-- 1 smaman BIOINFO 151K 17 juil. 11:55 transcriptInfo.tab
86 #-rw-r--r-- 1 smaman BIOINFO 171K 17 juil. 11:55 exonInfo.tab
87 #-rw-r--r-- 1 smaman BIOINFO 325K 17 juil. 11:55 sjdbList.fromGTF.out.tab
88 #-rw-r--r-- 1 smaman BIOINFO 272K 17 juil. 11:55 sjdbInfo.txt
89 #-rw-r--r-- 1 smaman BIOINFO 325K 17 juil. 11:55 sjdbList.out.tab
90 #-rw-r--r-- 1 smaman BIOINFO 11 17 juil. 11:55 chrName.txt
91 #-rw-r--r-- 1 smaman BIOINFO 9 17 juil. 11:55 chrLength.txt
92 #-rw-r--r-- 1 smaman BIOINFO 11 17 juil. 11:55 chrStart.txt
93 #-rw-r--r-- 1 smaman BIOINFO 20 17 juil. 11:55 chrNameLength.txt
94 #-rw-r--r-- 1 smaman BIOINFO 47M 17 juil. 11:55 Genome
95 #-rw-r--r-- 1 smaman BIOINFO 360M 17 juil. 11:55 SA
96 #-rw-r--r-- 1 smaman BIOINFO 1,5G 17 juil. 11:55 SAindex
97
98
99 #STAR --readFilesIn WTr1.fastq WTr2.fastq --genomeDir /work/smaman/TP_RNAseq/INDEX/ --sjdbGTFfile ITAG_pre2.3_gene_models_Ch6.gtf --outSAMtype BAM SortedByCoordinate --alignIntronMin 20 --alignIntronMax 1000000 --outFilterMismatchNmax 10 --outSAMtype BAM SortedByCoordinate --runThreadN 4 --outFileNamePrefix galaxyName --outSAMstrandField intronMotif --outFilterIntronMotifs RemoveNoncanonical --outFilterType BySJout --quantMode TranscriptomeSAM
100
101 #-rw-r--r-- 1 smaman BIOINFO 45M 26 mars 2015 ITAG2.3_genomic_Ch6.fasta
102 #-rw-r--r-- 1 smaman BIOINFO 1,6M 26 mars 2015 ITAG_pre2.3_gene_models_Ch6.gtf
103 #-rw-r--r-- 1 smaman BIOINFO 29 26 mars 2015 ITAG2.3_genomic_Ch6.fasta.fai
104 #-rw-r--r-- 1 smaman BIOINFO 614 17 juil. 10:20 WTr1.fastq
105 #-rw-r--r-- 1 smaman BIOINFO 589 17 juil. 10:20 WTr2.fastq
106 #-rw-r--r-- 1 smaman BIOINFO 14K 17 juil. 11:55 Log.out
107 #-rw-r--r-- 1 smaman BIOINFO 35K 17 juil. 12:03 galaxyNameAligned.toTranscriptome.out.bam
108 #-rw-r--r-- 1 smaman BIOINFO 637 17 juil. 12:03 galaxyNameAligned.sortedByCoord.out.bam +++++++++
109 #-rw-r--r-- 1 smaman BIOINFO 0 17 juil. 12:03 galaxyNameSJ.out.tab ++++++++++++++++
110 #-rw-r--r-- 1 smaman BIOINFO 246 17 juil. 12:03 galaxyNameLog.progress.out
111 #-rw-r--r-- 1 smaman BIOINFO 1,7K 17 juil. 12:03 galaxyNameLog.final.out +++++++++++++++
112 #-rw-r--r-- 1 smaman BIOINFO 16K 17 juil. 12:03 galaxyNameLog.out
113
114
115
116 #workspace
117 my $debug = 0; #Mode debug
118 if ($debug == 0)
119 {
120 print STDOUT "Debug mode OK \n";
121 }
122 else
123 {
124 $PATH = dirname($outputfile);
125 print STDOUT "No debug \n";
126 }
127
128
129 #Récuperer le numero (unique) de l'output afin, si besoin, de créer un répertoire de travail unique dans /work/galaxy-dev/workspace
130 my ($nb) = ($outputfile=~/dataset_(\d+)\.\S+$/);
131
132 #Repertoire de sortie cree par le script, verif des droits d'ecriture sur ce repertoire de sortie
133 `cd $PATH/; mkdir $nb/; chmod -R 777 $nb/; cd $nb/;`;
134 my $dirresults= "$PATH/".$nb;
135
136 print STDOUT "Job working directory : $dirresults \n";
137
138
139 if ($refselector eq "ownfasta"){
140 my $cmdSTARindex="(cd $dirresults/; mkdir INDEX/; chmod 777 INDEX/; $STAR --runThreadN $Nthreads --runMode genomeGenerate --genomeDir $dirresults/INDEX --genomeFastaFiles $refownfastaref --sjdbGTFfile $refowngtf --sjdbOverhang 100) >& ./out_Starindex.log 2>&1";
141 system $cmdSTARindex;
142 #Info pour les biologistes
143 print STDOUT "STAR Genome Generate : \n\n $cmdSTARindex \n\n ";
144 $genome_path = "$dirresults/INDEX/";
145 }
146
147 my $addcuff;
148 if ($cufflinks eq "cuff"){
149 $addcuff="--outSAMstrandField intronMotif --outFilterIntronMotifs RemoveNoncanonical --outFilterType BySJout --quantMode TranscriptomeSAM ";
150 }else{
151 $addcuff="";
152 }
153
154
155 my $cat;
156 if ($reads_selector eq "single"){
157
158 my $in;
159 if ($compress eq "compress"){
160 #Si besoin, recupération du fichier de configuration avec modification de l extension
161 `ln -s $input_read $dirresults/input_read.fastq.gz;`;
162 $in = "$dirresults/input_read.fastq.gz";
163 $cat="--readFilesCommand zcat";
164 }else
165 {`ln -s $input_read $dirresults/input_read.fastq;`;
166 $in = "$dirresults/input_read.fastq";
167 $cat="";}
168
169 if ($orientation eq "No"){
170 $cmd1 = "(cd $dirresults; $STAR --runThreadN $Nthreads --genomeDir $genome_path --readFilesIn $in --outSAMtype BAM SortedByCoordinate --alignIntronMin $alignIntronMin --alignIntronMax $alignIntronMax --outFilterMismatchNmax $outFilterMismatchNmax $cat --outFileNamePrefix $nb $addcuff) >& ./out_Star.log 2>&1";
171 system $cmd1;
172 #Info pour les biologistes
173 print STDOUT "STAR command run on cluster without oriented reads : \n\n $cmd1 \n\n ";
174 }
175 else
176 {
177 $cmd2 = "(cd $dirresults; $STAR --runThreadN $Nthreads --genomeDir $genome_path --readFilesIn $in --outSAMtype BAM SortedByCoordinate --alignIntronMin $alignIntronMin --alignIntronMax $alignIntronMax --outFilterMismatchNmax $outFilterMismatchNmax $cat --outFileNamePrefix $nb $addcuff) >& ./out_Star.log 2>&1";
178 system $cmd2;
179 #Info pour les biologistes
180 print STDOUT "STAR command run on cluster with oriented reads : \n\n $cmd2 \n\n
181 Instead, you need to run Cufflinks with the library option --library-type options. For example, cufflinks <…> -library-type fr-firststrand should be used for the “standard” dUTP protocol. This option has to be used only for Cufflinks runs and not for STAR runs.\n\n";
182 }
183 }else{
184
185
186 my $in1;
187 my $in2;
188 if ($compress eq "compress"){
189 #Si besoin, recupération du fichier de configuration avec modification de l extension
190 `ln -s $Read1fastqgz $dirresults/Read1.fastq.gz; ln -s $Read2fastqgz $dirresults/Read2.fastq.gz;`;
191 $in1="$dirresults/Read1.fastq.gz";
192 $in2="$dirresults/Read2.fastq.gz";
193 $cat="--readFilesCommand zcat";
194 }else
195 {`ln -s $Read1fastqgz $dirresults/Read1.fastq; ln -s $Read2fastqgz $dirresults/Read2.fastq;`;
196 $in1="$dirresults/Read1.fastq";
197 $in2="$dirresults/Read2.fastq";
198 $cat="";}
199
200
201 if ($orientation eq "No"){
202 $cmd3 = "(cd $dirresults; $STAR --runThreadN $Nthreads --genomeDir $genome_path --readFilesIn $in1 $in2 --outSAMtype BAM SortedByCoordinate --alignIntronMin $alignIntronMin --alignIntronMax $alignIntronMax --outFilterMismatchNmax $outFilterMismatchNmax $cat --outFileNamePrefix $nb $addcuff) >& ./out_Star.log 2>&1";
203 system $cmd3;
204 #Info pour les biologistes
205 print STDOUT "STAR command run on cluster without oriented reads : \n\n $cmd3 \n\n ";
206 }
207 else
208 {
209 $cmd4 = "(cd $dirresults; $STAR --runThreadN $Nthreads --genomeDir $genome_path --readFilesIn $in1 $in2 --outSAMtype BAM SortedByCoordinate --alignIntronMin $alignIntronMin --alignIntronMax $alignIntronMax --outFilterMismatchNmax $outFilterMismatchNmax $cat --outFileNamePrefix $nb $addcuff) >& ./out_Star.log 2>&1";
210 #Info pour les biologistes
211 system $cmd4;
212 print STDOUT "STAR command run on cluster with oriented reads : \n\n $cmd4 \n\n
213 Instead, you need to run Cufflinks with the library option --library-type options. For example, cufflinks <…> -library-type fr-firststrand should be used for the “standard” dUTP protocol. This option has to be used only for Cufflinks runs and not for STAR runs.\n\n";
214 }
215
216
217 }
218
219 #Si besoin :
220 #TEST 1 : command ligne on vm-galaxy
221 #TEST 2 perl Galaxy file : perl script.pl path/to/tests/files/used/for/galaxy/perl/script out1
222
223 #Recuperation des fichiers par Galaxy
224 #-rw-r--r-- 1 smaman BIOINFO 35K 17 juil. 12:03 galaxyNameAligned.toTranscriptome.out.bam +++++
225 #-rw-r--r-- 1 smaman BIOINFO 637 17 juil. 12:03 galaxyNameAligned.sortedByCoord.out.bam +++++++++
226 #-rw-r--r-- 1 smaman BIOINFO 0 17 juil. 12:03 galaxyNameSJ.out.tab ++++++++++++++++
227 #-rw-r--r-- 1 smaman BIOINFO 1,7K 17 juil. 12:03 galaxyNameLog.final.out +++++++++++++++
228 my $bam = glob("$dirresults/*$nb*Aligned.sortedByCoord.out.bam");
229 if (! -e $bam){print STDERR "Aligned.sortedByCoord.out.bam file not found. \n";}else{`cp -a $bam $outputfile`;}
230 my $bamT = glob("$dirresults/*$nb*Aligned.toTranscriptome.out.bam");
231 if (! -e $bamT){print STDERR "Aligned.toTranscriptome.out.bam file not found. \n";}else{`cp -a $bamT $outputfileT`;}
232 my $logSJ = glob("$dirresults/$nb*SJ.out.tab");
233 if (! -e $logSJ){print STDERR "SJ.out.tab log file not found. \n";}else{`cp -a $logSJ $outputlogSJ`;}
234 my $logfinal = glob("$dirresults/$nb*Log.final.out");
235 if (! -e $logfinal){print STDERR "Log.final.out log file not found. \n";}else{`cp -a $logfinal $outputlogfinal`;}
236