[Bio] / FigKernelScripts / rs_pipeline.pl Repository:
ViewVC logotype

View of /FigKernelScripts/rs_pipeline.pl

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (download) (as text) (annotate)
Mon Jun 3 16:44:35 2013 UTC (6 years, 5 months ago) by overbeek
Branch: MAIN
CVS Tags: rast_rel_2014_0729, rast_rel_2014_0912, HEAD
the rna_seq pipeline

use strict;
use Data::Dumper;
use File::Basename;
use Getopt::Long;
use SeedEnv;
use PG;

my $usage = "usage: rs_pipeline -d DataDir\n";
my $dataD;

my $rc  = GetOptions('d=s' => \$dataD);

if ((! $rc) || 
    (! -d $dataD))
{ 
    print STDERR $usage; exit ;
}

opendir(REF,"$dataD/ReferenceGenomes")  || die "Where are the reference genomes?";
my @refs = grep { $_ !~ /^\./ } readdir(REF);
closedir(REF);

opendir(SETS,"$dataD/Sets") || die "Where are the $dataD/Sets?";
my @sets = grep { $_ !~ /^\./ } readdir(SETS);
closedir(SETS);

foreach my $set (@sets)
{
    opendir(RUN,"$dataD/Sets/$set") || die "cannot open set $set";
    my @runs = grep { $_ !~ /^\./ } readdir(RUN);
    closedir(RUN);

    foreach my $run (@runs)
    {
	foreach my $ref (@refs)
	{
	    my $outputD = "$dataD/Output/$set.$run.$ref";
	    my $fastq   = "$dataD/Sets/$set/$run/fastq.gz";
	    (-s $fastq) || die "Where is $fastq?";
	    my $fasta   = "$dataD/Sets/$set/$run/fasta";
	    if (! -s $fasta ) 
	    {
		my $i;
		&SeedUtils::run("zcat $fastq | perl -e '$i=0;while(<>){if(/^\@/&&$i==0){s/^\@/\>/;print;}elsif($i==1){print;$i=-3}$i++;}' > $fasta");
	    }

	    if (! -s "$outputD/blast")
	    {
		die "Make blast file";
	    }

	    if (! -s "$outputD/filtered.blast")
	    {
		&SeedUtils::run("perl rs_get_rnas.pl -d $dataD -r $ref < $outputD/blast > $outputD/filtered.blast");
	    }
	    if (! -s "$outputD/filtered.blast.profile")
	    {
		&SeedUtils::run("perl rs_create_ref_profile.pl < $outputD/filtered.blast > $outputD/filtered.blast.profile");
		&SeedUtils::run("perl rs_create_ref_profile.pl < $outputD/blast > $outputD/blast.profile");
	    }
	    if (! -s "$outputD/peg.counts")
	    {
		$ref =~ /^(\d+\.\d+)/;
		my $g = $1;
		&SeedUtils::run("perl rs_get_peg_counts.pl -d $dataD -g $dataD/ReferenceGenomes/$ref/$g -p $outputD/filtered.blast.profile > $outputD/peg.counts");
	    }

	    if (! -s "$outputD/intergenic.counts")
	    {
		$ref =~ /^(\d+\.\d+)/;
		my $g = $1;
		&SeedUtils::run("perl rs_get_intergenic_counts.pl -d $dataD -g $dataD/ReferenceGenomes/$ref/$g -p $outputD/filtered.blast.profile > $outputD/intergenic.counts");
	    }
	    print STDERR "finished $outputD\n";
	}
    }
}

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3