add option -s to specify strand information

xulab-nwafu · Sep 3, 2017 · e57a4b9 · e57a4b9
1 parent 17c96aa
commit e57a4b9
Show file tree

Hide file tree

Showing 4 changed files with 32 additions and 9 deletions.
diff --git a/README.md b/README.md
@@ -7,7 +7,7 @@ N.B.: the scripts `dnaseq` and `rnaseq` are work with the hardware and software
 HIgh-Throughput Sequencing (hits) pipeline
 ------------------------------------------
 ```
-hits 0.3.0
+hits 0.4.0
 
 Usage:
 
@@ -17,6 +17,9 @@ Usage:
     -d  directory that containing FASTQ file(s), .fq/.fastq(.gz)
     -g  reference genome
     -a  annotation file in GTF format [work with -p R]
+    -s  strand-specific information [work with -p R]
+        1 - stranded [default]
+        0 - unstranded
     -o  output directory
     -t  number of alignment threads
     -h  help
@@ -25,7 +28,7 @@ Usage:
 Examples:
 
     hits -p D -d fastq_dir -g genome.fa -o output_dir -t 8
-    hits -p R -d fastq_dir -g genome.fa -a gene.gtf -o output_dir -t 8
+    hits -p R -d fastq_dir -g genome.fa -a gene.gtf -s 1 -o output_dir -t 8
 
 Reporting Bugs:
 

diff --git a/hits b/hits
@@ -10,7 +10,7 @@ use warnings;
 use Getopt::Std;
 
 my $appname = 'hits';
-my $version = '0.3.0';
+my $version = '0.4.0';
 my %options = ();
 
 init();
@@ -107,9 +107,13 @@ sub rnaseq {
 	}
 
 	print STDERR "Runing RNAseq pipeline\n";
+	my $rnaseq_strand_option = '';
+	if ($options{'s'} == 1) {
+		$rnaseq_strand_option = '--rna-strandness ';
+	}
 	foreach my $fq (sort keys %{$fqs}) {
 		print STDERR "Mapping $fq ...\n";
-		my $hisat2_options = "--rna-strandness RF --min-intronlen 20 --max-intronlen 4000 --no-unal --no-temp-splicesite " .
+		my $hisat2_options = $rnaseq_strand_option .  "--min-intronlen 20 --max-intronlen 4000 --no-unal --no-temp-splicesite " .
 		"--novel-splicesite-outfile $dir/bam/$fq.splicesite --novel-splicesite-infile $dir/bam/$fq.splicesite";
 		my $map_cmd = "hisat2 -p $options{'t'} $hisat2_options -x " . $options{'g'} .
 			' -1 ' . $fqs->{$fq}->{'1'} .
@@ -231,7 +235,7 @@ sub init {
 	}
 
 	# Options
-	getopts("p:d:g:a:o:t:hv", \%options) or die "$!\n" . usage();
+	getopts("p:d:g:a:s:o:t:hv", \%options) or die "$!\n" . usage();
 
 	# Check help
 	if ($options{'h'}) {
@@ -274,14 +278,27 @@ sub init {
 		}
 	}
 
-	# In RNAseq, must specify option -a
+	# Check important RNA-seq options
 	if (uc($options{'p'}) eq 'R') {
+		# In RNAseq, must specify option -a
 		unless ( defined($options{'a'}) ) {
 			print STDERR "Aborted: No annotation file specified!\n";
 			usage();
 		}
+		# Defalut: stranded RNA-seq
+		unless (defined($options{'s'})) {
+			print STDERR "No strand information of RNA-seq specified, set to defualt (1-stranded)\n";
+			$options{'s'} = 1;
+		} else {
+			unless ($options{'s'} == 1 or $options{'s'} == 0) {
+				print STDERR "Aborted: Unknown strand information of RNA-seq specified!\n";
+				print STDERR "required: 1 or 0, specified as $options{'a'}\n";
+				usage();
+			}
+		}
 	}
 
+
 	return;
 
 }
@@ -317,6 +334,9 @@ Usage:
     -d  directory that containing FASTQ file(s), .fq/.fastq(.gz)
     -g  reference genome
     -a  annotation file in GTF format [work with -p R]
+    -s  strand-specific information [work with -p R]
+        1 - stranded [default]
+        0 - unstranded
     -o  output directory
     -t  number of alignment threads
     -h  help
@@ -325,7 +345,7 @@ Usage:
 Examples:
 
     $appname -p D -d fastq_dir -g genome.fa -o output_dir -t 8
-    $appname -p R -d fastq_dir -g genome.fa -a gene.gtf -o output_dir -t 8
+    $appname -p R -d fastq_dir -g genome.fa -a gene.gtf -s 1 -o output_dir -t 8
 
 Reporting Bugs:
 

diff --git a/rnaseq b/rnaseq
@@ -26,5 +26,5 @@ fi
 dir_name=`dirname $1`
 base_name=`basename $1`
 
-nohup hits -p R -g ~/db/Fusarium_graminearum/current/hisat2/FG.RR.*.genome.fa -a ~/db/Fusarium_graminearum/current/FG.RR.*.gene.gtf -d $1 -o $dir_name/$base_name.outdir -t 8 1>$dir_name/$base_name.rnaseq_log.txt 2>&1 &
+nohup hits -p R -g ~/db/Fusarium_graminearum/current/hisat2/FG.RR.*.genome.fa -a ~/db/Fusarium_graminearum/current/FG.RR.*.gene.gtf -s 1 -d $1 -o $dir_name/$base_name.outdir -t 8 1>$dir_name/$base_name.rnaseq_log.txt 2>&1 &
 echo "Job submitted, outputs will be in $dir_name/$base_name.outdir"
diff --git a/rnaseq.ss b/rnaseq.ss
@@ -26,5 +26,5 @@ fi
 dir_name=`dirname $1`
 base_name=`basename $1`
 
-nohup hits -p R -g ~/db/Sclerotinia_sclerotiorum/current/hisat2/Sclerotinia_sclerotiorum.genome.fa -a ~/db/Sclerotinia_sclerotiorum/current/Sclerotinia_sclerotiorum.gene.gtf -d $1 -o $dir_name/$base_name.outdir -t 8 1>$dir_name/$base_name.rnaseq_log.txt 2>&1 &
+nohup hits -p R -g ~/db/Sclerotinia_sclerotiorum/current/hisat2/Sclerotinia_sclerotiorum.genome.fa -a ~/db/Sclerotinia_sclerotiorum/current/Sclerotinia_sclerotiorum.gene.gtf -s 1 -d $1 -o $dir_name/$base_name.outdir -t 8 1>$dir_name/$base_name.rnaseq_log.txt 2>&1 &
 echo "Job submitted, outputs will be in $dir_name/$base_name.outdir"