thanks for reply. is BEDOPS graphics or command line operation? can you tell me, how I can give input of these three files to get compare in graphics.
I improved my script for comparing three files, now it is comparing file 1 to 2, file 1 to 3 but it is not able to compare file 2 to 3.
file 1
AT4G01510.1 1 6993 7241
AT1G01020.2 1 7320 8668
AT1G01050.1 1 31388 32672
file 2
AT1G43722.1 1 13306624 13307584
AT3G30520.1 1 13307748 13309520
AT5G59280.1 1 13330509 13331505
file 3
AT1G01010.1 1 3760 5627
AT1G01020.1 1 6918 7232
AT1G01020.1 1 8236 8666
#!/usr/bin/perl -w
use strict;
use Getopt::Std;
use vars qw ($opt_s $opt_p $opt_t);
getopts ('s:p:t:');
if(! $opt_s || !$opt_p || !$opt_t){
print "Usage: $0\n";
print "-s file1 output file \n";
print "-p file2 output file \n";
print "-t file3 output file \n";
exit;
}
my $tablefile1 = $opt_s;
my $tablefile3 = $opt_p;
my $tablefile2 = $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 $tablefile1..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\tfile1-ID\tStart\tEnd\tfile2-ID\tStart\tEnd\tStart-diff\tend-diff\tfile3-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;
}
1 answer
Looks like you are trying to do genomic set operations, for which you can use BEDOPS tools to accomplish this task easily and quickly.
Rearrange the columns in your files into BED-formatted files with awk, sort the BED files with sort-bed, then use bedmap to map one BED file against another.
Use the --echo-map-* options to decide what data you want from overlapping elements. You can add the --skip-unmapped option to effectively filter unmapped results, as well as adjust overlap stringency.
See the bedmap documentation for more details, or perhaps ask a second question with more detail (or revise this one) and I'll try to help.
But assuming you gave 10k or 1M genomic intervals in any/all your say A, B, C files, how exactly will the graphical comparison look like? Total bases overlapping, number of intervals overlapping by 1 base or something else?
Also after being recently bitten by a data file in which somebody got an idea of storing for some rows multiple values separated by comas, I am a bit paranoid about awk-ing instead of splitting lines say in Python, checking the number of columns, making sure it is not an empty one, etc. But I strongly second Alex opinion how to proceed once you got 2 BED files. Since there are well tested and optimized tools for processing/intersecting genomic intervals. one should not waste time reinventing them.
overlap concept algorithm is described in script in subroutine part.
I have no sufficient experience about BEDOPS. one of following application performing compare of multiple files, but it is not overlap concept. it is comparing exactly lines.
exactly it is genomic set operation and overlap concept..
bedtools intersect -a file1.txt -b file2.txt > output.txt
I am using this command to compare two files for overlap finding, the output file has only one input file data. How can I get overlap from both files in output file?
Try running:
bedtools intersect
without any arguments (it will print the options). One of them is:
-wo Write the original A and B entries plus the number of base
pairs of overlap between the two features
I am following the command for two files comparing with query, hence not getting results for both two files with query.
bedtools intersect -wa -wb -a query.bed -b file1.bed file2.bed >out2.txt
And also how I can specify number of base pairs of overlap in the command
thanks
bedtools intersect -wo -a query.bed -b file1.bed | filter_script.py > filtered_out.tsv
filter_script.py:
import sys
min_overlap = 10 #change the number
for line in sys.stdin:
sl = line.split()
overlap = atoi(sl[-1])
if overlap >= min_overlap:
print line,
In Python whitespace is significant, so you must keep the indentation for code blocs.
Log in to answer this question.
Could you say what your output should look like? That helps a lot.
Instead of assuming that anyone will bother going through your perl script, you could clearly state what exactly you want to compare..
I need to compare [4],[5],[6] columns in file1. [6], [7], [8] columns in file2 and [4], [5] and [6] in file3 and output looks like
Please be more specific. Are you comparing column 4 of file 1 to all the other columns of the other two files or what? Are you expecting equal values? Different values? etc.
Essentially you can:
So then you have a number of columns (sum of f1, f2 and f3) at stdout and you can pipe it to awk:
Where you replace
$numbersand with the column numbers that you want to compare (so e.g. the number of column 1 of f2 is number of columns of f1 + 1). Also, you might have to sort based on multiple columns. I wouldn't know.. the info you have given us is very vague..Hi, my main concern is to compare start and stop position number in above three files like column [5], [6] in file 1 to column [7], [8] in file 2 and column [5], [6] in file 1 to column [5], [6] in file 3. column [7], [8] in file 2 to column [5], [6] in file 3.
it is overlap concept not exact match number, please have a look in script for compare two file, the following concept I applied to compare two file as overlap
As others here have said, this question should be revised to clearly state what you are trying to accomplish, but in case you just need something quick and visual to compare 3 files at once, I'd like to add that there are several GUI options for 3-way ascii file comparisons such as kdiff3 and meld.