| Previous changeset 5:873b49ef0e14 (2017-03-29) Next changeset 7:61e733c91a6b (2017-03-29) |
|
Commit message:
Uploaded |
|
added:
compStrains.pl |
| b |
| diff -r 873b49ef0e14 -r 90edb3f47005 compStrains.pl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/compStrains.pl Wed Mar 29 13:58:53 2017 -0400 |
| [ |
| b'@@ -0,0 +1,252 @@\n+#!/usr/bin/perl -w\n+\n+#Margaret Antonio 16.01.13\n+\n+#DESCRIPTION: Takes two aggregate.pl outputs and compares them using mean difference, pval for each\n+#gene. Can compare, for example, 19F in glucose and TIGR4 in glucose.\n+#DIFFERENT GENOMES (ie. diff. strains).\n+#Requires CONVERSION FILE\n+\n+#USAGE: perl compStrains.pl -c <conversion.csv> <options>\n+ #[<aggregateFile1.csv aggregateFile2.csv> OR -indir <indir/>]\n+\n+use Data::Dumper;\n+use strict;\n+use Getopt::Long;\n+use warnings;\n+use File::Path;\n+use File::Basename;\n+use Statistics::Distributions;\n+\n+#ASSIGN INPUTS TO VARIABLES USING FLAGS\n+our ($indir,$h,$out,$sortkey,$round,$l,$cfile);\n+GetOptions(\n+\'d:s\' => \\$indir,\n+\'h\' => \\$h,\n+\'o:s\' =>\\$out,\n+\'s:i\' => \\$sortkey,\n+\'r:i\'=> \\$round,\n+\'l1:s\'=> \\$l1,\n+\'l2:s\'=> \\$l2,\n+\'c:s\'=> \\$cfile,\n+);\n+\n+sub print_usage() {\n+ print "\\n";\n+ print "\\n##################################################################\\n";\n+ print "compStrains.pl: compare genes from a tn-seq experiment\\n";\n+ print "\\tfor two DIFFERENT strains/genomes using aggregate files\\n";\n+ \n+ print "\\nDESCRIPTION: Takes two aggregate.pl outputs and compares them by\\n";\n+ print "calculating the difference in mean fitness.\\n";\n+ \n+ print "Example: two strains tested under same condition.\\n";\n+ print "Note: For same strains (genomes), use compGenes.pl\\n";\n+ \n+ print "\\nUSAGE:\\n";\n+ print "perl compStrains.pl -c conversion.csv -d inputs/\\n";\n+ \n+ print "\\nREQUIRED:\\n";\n+ print " -d\\tDirectory containing all input files (files from\\n";\n+ print " \\taggregate fitness script)\\n";\n+ print " \\tOR\\n";\n+ print " \\tIn the command line (without a flag), input the name(s) of\\n";\n+ print " \\ttwo files containing aggregate gene fitness values. \\n\\n";\n+ print " -c\\tConversion file: two columns with homologs for both organisms\\n";\n+\n+ print "\\nOPTIONAL:\\n";\n+ print " -h\\tPrints usage and exits program\\n";\n+ print " -o\\tOutput file for comparison data. Default: label1label2.csv\\n";\n+ print " -s\\tSort output by this index of the file (indices begin at 0).\\n";\n+ print " \\tDefault: by mean\\n";\n+ print " -r\\tRound final output numbers to this number of decimals\\n";\n+ print " -l\\tLabels for input files. Default: filenames\\n";\n+ print " \\tTwo strings, comma separated (i.e. -l expt1,expt2).\\n";\n+ print " \\tOrder should match file order.\\n";\n+\tprint " \\n~~~~Always check that file paths are correctly specified~~~~\\n";\n+ print "\\n##################################################################\\n";\n+}\n+if ($h){\n+ print_usage();\n+ exit;\n+}\n+if (!$indir and (scalar @ARGV==0)){\n+\tprint "\\nERROR: Please correctly specify input files or directory\\n";\n+ print_usage();\n+\tprint "\\n";\n+\texit;\n+}\n+if (!$cfile){\n+\tprint "\\nERROR: Please correctly specify the required conversion file\\n";\n+ print_usage();\n+\tprint "\\n";\n+\texit;\n+}\n+\n+#THE @files ARRAY WILL CONTAIN INPUT FILE NAMES, EXTRACTED FROM A DIRECTORY (-indir) OR ARGV\n+my @files;\n+if ($indir){\n+ my $directory="$indir";\n+ opendir(DIR, $directory) or (print "Couldn\'t open $directory: $!\\n" and print_usage() and exit);\n+ my @direct= readdir DIR;\n+ my $tail=".csv";\n+ foreach (@direct){\n+ if (index($_, $tail) != -1){\n+ $_=$indir.$_;\n+ push (@files,$_);\n+ }\n+ }\n+ closedir DIR;\n+}\n+else{\n+ @files=@ARGV;\n+}\n+\n+#GET LABELS: USE (-l) OR USE FILNEAMES AS LABELS FOR COLUMNS IN OUTPUT FILE\n+\n+my @labels;\n+\n+my @labels = ($l1,$l2);\n+#if ($l){\n+# @labels=split(\',\',$l);\n+#}\n+#else{\n+# foreach (@files){\n+# my $filename=basename($_);\n+# my @temp=split(\'\\\\.\',$filename);\n+# my $colName=$temp[0];\n+# push (@labels,$colName);\n+ # }\n+#}\n+\n+#CHECK IF REQ. VARIABLES WERE DEFINED USING FLAGS. IF NOT THEN USE DEFAULT VALUES\n+\n+if (!$out) {$out="comp.".$labels[0].$labels[1].".csv"}\n+if (!$round){$round=\'%.4f\'}\n+\n+#OPEN INPUTTED AGGREGATE GENE FILES AND '..b'SH %ONE AND FILE2 GOES INTO HASH %TWO.\n+\n+#FILE1 OPENING ---> %one WHERE KEY:VALUE IS GENE_ID:(GENE_ID,INSERTIONS,MEAN,ETC.)\n+my @header;\n+my %one;\n+\n+open (F1,\'<\',$files[0]);\n+\n+#STORE COLUMN NAMES (FIRST LINE OF FILE1) FOR HEADER AND APPEND LABELS\n+my $head=<F1>; #the header in the file\n+my @cols=split(\',\',$head);\n+@cols=@cols[0,1,2,3,4,5,6]; #get rid of blank columns\n+for (my $j=0;$j<scalar @cols;$j++){\n+ $cols[$j]=$cols[$j].\'-\'.$labels[0]; #mark each column name with file it comes from\n+}\n+push (@header,@cols);\n+\n+while (my $line=<F1>){\n+ chomp $line;\n+ my @info=split(",",$line);\n+ #Only keep the first 7 columns (Ones about blanks aren\'t needed for comparisons)\n+ @info=@info[0,1,2,3,4,5,6];\n+ #Sometimes genes that don\'t have a gene name can\'t be blank, so fill with NA\n+ if (!$info[5]){\n+ $info[5]="NA";\n+ }\n+ #If there are no insertions in the column "total", then make it =0 rather than blank\n+ if (!$info[6]){\n+ $info[6]=0;\n+ }\n+ $one{$info[0]}=\\@info;\n+}\n+close F1;\n+\n+#FILE2 OPENING ---> %two WHERE KEY:VALUE IS GENE_ID:(GENE_ID,INSERTIONS,MEAN,ETC.)\n+\n+my %two;\n+open (F2,\'<\',$files[1]);\n+\n+#STORE COLUMN NAMES (FIRST LINE OF FILE2) FOR HEADER AND APPEND LABELS\n+$head=<F2>; #the header in the file\n+@cols=split(\',\',$head);\n+@cols=@cols[0,1,2,3,4,5,6]; #get rid of blank columns\n+for (my $j=0;$j<scalar @cols;$j++){\n+ $cols[$j]=$cols[$j].\'-\'.$labels[1]; #mark each column name with file it comes from\n+}\n+push (@header,@cols);\n+\n+while (my $line=<F2>){\n+ chomp $line;\n+ my @info=split(",",$line);\n+ @info=@info[0,1,2,3,4,5,6];\n+ if (!$info[5]){\n+ $info[5]="NA";\n+ }\n+ if (!$info[6]){\n+ $info[6]=0;\n+ }\n+ $two{$info[0]}=\\@info;\n+}\n+close F2;\n+\n+\n+#READ CONVERSION FILE INTO ARRAY.\n+#Conversion file must have strain 1 for file 1 in column 1 (index 0) and\n+ #strain 2 for file 2 in column 2 (index 1)\n+ #conversion file must be tab delimited with no NA fields\n+#If homologs (exist then take info from hashes (%one and %two) by referring to gene_id in KEY\n+\n+my @all; #store all homologs in this hash\n+open (CONV,\'<\',$cfile);\n+while (my $line=<CONV>){\n+ chomp $line;\n+ my @genes=split("\\t",$line); #Array @genes will contain two genes (SP_0000,SPT_0000)\n+ if (scalar @genes==2 and $genes[0] ne "" and $genes[1] ne ""){\n+ my @info;\n+ my @oneArray=@{$one{$genes[0]}};\n+ my @twoArray=@{$two{$genes[1]}};\n+ push (@info,@oneArray,@twoArray);\n+ my $diff=sprintf("$round",($info[1]-$info[8]));\n+ my $total1=$info[6];\n+ my $total2=$info[13];\n+ my $sd1=$info[3];\n+ my $se1=$info[4];\n+ my $sd2=$info[10];\n+ my $se2=$info[11];\n+ my $df=$total1+$total2-2;\n+ my $tdist;\n+ my $pval;\n+ #TDIST, PVAL calculations with fail if standard dev, error, or counts are not real numbers\n+ #or if 0 ends up in denominator\n+ if ($se1 eq "X" or $se2 eq "X" or $sd1 eq "X" or $sd2 eq "X" or $total1==0 or $total2==0 or $sd1==0 or $sd2==0){\n+ ($tdist,$pval)=("NA","NA");\n+ }\n+ else{\n+ $tdist=sqrt((($diff)/(sqrt((($sd1**2)/$total1)+(($sd2**2)/$total2))))**2);\n+ $pval=Statistics::Distributions::tprob($df,$tdist);\n+ }\n+ push (@info,$diff,$df,$tdist,$pval);\n+ push (@all,\\@info);\n+ }\n+}\n+close CONV;\n+\n+#SORT THE HOMOLOGS BY THE SORTKEY OR BY DEFAULT DIFFERENCE IN MEAN FITNESSES\n+if (!$sortkey){\n+ $sortkey=14; #for mean difference\n+}\n+my @sorted = sort { $b->[$sortkey] <=> $a->[$sortkey] } @all;\n+\n+#FINISH THE HEADER BY ADDING COLUMN NAMES FOR MEAN-DIFF, DOF, TDIST, AND PVALUE\n+my $field="MeanDiff(".$labels[0].\'.\'.$labels[1].")";\n+push (@header,$field,"DOF","TDIST","PVALUE");\n+\n+#PRINT MATCHED HOMOLOG INFORMATION INTO A SINGLE OUTPUT FILE\n+open OUT, \'>\',"$out";\n+print OUT (join(\',\',@header),"\\n");\n+foreach (@sorted){\n+ my @woo=@{$_};\n+ print OUT join(\',\',@woo),"\\n";\n+ }\n+\n+close OUT;\n+\n+\n' |