| Previous changeset 0:bb1dbc0a1763 (2017-03-28) Next changeset 2:3ed885628c9f (2017-03-28) |
|
Commit message:
Uploaded |
|
added:
dataOverview.pl |
| b |
| diff -r bb1dbc0a1763 -r b66f4a551e25 dataOverview.pl --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/dataOverview.pl Tue Mar 28 21:56:04 2017 -0400 |
| [ |
| b'@@ -0,0 +1,392 @@\n+#!/usr/bin/perl -w\n+\n+#Margaret Antonio 16.08.29\n+\n+#use strict;\n+use Getopt::Long;\n+use Bio::SeqIO;\n+use autodie;\n+no warnings;\n+\n+\n+\n+#AVAILABLE OPTIONS. WILL print OUT UPON ERROR\n+sub print_usage() {\n+\n+ print "\\n###############################################################\\n";\n+ print "dataOverview: outputs basic statistics for tn-seq library files \\n\\n";\n+\n+ print "USAGE:\\n";\n+ print "perl dataOverview.pl -i inputs/ -f genome.fasta -r genome.gbk\\n";\n+ \n+ print "\\nREQUIRED:\\n";\n+ print " -d\\tDirectory containing all input files (results files from\\n";\n+ print " \\tcalc fitness script)\\n";\n+ print " \\t OR\\n";\n+ print " \\tIn the command line (without a flag), input the name(s) of \\n";\n+ print " \\tthe files containing fitness values for individual \\n\\tinsertion mutants\\n";\n+ print " -f\\tFilename for genome sequence, in fasta format\\n";\n+ print " -r\\tFilename for genome annotation, in GenBank format\\n";\n+\n+ print "\\nOPTIONAL:\\n";\n+ print " -h\\tprint OUT usage\\n";\n+ print " -c\\tCutoff average(c1+c2)>c. Default: 15\\n";\n+ print " -o\\tFilename for output. Default: standard output\\n";\n+ print " \\n~~~~Always check that file paths are correctly specified~~~~\\n";\n+ print " \\n###############################################################\\n";\n+\n+}\n+\n+# print "What\'s on the commandline: ", $ARGV;\n+\n+sub get_time() {\n+ my ($sec, $min, $hour, $mday, $mon, $year, $wday, $yday, $isdst) = localtime(time);\n+ return "$hour:$min:$sec";\n+ }\n+sub mean {\n+ my $sum=0;\n+ foreach my $n(@_){\n+ \t$sum+=$n;\n+ \t}\n+ my $total=scalar @_;\n+ my $mean=$sum/$total;\n+ return $mean; \t\n+}\n+sub minmax{\n+ my @unsorted=@_;\n+\tmy @sorted = sort { $a <=> $b } @unsorted;\n+\tmy $min = $sorted[0];\n+\tmy $max = $sorted[scalar @sorted -1];\n+\treturn ($min, $max);\n+\t}\n+sub uniq{\n+ my @input=@_;\n+ my @unique = do { my %seen; grep { !$seen{$_}++ } @input };\n+ }\n+\n+#ASSIGN INPUTS TO VARIABLES\n+our ($cutoff,$fastaFile, $outfile,$help,$ref,$weight_ceiling);\n+GetOptions(\n+\'r:s\' => \\$ref,\n+\'f:s\' => \\$fastaFile,\n+\'c:i\'=>\\$cutoff,\n+\'o:s\' => \\$outfile,\n+\'h\'=> \\$help,\n+\'w:i\' => \\$weight_ceiling,\n+);\n+\n+# Set defaults\n+#if (!$weight_ceiling){$weight_ceiling=50;}\n+#if (!$cutoff){$cutoff=10;}\n+\n+# If help option is specified or required files are not specified:\n+\n+if ($help) {\n+ print print_usage();\n+\tprint "\\n";\n+\texit;\n+}\n+\n+if (!$fastaFile or !$ref){\n+\tprint "\\nERROR: Please correctly specify reference genome fasta and genbank files\\n";\n+\tprint "Most genomes (in fasta and gbk format) are available at NCBI\\n";\n+ print print_usage();\n+\tprint "\\n";\n+\texit;\n+}\n+# Redirect STDOUT to log.txt. Anything print OUTed to the terminal will go into the log file\n+if (! $outfile){\n+\t$outfile="summary.txt";\n+}\n+\n+open OUT, ">",$outfile;\n+\n+#Not sure if I\'ll need this but sometimes funky data inputs have hidden characters\n+sub cleaner{\n+ my $line=$_[0];\n+ chomp($line);\n+ $line =~ s/\\x0d{0,1}\\x0a{0,1}\\Z//s;\n+ return $line;\n+}\n+\n+\n+#Get the input files out of the input directory, or take off of command line\n+\n+my @files=@ARGV;\n+foreach my $f(@files){\n+\t#print $f;\n+\t}\n+my $num=(scalar @files);\n+\n+#print OUT "Gathering data overview for Tn-Seq experiment\\n\\n";\n+#print OUT "Begin time: ",get_time(),"\\n\\n";\n+\n+#CREATE AN ARRAY OF DATA FROM INPUT CSV FILE(S). \n+#These are the "results" files from calc_fitness.pl. Insertion location, fitness, etc.\n+#Go through each file from the commandline (ARGV array) and read each line as an array\n+#into select array if values satisfy the cutoff\n+\n+\n+#Store ALL insertion locations in this array. Later, get unique insertions\n+my @insertions_all;\n+#Store all genes with valid insertions here\n+my @genes_insertions;\n+#all lines that satisfied cutoff\n+my @unsorted;\n+#array to hold all positions of insertions. Going to use this later to match up with TA sites\n+my @insertPos;\n+\n+#Markers\n+my $rows=-1;\n+my $last=0;\n+\n+print O'..b'GC content of this genome\\n";\n+print OUT "$ATcontent%\\tAT content of this genome\\n";\n+\n+print OUT "$satAll%\\tSaturation of TA sites before cutoff filter (allInsertions/TAsites)\\n";\n+print OUT "$sat%\\tSaturation of TA sites after cutoff filter (validInsertions/TAsites)\\n";\n+print OUT "$inscov%\\tGenome coverage by insertions (validInsertions/genomeSize)\\n";\n+print OUT "$tacov%\\tGenome coverage by TA sites (TAsites/genomeSize)\\n";\n+print OUT "$lg_dist_ta\\tLargest distance between TA sites\\n";\n+print OUT "$lg_dist_ins\\tLargest distance between insertions\\n";\n+print OUT "\\n\\nOpen Reading Frames\\n\\n";\n+\n+#Store everything to be print OUTed in array\n+my @table;\n+\n+#Find open reading frames from fasta file\n+local $_ = $fasta;\n+my @orfSize;\n+my @allc; #numbers of TAs in the ORFS here.\n+my $blank=0; #ORFS that don\'t have any TA sites.\n+my $orfCount=0; #keep track of the number of ORFs found.\n+my $minSize=0; \n+#Read somewhere that 99 is a good min but there is an annotated 86 bp gene for 19F\n+while ( /ATG/g ) {\n+ my $start = pos() - 3;\n+ if ( /T(?:AA|AG|GA)/g ) {\n+ my $stop = pos;\n+ my $size=$stop - $start;\n+ if ($size>=$minSize){\n+\t\t push (@orfSize,$size);\n+\t\t my $seq=substr ($_, $start, $stop - $start); \n+\t\t my @ctemp = $seq =~ /$x/g;\n+\t\t my $countTA = @ctemp;\n+\t\t if ($countTA==0){$blank++}\n+\t\t push (@allc,$countTA); \n+\t\t $orfCount++; \n+\t }\n+\t}\n+}\n+\n+print OUT "\\nORFs based on Fasta sequence and start (ATG) and end (TAA,TAG,TGA) codons\\n";\n+push (@table,["Set minimum size for an ORF",$minSize]);\n+print OUT "$orfCount\\tTotal number of ORFs found\\n";\n+my ($minORF, $maxORF) = minmax(@orfSize);\n+print OUT "$minORF\\tSmallest ORF\\n";\n+print OUT "$maxORF\\tLargest ORF\\n";\n+my ($mintaORF,$maxtaORF) = minmax(@allc);\n+print OUT "$mintaORF\\tFewest # TA sites in an ORF\\n";\n+print OUT "$maxtaORF\\tGreatest # TA sites in an ORF\\n";\n+print OUT "$blank\\tNumber of ORFs that don\'t have any TA sites\\n";\n+\n+\n+print OUT "\\nGenes using the genbank annotation file\\n\\n";\n+###Get genbank file. Find all start and stop for genes\n+#See how many insertions fall into genes vs intergenic regions\n+#Get array of coordinates for all insertions then remove insertion if it is\n+#within a gene region\n+my $gb = Bio::SeqIO->new(-file => $ref, -format => \'genbank\');\n+my $refseq = $gb->next_seq;\n+\n+#store number of insertions in a gene here\n+my @geneIns;\n+my @allLengths;\n+my $blankGene=0; #Number of genes that don\'t have any insertions in them\n+my @genomeSeq=split(\'\',$fasta);\n+\n+\n+#keep a copy to remove insertions that are in genes\n+my @insertPosCopy=@insertPos;\n+\n+my @features = $refseq->get_SeqFeatures(); # just top level\n+foreach my $feature ( @features ) {\n+\tif ($feature->primary_tag eq "gene"){\n+\t\tmy $start=$feature->start;\n+\t\tmy $end=$feature->end;\n+\t\tmy $length=$end-$start;\n+\t\tpush (@allLengths,$length);\n+\t\t#turn this into a for loop\n+\t\tmy $i=0;\n+\t\tmy $ins=0;\n+\t\tmy $current=$insertPos[$i];;\n+\t\twhile ($current<=$end && $i<scalar @insertPos){\n+\t\t\tif ($current>=$start){\n+\t\t\t\tsplice(@insertPosCopy, $i, 1);\n+\t\t\t\t$ins++;\t\t\t\n+\t\t\t}\n+\t\t\t$current=$insertPos[$i++];\n+\t\t}\n+\t\tif ($ins==0){$blankGene++}\n+\t\tpush (@geneIns,$ins);\n+\t}\n+}\n+my $avgLength=sprintf("%.2f",mean(@allLengths));\n+\n+my ($minLength, $maxLength) = minmax @allLengths;\n+my $avgInsGene=sprintf("%.2f",mean(@geneIns));\n+\n+\n+\n+\n+\n+my ($minInsGene, $maxInsGene) = minmax @geneIns;\n+my $nonGeneIns=scalar @insertPosCopy;\n+my $totalIns=scalar @insertPos;\n+my $percNon=sprintf("%.2f",($nonGeneIns/$totalIns)*100);\n+\n+print OUT "Length of a gene\\n";\n+print OUT "$avgLength\\tAverage","\\n$minLength\\tMininum","\\n$maxLength\\tMaximum\\n";\n+print OUT "\\nFor insertions in a gene:\\n";\n+print OUT "$avgInsGene\\tAverage","\\n$minInsGene\\tMininum","\\n$maxInsGene\\tMaximum\\n";\n+print OUT "Number of genes that do not have any insertions: ",$blankGene,"\\n";\n+print OUT "\\n$nonGeneIns\\tInsertions that are not in genes","\\n$percNon% of all insertions\\n";\n+#How many insertions are in genes and how many are in non-gene regions?\n+\n+\n+\n' |