SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Illumina fragment length distribution delphi_ote Genomic Resequencing 3 05-18-2012 02:59 AM
Plotting length distribution of contigs? kga1978 Bioinformatics 2 01-16-2012 02:22 AM
Demultilplexing Sequences according to barcode --- Homemade script Giorgio C Bioinformatics 9 10-14-2011 08:17 AM
Script for calculating the average distribution of tags per gene from ChIP-seq?? Giles General 1 09-12-2011 03:25 AM
Compute paired-end distance distribution? krobison Bioinformatics 11 11-12-2009 09:30 AM

Reply
 
Thread Tools
Old 11-29-2011, 07:19 AM   #1
Giorgio C
Member
 
Location: ITALY

Join Date: Oct 2010
Posts: 89
Default < Script to compute distribution length of sequences >

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
Giorgio C is offline   Reply With Quote
Old 11-29-2011, 07:40 AM   #2
kmcarr
Senior Member
 
Location: USA, Midwest

Join Date: May 2008
Posts: 1,155
Default

Quote:
Originally Posted by Giorgio C View Post
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.

Code:
#!/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.
kmcarr is offline   Reply With Quote
Old 11-29-2011, 07:43 AM   #3
Giorgio C
Member
 
Location: ITALY

Join Date: Oct 2010
Posts: 89
Default

Amazing !!! Heartfelt thanks !!!
Giorgio C is offline   Reply With Quote
Old 11-29-2011, 08:01 AM   #4
Giorgio C
Member
 
Location: ITALY

Join Date: Oct 2010
Posts: 89
Default

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
Giorgio C is offline   Reply With Quote
Old 11-29-2011, 10:40 AM   #5
kmcarr
Senior Member
 
Location: USA, Midwest

Join Date: May 2008
Posts: 1,155
Default

From the module documentation for the frequency_distribution method:
Quote:
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.
kmcarr is offline   Reply With Quote
Old 11-30-2011, 01:10 AM   #6
Giorgio C
Member
 
Location: ITALY

Join Date: Oct 2010
Posts: 89
Default

Infact it wasn't clear to me but now works fine for my purpose; Thank you very much again you have been very helpful !
Giorgio C is offline   Reply With Quote
Old 11-30-2011, 01:47 AM   #7
NicoBxl
not just another member
 
Location: Belgium

Join Date: Aug 2010
Posts: 263
Default

other bioperl script

Code:
#!/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);
NicoBxl is offline   Reply With Quote
Old 11-30-2011, 06:05 AM   #8
Giorgio C
Member
 
Location: ITALY

Join Date: Oct 2010
Posts: 89
Default

Thanks you, i'm trying it !!!
Giorgio C is offline   Reply With Quote
Old 08-23-2012, 03:29 AM   #9
mslider
Junior Member
 
Location: france

Join Date: Sep 2010
Posts: 24
Default problem to setup @bins

--hi

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

thank you --

Laurent --
mslider is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 07:52 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2018, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO