Finding Fasta Coordinates That Do Not Overlap

Question

[ashtx]

Hi all, I am very sorry if this question is very trivial. I have a list of fasta coordinates for the multiple sequence fasta file. For example

Id            start position       End position
1398         4                             8
1398         5                             10
1398         12                            15 
1756         25                            30
1756         28                            35

Is it possible to convert the top table as given below?

ID          start Position        End position
1398        4                               10
1398        12                              15
1756        25                              35

Here is my actual part of the raw data

So basically I am looking for a way that gives me non-overlapping fasta coordinates from the over lapping coordinates. I saw a very similar question on this link : http://www.biostars.org/post/show/7825/how-to-get-non-overlapping-coordinates-from-a-list-that-contains-overlapping-coordinates/ but it does not help in my case as I have unique IDs instead of chromosome numbers. I tried to convert it to bed format using sortBed but it seems to give me this error "Error: malformed GFF entry at line 1. Start was greater than end. Exiting"

I got these columns from Blast results which are based on the contig assembly. Now I want to form a sequence based on those blast results. I hope this question is clear.

[Leonor Palmeira]

Benm's solution from http://www.biostars.org/post/show/7825/how-to-get-non-overlapping-coordinates-from-a-list-that-contains-overlapping-coordinates/ works perfectly on your dataset. Did you try it?

[ashtx]

I tried but it should not report/output sub-region within the region. It means that, we still get region 5-10 even if it covered under region 4-12. According to his solution, it should take that into account but it does not seem so. I might be missing some tricks.

[Leonor Palmeira]

Could you give an example (editing your question for instance) of what output you get with Benm's solution? I get exactly the solution you describe in your question...

[ashtx]

Benm's solution works for some coordinates but then on some parts it does not work. Here is the input file http://pastie.org/4221091 Here is the output file http://pastie.org/4221095

[Vikas Bansal]

I think I understood. I edited my code and added if ($data[2] < $new[2]) . I think now it should work on your sorted file and will not report so called underlaps. Please have a look.

[ashtx]

Vikas, your solution works for some coordinates just as Benm's solution. Please see my input and output file in the previous reply. Thanks again.

[Vikas Bansal]

There are 2 problems in your input file because of which you are not getting desired results with my code. One is that the file is not sorted and for that I have already mentioned to use sort command on first 2 columns. Second is that you have negative strands also. Do you want to treat them separately or you want to merge them also.

Eg- If you have

829   4    10
829   9     3

You want to merge them or want to keep them separately?

[ashtx]

I want to merge them.

[Vikas Bansal]

Please see my edit and let me know if it works.

[ashtx]

It works pretty well. Thank you.

[ashtx]

My textfile does not contain any header line. I am using it to illustrate my question. Sorry about the confusion.

Answer 1

[Obi Griffith]

If you create a tab-delimited file from your first table, you can use the following R code to create a tab-delimited result file that looks like your second table:

library("GenomicRanges")
library("IRanges")

datadir="/your/data/dir/"
setwd(datadir)

#Import regions to exclude
Data2Mergefile="Data2Merge.txt"
Data2Merge=read.table(file=Data2Mergefile, sep="\t", header=TRUE)

#Create GenomicRanges object for manipulation
ranges2merge = GRanges(seqnames = Rle(Data2Merge[,"Id"]),
ranges = IRanges(Data2Merge[,"start position"], end = Data2Merge[,"End position"]), 
strand = Rle(strand("*"), length(rownames(Data2Merge) ) ) )

#Simplify ranges to merge overlaps
reduced=reduce(ranges2merge)

#Reformat and print to file
results=data.frame(Id=as.vector(seqnames(reduced)), Start=start(reduced), End=end(reduced))
write.table(results, file="merged.txt", sep="\t", row.names=FALSE, quote=FALSE)

[ashtx]

This works but I had some coordinates in reverse direction(i.e end position having lower numerical value than the start position) and I got this error, "Error in .Call2("solveuserSEW0", start, end, width, PACKAGE = "IRanges") : solving row 1: negative widths are not allowed".

To get around this problem, I sorted my file such that start position had lower numerical value than the end position. Thanks a lot.

[Obi Griffith]

Yes. by convention it is best if coordinates are specified in same direction unless you have a concept of +/- strand for your data. Then I believe + strand coordinates could be forward and - strand in reverse and GenomicRanges/IRanges would know what to do with it.

Answer 2

[Vikas Bansal]

EDIT - After OP's input file and comments.

First we will swap the start and end position, if start position is greater than end position and then we will sort this file.

Say this perl script is strand.pl -

use strict;
use warnings;
open my $file, '<',$ARGV[0] or die "Unable to open input file: $!";
my @data;

while (<$file>) {
@data = split;
  if ($data[1] > $data[2]) {

    my $temp = $data[1];
    $data[1] = $data[2];
    $data[2] = $temp;
  }

  print join("\t",@data),"\n";
}

Use -

perl strand.pl input.file |  sort -k1,1n -k2,2n > sorted.file

Now we will merge the overlapping positions.

Then use this perl script say merge.pl -

use strict;
use warnings;
open my $file, '<',$ARGV[0] or die "Unable to open input file: $!";
my @data;

while (<$file>) {

  if (not @data) {
    @data = split;
    next;
  }

  my @new = split;

  if ($new[0] == $data[0] and $new[1] <= $data[2] + 1) {

     if ($data[2] < $new[2]){
     $data[2] = $new[2];
      }
  }
  else {
    print join("\t", @data), "\n";
    @data = @new;
  }

  print join("\t", @data), "\n" if eof $file;

}

Use -

perl merge.pl sorted.file > result.file

[ashtx]

Sorry it does not work but I would love to see the solution using Perl.

[Vikas Bansal]

Can you please tell me, what error did you get while running Perl script on your sorted file?

[Vikas Bansal]

I have edited my code. There was error while pasting the code. I don't know in "while" condition "DATA" became "[HTML]". But now I have changed it to file.

[ashtx]

Now the code works but it does not solve my problem of getting rid of underlaps. For example if the coordinates for a particular ID are 1to147, it should get rid of other coordinates that are from 20-50 because I have already have those covered in 1-147. Thanks.

Answer 3

[Rm]

if I understood your Question correctly: 1.Add hash to the header line.

2.check to see if start positions ($2) is less than End position($3):

you can do that simply by using awk:

awk ' BEGIN {OFS = FS = "\t"};{if($2 < $3)print ; else print $1, $3, $2}'