< Script to compute distribution length of sequences >

[Giorgio C]

Hi all,

i have a question for you:

Is there any tool to identify in a fasta file the read lenght distribution ?? I mean for exemple talking about miRNAs: " How many reads have 18 nucleotides, how many 19 nt, 20 nt....until 28 nt !!
Can you suggest me something also a little script please ???
At the end i need to do a graph showing this information.

Thank you very much

[kmcarr]

Originally Posted by Giorgio C

Hi all,

i have a question for you:

Is there any tool to identify in a fasta file the read lenght distribution ?? I mean for exemple talking about miRNAs: " How many reads have 18 nucleotides, how many 19 nt, 20 nt....until 28 nt !!
Can you suggest me something also a little script please ???
At the end i need to do a graph showing this information.

Thank you very much

In perl you can use the Statistics::Descriptive::Full module. Here is an example using the range you specified. It will read from a fasta file (named on the command line), fill an object with length data for each read and then print out a summary of that data, including a table of read length distribution suitable for plotting a histogram.

#!/usr/bin/perl

use warnings;
use strict;
use Statistics::Descriptive;

my $stat = Statistics::Descriptive::Full->new();
my (%distrib);
my @bins = qw/18 19 20 21 22 23 24 25 26 27 28/;

my $fastaFile = shift;
open (FASTA, "<$fastaFile");
$/ = ">";

my $junkFirstOne = <FASTA>;

while (<FASTA>) {
	chomp;
	my ($def,@seqlines) = split /\n/, $_;
	my $seq = join '', @seqlines;
	$stat->add_data(length($seq));
}


%distrib = $stat->frequency_distribution(\@bins);

print "Total reads:\t" . $stat->count() . "\n";
print "Total nt:\t" . $stat->sum() . "\n";
print "Mean length:\t" . $stat->mean() . "\n";
print "Median length:\t" . $stat->median() . "\n";
print "Mode length:\t" . $stat->mode() . "\n";
print "Max length:\t" . $stat->max() . "\n";
print "Min length:\t" . $stat->min() . "\n";
print "Length\t# Seqs\n";
foreach (sort {$a <=> $b} keys %distrib) {
	print "$_\t$distrib{$_}\n";
}

You can adjust the range and spacing of the distribution plot by adjusting the @bins array. See the documentation for the module.

[Giorgio C]

Amazing !!! Heartfelt thanks !!!

[Giorgio C]

Do you know how to modify the @bins array to know the total number of reads > of 28 nt and the tolat number of reads < of 18 nt....Is there possibile in one time ???

Thanks

[kmcarr]

From the module documentation for the frequency_distribution method:

The number of entries for a particular partition key are the number of items which are greater than the previous partition key and less then or equal to the current partition key.

This means that the bin for the lowest value (18 in this case) will count all occurrences which are ≤ 18 (so the program as written really won't return what you asked for in this case). Also for this case any sequences > 28nt in length will not be counted. Since you know that the value for a sequence length must be an integer you can add 17 as the first element of the list; this will count all reads where 0 < read length ≤ 17 (i.e. less than 18). To count any reads > 28nt you will need the last element in your @bin array to be greater than or equal to the longest possible read length in your data set; this could be 100, 1,000 or 1,000,000, it doesn't matter so long as the interval is large enough to capture any possible read length > 28nt in your data set.

[Giorgio C]

Infact it wasn't clear to me but now works fine for my purpose; Thank you very much again you have been very helpful !

[NicoBxl]

other bioperl script

#!/usr/bin/perl

#
# distribution length of a fasta file
#

use warnings;
use strict;
use Bio::SeqIO;

my $i = $ARGV[0];
my $o = $ARGV[1];

open (OUT, '>',$o);

my $seq_in  = Bio::SeqIO->new( -format => 'fasta',-file => $i);

my %length = ();
while( my $seq = $seq_in->next_seq() ) {
	my $seq1 = $seq->seq;	
	chomp($seq1);
	#print($seq1);
	my $l = length($seq1);
	#print $l."\n";
	my $tmp =  $length {$l};
	$length {$l} = $tmp+1;  
}
my @x = keys(%length);
my $sum = 0;
foreach my $k (@x) {
   $sum = $sum + $length{$k};
}

@x = sort{$a <=> $b}@x;
print(@x);
foreach my $k (@x) {
   print OUT $k."\t".$length{$k}."\t".$length{$k}*100/$sum."\n";
}

close (OUT);

[Giorgio C]

Thanks you, i'm trying it !!!

[mslider]

--hi

how to adjust @bins to print the report of all the length distribution ?

thank you --

Laurent --