0

I have a file with 4 columns (file1.txt):

chr1    1156    G       G
chr1    1157    A       A
chr1    1165    T       T
chr1    1173    C       C
chr1    1175    G       G
chr1    1178    T       T
chr1    1181    C       C
chr1    1186    G       G

and a second file (file2.txt) with ranges, 2 columns:

1100    1160
1170    1180

I want to extract rows from file1 for which column 2 falls within ranges of file2. The desired output for the above example would be:

chr1    1156    G       G
chr1    1157    A       A
chr1    1173    C       C
chr1    1175    G       G
chr1    1178    T       T

Based on similar post I've tried, but it does not produced an output:

awk 'NR==FNR{ range[$1,$2]; next }{for(x in range) {split(x, check, SUBSEP); if($2>=check[1] && $2<=check[2]) print}} ' file2.txt file1.txt > output.txt

I have also tried the following with the same luck:

awk 'NR == FNR {ref[$1][$2]} if ($1 <= key && key <= $2) sum += ref[$2][key] print $0, sum} file2.txt file1.txt > output.txt

If anyone has any suggestion it would be appreciated it.

2
  • (1) “<a number> falls within the range of <two other numbers>” is notoriously vague / ambiguous. I see from the code that you tried that you want the ranges to be treated as inclusive, so, for example, if a line in file1.txt contains 1100, 1160, 1170 or 1180, you want it included in the output. When you say “between”, “within the range” or something similar, you should explain explicitly what you mean, demonstrate what you mean (by using boundary values like 1100, 1160, 1170 and 1180 in your sample data), or (ideally) both. … (Cont’d) Dec 22, 2021 at 19:59
  • (Cont’d) …  (2) You should either state that the ranges in file2.txt are discontiguous or specify what you want to happen if an input value (in file1.txt) falls within multiple ranges. (3) Please try to make your examples consistent — do the file names end with .txt or not? Dec 22, 2021 at 19:59

5 Answers 5

1

The following awk program should do:

awk 'NR==FNR{rng++;start[rng]=$1;end[rng]=$2;next}
     {for (i=1;i<=rng;i++) if (($2>=start[i])&&($2<=end[i])) {print; next}}' file2.txt file1.txt

It works as follows:

  • While parsing the first input file (file2.txt), indicated by NR, the global line-counter, being equal to FNR, the per-file line-counter, we register range start and end numbers in two arrays start and end (while counting the number of ranges in a counter rng). After that, we immediately skip execution to the next line.
  • While processing file1.txt (NR is now larger than FNR), we check, for each line, if column 2 falls in any range specified by corresponding entries in the start and end arrays, respectively. If so, we print the current line and again skip execution to the next line.
10
  • Hi AdminBee. Thanks for the suggestion. It does not seem to do what I need. The command runs, and got stuck printing the output file after an hour but it hasn't finished. I tried with a second file of ranges but it seems to produce the same results, and the ranges definitely do not overlap
    – Claudius
    Dec 3, 2021 at 15:20
  • @Claudius Interesting. Would it work on the example input you showed in your edited question? In my tests it worked, so at the moment I cannot explain the behavior you see. Question: You say it gets "stuck printing the output file" - how do you generate that file, simply by redirection ( > output.txt)?
    – AdminBee
    Dec 3, 2021 at 15:21
  • I made up the ranges for the example and it does actually! Then, there might be something wrong with my ranges file? Perhaps format? I prepared the ranges file using cut -f both columns. Perhaps I should use awk and specify output to be \t delimited
    – Claudius
    Dec 3, 2021 at 15:26
  • @Claudius It's possible. Although it is probably not the problem here: Was any of the files edited on a windows system? If so, the line endings may be wrong, and you have to run the file through dos2unix to correct that. Other than that: do you get any output at all which might indicate where the problem starts?
    – AdminBee
    Dec 3, 2021 at 15:30
  • Yes I get a nice output that make sense based on the ranges of the example. I'll double check my files. Many thanks because I know your command works.
    – Claudius
    Dec 3, 2021 at 15:32
0

both files are thousand rows long.

So, a couple thousand lines is really, for any computer of the last 30 years, not much data at all. Efficiency doesn't matter for you. (rough computation: 32 bytes per line in the first file, 16 per line in the second, so 48 bytes per line in total, say your computer can spare 2 GB of RAM before even blushing slightly, you can read 44 million lines into RAM without having to worry.)

Since this looks like genomics / bioinformatics, I think you stand a high chance to come into contact with python sooner or later, anyways.

#!/usr/bin/env python3
file1 = open("file1.txt", "r", encoding="ascii")
file2 = open("file2.txt", "r", encoding="ascii")

lines1 = file1.readlines()
lines2 = file2.readlines()

file1.close()
file2.close()

for dataline, rangeline in zip(lines1, lines2):
  splitrange = rangeline.split()
  lower = int(splitrange[0])
  upper = int(splitrange[1])
  
  ignore, valuestring, nucleotide1, nucleotide2 = dataline.split()
  value = int(valuestring)
  if lower <= value and value <= upper:
    print(dataline)

and that's it.

Is it as concise as doing it in awk? Certainly not. Is this as fast as it could be? No, not at all (but that won't matter). Will you stand a chance to remember what this did in a week? Most probably.

And as said, chances are you're going to do other stuff that's not much of what AWK is designed to do, anyways, so Python might be a natural tool. Almost certainly, learning about Biopython is a good idea.

2
  • Marcus. Many thanks for the suggestion. Just to clarify. Where would you specify the columns of each file in .py script? I don't see that info in the script. However, I tried to run it (after changing the file names) and got the following error: [s_menb@claudius regions]$ python3 subset_lines.py Traceback (most recent call last): File "subset_lines.py", line 17, in <module> value = int(valuestring) ValueError: invalid literal for int() with base 10: 'site'
    – Claudius
    Dec 3, 2021 at 13:22
  • split() splits a string into a Python list, and [i] selects the ith element. You can also "deconstruct" a four element list into four variables by saying a, b, c, d = four_element_list. Dec 3, 2021 at 13:37
0

Use two invocations of awk:

<file2.txt awk '{ print "$2 >= " $1 " && $2 <= " $2 }' |
           awk -f - file1.txt
5
  • Thor, thanks. Can't make this to work: <file2.txt awk '{ print "$2 >= " $1 " && $2 <= " $2 }' | awk -f - file1.txt > exon_GFE_2.txt Error: awk: awkgram.y:5020: find_line: Assertion `lineno > 0' failed. Shouldn't the command start with awk?
    – Claudius
    Dec 3, 2021 at 14:07
  • @Claudius: <filename should work in all bourne-like shells. What version of awk and shell are you using?
    – Thor
    Dec 3, 2021 at 15:00
  • awk 4.0.2 and shell version 4.2.46(2). I'm running this straight in the command line. Would it make any difference to run it as a bash?
    – Claudius
    Dec 3, 2021 at 15:08
  • 1
    That shell version sounds like you're already using bash, e.g. centos.pkgs.org/7/centos-x86_64/…, and that awk version sounds like it's an older version of GNU awk (supported by your use of gawk-only arrays-of-arrays syntax in your sample code - ref[$1][$2]). That error message sounds like you're hitting a bug in your version of gawk (see the assertion at gitlab.redox-os.org/redox-os/gawk/-/blob/…), get a newer one as we're now on gawk 5.1.1.
    – Ed Morton
    Dec 4, 2021 at 15:00
  • Having said that - creating an awk script from your input data and then asking awk to execute it is pretty dangerous and fragile as what the shell and/or awk can do in that situation will greatly depend on your input data so I personally wouldn't do it.
    – Ed Morton
    Dec 4, 2021 at 15:07
0

Using python we first generate a list of range objects from file2.txt and store in LoR variable. The idea being to test the presence of the target field if file1.txt (in our case the target field is second).

python3 -c 'import re, sys

## unpack command line arguments
if1,if2,col = sys.argv[1:]

idx = int(col) - 1
fx = lambda x: n in x
d2u = re.compile(r"\r?\n")

with open(if1) as f1, open(if2) as f2:
  LoR = [range(int(a),1+int(b))
    for _ in f2
    for l in [re.sub(d2u,"",_)]
    for a,b in [l.split()]]

  for l in f1:
    _ = re.sub(d2u,"",l)
    F = _.split()
    if len(F) < int(col): continue
    n = int(F[idx])
    if any(filter(fx,LoR)):
      print(_)
' file1.txt file2.txt 2

Same approach as above but with Perl:

perl -MList::Util=any -lane '
BEGIN {
  @LoR = map { [/(\S+)\s+(\S+)/] } <STDIN>;
  $idx = -1+pop;
}
s/\r$//; 
my $n = $F[$idx];
print if any { $_->[0] <= $n and $n <= $_->[1] } @LoR;
' file1.txt < file2.txt 2
7
  • Thanks for the code guest_7. I used the perl one. The output stops printing after a while but the run has not finished. I wonder if the command looks in the ranges file (file2.txt) for exact matches or not. For example the last printed line stopped in "1067833" value of column 2 (file1.txt). But in my file1.txt there seems to be another line with value "11067833". in my example file it would be like having 1156, but also 11156, 11560, 21156 and so on...which is the case for many many rows.
    – Claudius
    Dec 6, 2021 at 8:55
  • Without a representative data set we can only do so much. You should numeric sort on column 1 then column2 the ranges file (file2.txt) Validate it beforehand that it has the format as shown: all lines have exactly two fields, both numeric. The ranges mustn't overlap. What is number of lines of file2.txt ? Try the Python code to see if that helps.
    – guest_7
    Dec 6, 2021 at 17:31
  • Add the following line as the last line in the Perl code: BEGIN{use Data::Dumper;print Dumper(\@LoR);} then see whether the range file is correct.
    – guest_7
    Dec 6, 2021 at 18:18
  • Thanks for the suggestion. My ranges file2.txt has 303484 lines. I tried adding the additional line to the perl code but gets stuck again. I sorted files2.txt and there are some ranges contained in bigger ranges: 1 39/ 1 41/ 1 1006/ 1 256/ 1 393/ 1 66/ Would this be considered overlapping? I thought the code reads file2.txt line by line.
    – Claudius
    Dec 7, 2021 at 8:24
  • That is the problem. The overlapping ranges make the code spend unnecessary cycles leading to long times. You need to prune the range file properly like as shown in your example. They must NOT overlap.
    – guest_7
    Dec 7, 2021 at 11:54
-2
cat file2 |while read line ; do col1=$(echo $line| awk '{print $1}'); col2=$(echo $line | awk '{print $2}'); cat file1|while read fine; do echo $fine |awk -v col1="$col1" -v col2="$col2" '$2 >=col1 && $2 <col2'; done; done

output

chr1    1156    G       G
chr1    1157    A       A
chr1    1173    C       C
chr1    1175    G       G
chr1    1178    T       T
3
  • reason for downvote Dec 8, 2021 at 18:59
  • (1) Early versions of the question said “both files are thousand rows long.”  Even if you didn’t see that version of the question, you should have seen that statement quoted in Marcus Müller’s answer.  Well, if file1 has 14000 lines and file2 has 3000 lines, then your answer runs awk 42,006,000 times! (2) Using awk to break an input line into fields?  Why not say read col1 col2? … (Cont’d) Dec 22, 2021 at 20:00
  • (Cont’d) …  (3) If you look closely at the question, and read the code that the OP tried, you can see that he wanted ranges to be treated as inclusive, so, for example, if a line in file1 contains 1100, 1160, 1170 or 1180, it should be included in the output.  But your answer would get 1100 and 1170 but not 1160 and 1180, because you test $2 < col2 rather than $2 <= col2. Dec 22, 2021 at 20:00

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .