CGI/Perl Guide | Learning Center | Forums | Advertise | Login
Site Search: in

  Main Index MAIN
INDEX
Search Posts SEARCH
POSTS
Who's Online WHO'S
ONLINE
Log in LOG
IN

Home: Perl Programming Help: Beginner:
Compare two datasets

 



JackieZhou
New User

May 13, 2013, 1:58 PM

Post #1 of 7 (661 views)
Compare two datasets Can't Post

Hi,

I have two datasets about chromosome 1. data1 is a bed file, which has 3 columns, 1st column has 'chromosome', 2nd column has 'start coordinate', 3rd column has 'end coordinate':
------------------------------
chr1 100 600
chr1 800 1000
chr1 1200 3200
-----------------------------

Data2 has two columns, 1st column has the genomic coordinate information (indicates the base position. i.e., 5 means this is the 5th base on chr1), 2nd column is phastCons for each base.
-----------------------------------
755 0.7
756 0.9
767 0.3
800 0.2
801 0.1
-----------------------------------
So, for data2, the bases are not necessarily consecutive, because some regions may not get phastCons scores. As in the sample data I put above, the region between 767 and 800 are missing.

What I need to do now, is to the get PhastCons scores for the regions in data1. For example, in the region chr1:100 -- chr1:600, for these bases that have a phastCons, what is the maximum, min, sum and mean scores of the bases in this region?

I have written a program to do it, theoretically, it should work, but in practical, it is not going to, because it is super slow.

Can any one tell me how to do it efficiently?

Thank you!

perlChaser


Chris Charley
User

May 13, 2013, 4:32 PM

Post #2 of 7 (652 views)
Re: [JackieZhou] Compare two datasets [In reply to] Can't Post

Are all chromosomes going to be the same name (i.e. chr1 or chrY, etc)? How many regions mght there be?


perlchaser
Novice

May 13, 2013, 5:28 PM

Post #3 of 7 (644 views)
Re: [Chris Charley] Compare two datasets [In reply to] Can't Post


In Reply To
Are all chromosomes going to be the same name (i.e. chr1 or chrY, etc)? How many regions mght there be?


Thank you for replying. Yes, in this case, both files are about the same chromosome. there are around 20000 regions.


Chris Charley
User

May 13, 2013, 8:21 PM

Post #4 of 7 (632 views)
Re: [perlchaser] Compare two datasets [In reply to] Can't Post

The cleanest solution I was able to come up with is below. In it, I treated the phastcons file as a database so I could query it when reading from the BED file (of ranges).

Code
#!/usr/bin/perl 
use strict;
use warnings;
use List::Util qw/ min max sum /;
use DBI;

my $bed_file = 'sample.bed';
my $phastcon_file = 'phastcons.txt';

my $dbh = DBI->connect(qq{DBI:CSV:});
$dbh->{'csv_tables'}->{'scores'} = { 'file' => $phastcon_file,
'col_names' => [qw/number score/],
'sep_char' => " "};

my $sth = $dbh->prepare(<<SQL);
SELECT score
FROM scores
WHERE number >= ?
AND number <= ?
SQL

open BED, "<", $bed_file or die "Unable to open $bed_file. $!";
while (<BED>) {
my ($chr, $beg, $end) = split;
$sth->execute($beg, $end);
my $aref = $sth->fetchall_arrayref;
next unless my @scores = map $_->[0], @$aref;

printf "%s %s-%s: min: %.2f max: %.2f sum: %.2f mean: %.2f\n",
$chr, $beg, $end, min(@scores), max(@scores),
sum(@scores), sum(@scores) / @scores;
}
close BED or die "Unable to close $bed_file. $!";

__END__
*** sample.bed contents
chr1 100 600
chr1 700 1000
chr1 1200 3200

*** phastcons.txt contents
755 0.7
756 0.9
767 0.3
800 0.2
801 0.1
1100 .5
1101 .5
1200 .7

*** output
chr1 700-1000: min: 0.10 max: 0.90 sum: 2.20 mean: 0.44
chr1 1200-3200: min: 0.70 max: 0.70 sum: 0.70 mean: 0.70



perlchaser
Novice

May 13, 2013, 8:33 PM

Post #5 of 7 (626 views)
Re: [Chris Charley] Compare two datasets [In reply to] Can't Post

Hi Chris,

Thank you very much for the help. I will definitely try the code. Looks great!


BillKSmith
Veteran

May 13, 2013, 9:21 PM

Post #6 of 7 (624 views)
Re: [JackieZhou] Compare two datasets [In reply to] Can't Post

Here is a much different approach. I have no idea how fast it is. Memory could be a problem if there are too many regions. The number of bases affects time, but not memory.

Code
use strict; 
use warnings;
use Data::Dumper;
my @phastCons;
open my $SCORES, '<', 'file2.txt' or die "Cannot open file2.txt:$!";
while (<$SCORES>) {
chomp;
push @phastCons, [split /\s/, $_];
}
close $SCORES;

open my $BED, '<', 'file1.txt' or die "Cannot open file1.txt:$!";
open my $OUT, '>', 'file3.txt' or die "Cannot open file3.txt:$!";
while (<$BED>) {
chomp;
my @range = (split /\s/, $_)[1,2];
my $max=0;
my $min=99999;
my $sum=0;
my $count=0;
my $mean;
foreach (@phastCons) {
my ($pos, $score) = @$_;
if ($range[0] <= $pos and $pos <= $range[1]) {
$max = ($score > $max) ? $score : $max;
$min = ($score < $min) ? $score : $min;
$count++;
$sum += $score;
}
}
if ($count > 0){
$mean = $sum/$count;
print {$OUT} "$range[0] $range[1] $max $min $sum, $mean\n";
}
}
close $OUT;
close $BED;

Good Luck,
Bill


perlchaser
Novice

May 13, 2013, 9:51 PM

Post #7 of 7 (622 views)
Re: [BillKSmith] Compare two datasets [In reply to] Can't Post

Works perfectly. Thank you!

 
 


Search for (options) Powered by Gossamer Forum v.1.2.0

Web Applications & Managed Hosting Powered by Gossamer Threads
Visit our Mailing List Archives