countfasta.pl
Ce script écrit en perl permet de produire un histogramme des longueurs de séquences fasta et surtout de calculer quelques statistiques:
-nombre total de séquences
-longueur totale (taille cumulée de toutes les séquences)
-valeur de N25, N50 et N75
-Taux global de GC
#!/usr/bin/perl -w
# count_fasta.pl
# AUTHOR: Joseph Fass (modified from script by Brad Sickler)
# LAST REVISED: November 2010
#
# The Bioinformatics Core at UC Davis Genome Center
# https://bioinformatics.ucdavis.edu
# Copyright (c) 2009 The Regents of University of California, Davis Campus.
# All rights reserved.
use strict;
use POSIX;
use Getopt::Std;
my $usage = "\n\nusage: $0 [-i <interval size in # of residues>] <fasta file(s)>\n".
"Produces a histogram of sequence lengths and other measures to standard out.\n\n".
"-i # specify bin size for histogram (default 100)\n\n";
our($opt_i); # histogram interval
getopts('i:') or die $usage;
if (!defined($opt_i) or !($opt_i =~ m/^[0-9]+$/)) {$opt_i = 100;}
if( ( $#ARGV + 1 ) < 1 ) {
die $usage;
}
# Read in sequences from one or more fasta files
my @data_files;
for(my $i = 0; $i < ($#ARGV + 1); $i++){
$data_files[$i] = $ARGV[$i];
}
my $Id;
my %seq;
foreach my $file (@data_files){
open(FASTA, $file) or die"Can't open file $file\n";
while (<FASTA>) {
# if (/^>(.+)\s/) { $Id = $1; } # overwrites if id up to 1st whitespace is non-unique in file
if (/^>(.*)$/) { $Id = $1; }
# if (/(\w+)/) {
elsif (/^(\S+)$/) { $seq{$Id} .= $1 if $Id; }
}
close (FASTA);
}
# Count the number of sequences in the file and create a histogram of the distribution
my $n = 0;
my $int;
my $totalLength = 0;
my $gcCount = 0;
my %len= ();
my @seqLengths;
foreach my $id (keys %seq) {
push @seqLengths, length($seq{$id}); # record length for N50 calc's
$n++;
$int = floor( length($seq{$id})/$opt_i );
$totalLength += length($seq{$id});
$gcCount += ($seq{$id} =~ tr/gGcC/gGcC/);
if( !defined($len{$int}) ) {
$len{$int} = 1;
} else {
$len{$int}++;
}
}
# Calculate N25, N50, and N75 and counts
my $N25; my $N50; my $N75;
my $N25count=0; my $N50count=0; my $N75count=0;
my $frac_covered = $totalLength;
@seqLengths = reverse sort { $a <=> $b } @seqLengths;
$N25 = $seqLengths[0];
while ($frac_covered > $totalLength*3/4) {
$N25 = shift(@seqLengths);
$N25count++; $N50count++; $N75count++;
$frac_covered -= $N25;
}
$N50 = $N25;
while ($frac_covered > $totalLength/2) {
$N50 = shift(@seqLengths);
$N50count++; $N75count++;
$frac_covered -= $N50;
}
$N75 = $N50;
while ($frac_covered > $totalLength/4) {
$N75 = shift(@seqLengths);
$N75count++;
$frac_covered -= $N75;
}
# Print out the results
print "\n";
my @ints = sort { $a <=> $b } keys(%len);
for(my $i=$ints[0]; $i <= $ints[-1]; $i++) {
$len{$i} = 0 if(!defined($len{$i}));
printf "%d:%d \t$len{$i}\n", ( ($i*$opt_i), ($i*$opt_i+$opt_i-1) );
}
print "\nTotal length of sequence:\t$totalLength bp\n";
print "Total number of sequences:\t$n\n";
# not sure if these right wrt N25 and N75 ..
print "N25 stats:\t\t\t25% of total sequence length is contained in the ".$N25count." sequences >= ".$N25." bp\n";
print "N50 stats:\t\t\t50% of total sequence length is contained in the ".$N50count." sequences >= ".$N50." bp\n";
print "N75 stats:\t\t\t75% of total sequence length is contained in the ".$N75count." sequences >= ".$N75." bp\n";
print "Total GC count:\t\t\t$gcCount bp\n";
printf "GC %%:\t\t\t\t%.2f %%\n", ($gcCount/$totalLength * 100);
print "\n";