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

View of /FigKernelScripts/rs_get_rnas.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 gjoseqlib;
use Getopt::Long;

my $usage = "usage: rs_get_rnas -d DataDir -r ReferenceGenome\n";
my $dataD;
my $refG;

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

(-s "$dataD/ReferenceGenomes/$refG") || die "$refG is not a reference genome";

my %rrnas;

foreach $_ (`cat $dataD/ReferenceGenomes/$refG/rRNA`)
{
    chomp;
    my $end;
    my($contig,$beg,$strand,$ln) = split(/\t/,$_);
    if ($strand eq "+")
    {
	$end = $beg + ($ln - 1);
    }
    else
    {
	my $b1 = $beg - ($ln - 1);
	my $e1 = $beg;
	($beg,$end) = ($b1,$e1);
    }
    push(@{$rrnas{$contig}},[$beg,$end]);
}

my $tot = 0;
my $hits = 0;

while (defined($_ = <STDIN>))
{
    my @flds = split(/\s+/,$_);
    $tot = $tot + 1;
    my $got = 0;
    my $contig = $flds[1];
    my $locs = $rrnas{$contig};
    if (! $locs) { $locs = [] }
    foreach $_ (@$locs)
    {
	my($y1,$y2) = @$_;

	my $x1;
	my $x2;
	if ($flds[8] < $flds[9])
	{
	    $x1 = $flds[8];
	    $x2 = $flds[9];
	}
	else
	{
	    $x1 = $flds[9];
	    $x2 = $flds[8];
	}

#	    print "x1=$x1 x2=$x2 y1=$y1 y2=$y2\n";
	if (! (($x1 > $y2) || ($y1 > $x2)))
	{
	    $got = 1;
	}
    }
    if ($got == 1)
    {
	$hits = $hits+1;
    }
    else
    {
	print join("\t",@flds),"\n";
    }
}
#print "hits=$hits tot=$tot\n";

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3