--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/DEGseq_2.pl	Thu Sep 18 21:40:25 2014 -0400
@@ -0,0 +1,108 @@
+#!/usr/bin/perl -w
+#Filename:
+#Author: Tian Dongmei
+#Email: tiandm@big.ac.cn
+#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);
+}
+