Extract info from a column based on a range within another file

user_17566418 0
Viewed 75

I've tried looking around but what I found was looking up within the same file or combining columns with exact matches. Whereas I would not have exact matches and right now trying to combine these two codes is above my skill level. Basically I need to add an extra column to include the gene name based chromosome position and grabbing the gene name based on the range of the gene within another file. I know awk is my best bet possibly with FNR==NR.

File1 looks like this, where $1 is chromosome, $2 is position, the rest of the columns are sample coverage across that position:

chr1H   49525   47  41  60  74  93  34  117 
chr1H   49526   48  41  62  74  94  34  118 
chr1H   53978   48  40  61  73  94  33  117 
chr1H   53979   48  40  62  72  94  33  116

File2 looks like this, where $1 is the chromosome, $2 is the start of the gene $3 is the end of the gene and $4 is the gene name:

chr1H   49525   49772   gene1
chr1H   50194   50649   gene2   
chr1H   53978   54323   gene3   
chr1H   76743   77373   gene4

Either over writing or making a new file to end up with a file that looks like this:

chr1H   49525   47  41  60  74  93  34  117 gene1
chr1H   49526   48  41  62  74  94  34  118 gene1
chr1H   53978   48  40  61  73  94  33  117 gene3
chr1H   53979   48  40  62  72  94  33  116 gene3

Right now my code looks like this but I'm not sure how I specify the files (right now I've put in file1 or 2 so you know what my thinking is). So that the chromosomes match in both files and the position in the coverage file is between a range within start and end positions within the second file, then printing the entire line from file1 and the gene name from file2:

awk '{ if (file1$1 == file2$1 && file1$2 >= file2$2 && file1$2 <= file2$3) print file1$0, file2$4 }' file1 file2 > file3

Thanks for any help!

awk
3 Answers

You can do this fairly simply in awk by reading the range values from file2 into arrays indexed by gene name. That gives you a range by gene name to compare against the 2nd field in file1. You can do:

awk '
  NR == FNR {       # reading file2
    b[$4] = $2      # store b[] (begin) indexed by name
    e[$4] = $3      # store e[] (end) indexed by name
    next            # get next record
  } 
  {                 # for all file1 records
    for(i in b) {   # loop over values by gene name
      if ($2 >= b[i] && $2 <= e[i]) {   # if in range b[] to e[]
        printf "%s %s\n", $0, i         # output with gene name at end
        next                            # get next record
      }
    }
  }
' file2 file1

Example Use/Output

With the values shown in file1 and file2 you would have:

$ awk '
>   NR == FNR {       # reading file2
>     b[$4] = $2      # store b[] (begin) indexed by name
>     e[$4] = $3      # store e[] (end) indexed by name
>     next            # get next record
>   }
>   {                 # for all file1 records
>     for(i in b) {   # loop over values by gene name
>       if ($2 >= b[i] && $2 <= e[i]) {   # if in range b[] to e[]
>         printf "%s %s\n", $0, i         # output with gene name at end
>         next                            # get next record
>       }
>     }
>   }
> ' file2 file1
chr1H   49525   47  41  60  74  93  34  117  gene1
chr1H   49526   48  41  62  74  94  34  118  gene1
chr1H   53978   48  40  61  73  94  33  117  gene3
chr1H   53979   48  40  62  72  94  33  116  gene3
user_34221020

You can try this solution. Although I would suggest to use e.g. bedtools for stuff like this. You can encounter situations where intuition fails and it's good to have stress tested tools when that happens.

$ awk 'NR==FNR{ chr[NR]=$1; st[NR]=$2; en[NR]=$3; gene[NR]=$4; x=NR } 
    NR!=FNR{ for(i=1;i<=x;i++){ if($1==chr[i] && $2>=st[i] && $2<=en[i]){ 
    print $0,gene[i]; break } } }' f2 f1
chr1H   49525   47  41  60  74  93  34  117  gene1
chr1H   49526   48  41  62  74  94  34  118  gene1
chr1H   53978   48  40  61  73  94  33  117  gene3
chr1H   53979   48  40  62  72  94  33  116  gene3

Btw, this doesn't handle multiple matches. If you want to allow these you can remove the break statement and change the print to a printf.

user_94620950

You can use join to get half way there:

sort -k 2,2 file-1 > f1-sorted
sort -k 2,2 file-2 > f2-sorted
join -1 2 -2 2 -o 1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.4 f1-sorted f2-sorted

#gives
chr1H 49525 47 41 60 74 93 34 117 gene1
chr1H 53978 48 40 61 73 94 33 117 gene3
  • join joins files on a column, but it must be sorted
  • -1 2 -2 2 means join on the second column of file 1 and file 2
  • -o specifies output format for the columns (1.1 is file 1, column 1, 2.4 is file 2, column 4)

But you have an unusual requirement of matching the next immediate position number with the previous gene name. For that I would use this awk:

awk '
FNR==NR {name[$2]=$4}
FNR!=NR && name[$2] {n = name[$2]}
FNR!=NR && name[$2-1] {n = name[$2-1]}
FNR!=NR {print $0,n}' file-2 file-1

This gives your expected output exactly. You can also use FNR!=NR && n {print $0,n}, to only print a record if there's a match for the position number (column 2) in both files.

user_139196680
Related