annotate DEGseq_2.pl @ 0:07745c0958dd draft

Uploaded

#!/usr/bin/perl -w
#Filename:
#Author: Tian Dongmei
#Email: [e-mail]
#Date: 2009-05-06
#Modified:
#Description: ɾ³ýmatched reads 
my $version=1.00;

use strict;
use Getopt::Long;
use File::Basename;

my %opts;
GetOptions(\%opts,"i=s","outdir=s","column1:i","mark1=s","depth1:i","depth2:i","column2:i","mark2=s","h");
if (!(defined $opts{i} and defined $opts{outdir} and defined $opts{mark1} and defined $opts{mark2}) || defined $opts{h}) { #necessary arguments
&usage;
}

my $filein=$opts{'i'};
my $outputdir=$opts{'outdir'};
unless ($outputdir=~/\/$/) {$outputdir .="/";}
my $column1=defined $opts{column1} ? $opts{column1} : 3;
my $column2=defined $opts{column2} ? $opts{column2} : 4;
my $mark1=$opts{mark1};
my $mark2=$opts{mark2};
my $fileout=$outputdir."degseq.R";
my $log=$outputdir."LOG.txt";

open OUT,">$fileout"; #output file  
open LOG,">$log";
print LOG "JOB start!\t";
print LOG `date`;
print LOG "\n";
#my ($name,$dir);
#$name=basename($filein);
print OUT "library(DEGseq)\n";
print OUT "geneExpFile <- system.file(package=\"DEGseq\")\n";
print OUT "geneExpFile<-file.path(\"$filein\")\n";
print OUT "layout(matrix(c(1,2,3,4,5,6), 3, 2, byrow=TRUE))\npar(mar=c(2, 2, 2,2))\n";
print OUT "outputdir<-file.path(\"$outputdir\")\n";
print OUT "geneExpMatrix1 <- readGeneExp(file=geneExpFile, geneCol=1, valCol=c($column1))\n";
print OUT "geneExpMatrix2 <- readGeneExp(file=geneExpFile, geneCol=1, valCol=c($column2))\n";
if(defined $opts{'depth1'} && defined $opts{'depth2'}){
print OUT "DEGexp(geneExpMatrix1=geneExpMatrix1, geneCol1=1, expCol1=c(2), groupLabel1=\"$mark1\",geneExpMatrix2=geneExpMatrix2, geneCol2=1, expCol2=c(2), groupLabel2=\"$mark2\",depth1=$opts{depth1},depth2=$opts{depth2},outputDir=outputdir,method=\"MARS\")\n";
}
else{
print OUT "DEGexp(geneExpMatrix1=geneExpMatrix1, geneCol1=1, expCol1=c(2), groupLabel1=\"$mark1\",geneExpMatrix2=geneExpMatrix2, geneCol2=1, expCol2=c(2), groupLabel2=\"$mark2\",outputDir=outputdir,method=\"MARS\")\n";
}
close OUT;

print LOG "Prepare for DEGseq!\t";
print LOG `date`;
print LOG "\n";

system("R CMD BATCH $fileout");

wait;

my $outfile=$outputdir."result.txt";
open OUT ,">$outfile";
my $deg=$outputdir."output_score.txt";
open IN,"<$deg";
my %hash;
while (my $aline=<IN>) {
	chomp $aline;
	if($aline=~/^\"/){print OUT "#GeneID\tchromsome\tvalue1\tvalue2\ttag\n";next;}
	my @temp=split/\t/,$aline;
	#$hash{$temp[0].$temp[1].$temp[2]}=$temp[$#temp];
	#my @tmp=split/\|/,$temp[0];
	#my @tmp=split/\:/,$temp[0];
	#my @po=split/\-/,$tmp[1];
	print OUT $temp[0],"\t",$temp[1],"\t",$temp[2],"\t",$temp[-1],"\n";
}
close IN;

#open IN,"<$filein";
#while (my $aline=<IN>) {
#	chomp $aline;
#	my @temp=split/\t/,$aline;
#	if (defined $hash{$temp[0].$temp[2].$temp[3]}) {print OUT $aline,"\t",$hash{$temp[0].$temp[2].$temp[3]},"\n";
#	}
#}

print LOG "Finish all JOB !\t";
print LOG `date`;
print LOG "\n";

sub usage{
print <<"USAGE";
Version $version
Usage:
$0 -i -outdir -column1 -mark1 -column2 -mark2 -depth1 -depth2
options:
-i input file
-outdir output file dir
-column1 the first column for DEGseq
-mark1 the name of the column1
-depth1 depth for the first file,use for normalize
-column2 the second column for DEGseq
-mark2 the name of the column2
-depth2 depth for the second file,use for normalize

-h help
USAGE
exit(1);
}