SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
bedtools ninja_help Bioinformatics 0 08-25-2015 09:03 AM
Bedtools partition mcgreevy Bioinformatics 2 10-11-2013 06:17 AM
bedtools cmccabe General 2 10-31-2012 11:54 AM
BEDTools version 2.10.0 quinlana Bioinformatics 0 09-21-2010 05:25 PM
BEDTools article quinlana Bioinformatics 6 03-03-2010 08:53 AM

Reply
 
Thread Tools
Old 09-25-2015, 08:45 AM   #1
sekhwal
Junior Member
 
Location: Morden

Join Date: Sep 2015
Posts: 9
Default Bedtools question

I am new with bedtools..So I have a perl script to compare three file in overlap manner, however my script is not perfectly improved to compare such thing.
I need to compare file 1 to 2, file 1 to 3 and file 2 to 3, I expect a common file for result...

file1:
AT4G01510.1 1 6993 7241
AT1G01020.2 1 7320 8668
file2:
AT5G35640.1 1 6993 7246
AT5G23980.1 1 7320 8669
file3:
AT1G01010.1 1 6993 7267
AT1G01020.1 1 7320 8634

#########################
Script:

#!/usr/bin/perl -w
use strict;
use Getopt::Std;
use vars qw ($opt_s $opt_p $opt_t);
getopts ('s:t:');

if(! $opt_s || !$opt_p || !$opt_t){
print "Usage: $0\n";
print "-s file1 \n";
print "-p file2 \n";
print "-t file3 \n";
exit;
}

my $tablefile1 = $opt_s;
my $tablefile2 = $opt_p;
my $tablefile3 = $opt_t;



if(!$tablefile1){
die "Invalid data file name $tablefile1.\n";
}
if(!$tablefile2){
die "Invalid data file name $tablefile2.\n";
}
if(!$tablefile3){
die "Invalid data file name $tablefile3.\n";
}

open(IN,$tablefile1) || die "Can't open $tablefile1..exiting.\n";

my %hash = ();
my @coord_arr = ();
my @chrs = ();
while(<IN>){
chomp;
next if (/^\s*$/);
my @cols = split(/\t/, $_);

my $id = $cols[0];

my $chr = $cols[1];
my $s = $cols[2];
my $e = $cols[3];

if (exists($hash{$chr})) {
my $str = $hash{$chr};
$hash{$chr} = "$str/$id:$s:$e";
} else {
$hash{$chr} = "$id:$s:$e";
}

}
close IN;

open(IN,$tablefile2) || die "Can't open $tablefile2..exiting.\n";

while(<IN>){
chomp;
next if (/^\s*$/);
my @cols = split(/\t/, $_);

my $id = $cols[0];

my $chr = $cols[1];
my $s = $cols[2];
my $e = $cols[3];

if (exists($hash{$chr})) {
my $str = $hash{$chr};
$hash{$chr} = "$str/$id:$s:$e";
} else {
$hash{$chr} = "$id:$s:$e";
}

}
close IN;

open(IN,$tablefile3) || die "Can't open $tablefile3..exiting.\n";
my $found = 0;
print "Chr\tPSF-ID\tStart\tEnd\tPseudopipe-ID\tStart\tEnd\tStart-diff\tend-diff\tShiu-ID\tStart\tEnd\tstart-diff\tend-diff\n";
while(<IN>){
chomp;
next if (/^\s*$/);
my @cols = split(/\t/, $_);

my $id = $cols[0];

my $chr = $cols[1];
my $s = $cols[2];
my $e = $cols[3];

my $info_str = $hash{$chr};
next if (!$info_str);

my @pseudos = split(/\//, $info_str);
#my @pseudos2 = split(/\//, $info_str);

for (my $i = 0; $i < @pseudos; $i++) {
my ($id1, $s1, $e1) = split(/:/, $pseudos[$i]);

if (&is_overlap($s, $e, $s1, $e1)) {
print "$chr\t$id\t$s\t$e\t$id1\t$s1\t$e1\t";
my $d1 = $s1-$s;
my $d2 = $e1-$e;
print "$d1\t$d2\t";
}

}
for (my $j = 0; $j < @pseudos; $j++) {
my ($id2, $s2, $e2) = split(/:/, $pseudos[$j]);

if (&is_overlap2($s, $e, $s2, $e2)) {
print "$id\t$s2\t$e2\t";
my $d3 = $s2-$s;
my $d4 = $e2-$e;
print "$d3\t$d4\n";
}
}

}

close IN;

sub is_overlap {
my ($s, $e, $s1, $e1) = @_;

my $found = 0;

if ($s >= $s1 && $s < $e1 ||
$s1 >= $s && $s1 < $e ||
$s1 >= $s && $e1 <= $e ||
$s1 <= $s && $e1 >= $e) {
$found = 1;
}

return $found;
}

sub is_overlap2 {
my ($s, $e, $s2, $e2) = @_;

my $found = 0;

if ($s >= $s2 && $s < $e2 ||
$s2 >= $s && $s2 < $e ||
$s2 >= $s && $e2 <= $e ||
$s2 <= $s && $e2 >= $e) {
$found = 1;
}

return $found;
}
sekhwal is offline   Reply With Quote
Old 09-25-2015, 09:14 AM   #2
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 6,975
Default

@sekhwal: Please don't post unrelated questions in threads that have the right words in title.

I have moved your post to a new thread.

Look into bedops or bedtools intersect option for your needs.

Last edited by GenoMax; 09-25-2015 at 09:17 AM.
GenoMax 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 06:54 PM.


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