1

I have 2 genetic datasets. I filter file1 based on a column in file2. However, I also need to account for a second column in file2 and I'm not sure how to do this.

The condition for file 1 line extraction is that only lines that have a chromosome position either more than 5000 larger or less than 5000 smaller than any given chromosome position in the lines of file 2 which are on the same chromosome are selected. I run part of this in the code below, however this does not account for only filtering lines within the same chromosome number.

For example my data looks like:

File 1:

Variant   Chromsome   Chromosome Position  
Variant1      2             14000     
Variant2      1             9000              
Variant3      8             37000          
Variant4      1             21000     

File 2:

Variant  Chromosome  Chromosome Position  
Variant1     1                 10000                   
Variant2     1                 20000                   
Variant3     8                 30000                   

Expected output (of variants with a greater than +/-5000 position distance in comparison to any line of file 2 on the same chromosome):

Variant   Chromosome Position     Chromosome
Variant1    14000                  2
Variant3    37000                  8

#Variant1 at 14000, whilst within 5000 + of Variant1 at 10000 in file2 is on a different chromosome and therefore not compared and is kept.
#Variant3 is on the same chromosome as Variant4 in file1 but larger than 5000+ distance and is kept.

I code using the answer to my previous question (How to select rows from a file based on conditions of another file?) asking how to filter without chromosome consideration:

awk '
NR == FNR       {RGMIN[++IX] = $2 - 5000                        
                 RGMAX[IX]   = $2 + 5000                        
                 next
                }
FNR == 1        {print                                         
                 next
                }
                {PR = 1                                         

                 for (i=2; i<=IX; i++)  PR = PR * ($3 < RGMIN[i] || $3 > RGMAX[i])                                                                                          
                }
PR                                                             
' file2 file1 > newfile.txt

I have been trying to understand the syntax better to add in an additional line with an if statement for "if line of file1 has matching chromosome number as file2, then perform the 5000 +/- filtering within that chromosome number only", but I'm not how to incorporate it into what I've got if there's another way.

Edit:

Now trying with CHROM additions:

awk '
NR == FNR       {RGMIN[++IX] = $2 - 500000                        
                 RGMAX[IX]   = $2 + 500000  
                 CHROM[IX] = $2              
                 next
                }
FNR == 1        {print                                         
                 next
                }
                {PR = 1                                         
                 for (i=2; i<=IX; i++) if ($2 == CHROM[i]) PR = PR * ($2 < RGMIN[i] || $2 > RGMAX[i])                                              
                }
PR          
' file2 file1 > newfile.txt               
  • You seem to have mixed up $2 and $3 in the first block. I'm away from kbd right now and can't test - try to add sth. like CHROM[IX] = $2 after RGMAX, and if ($2 == CHROM[i]) in the for loop. – RudiC Feb 6 at 21:40
  • Please can you confirm the condition(s) for inclusion, and the examples? At the very least, there's a typo in the example explanation ("Variant3 is on the same chromosome as Variant4 in file1"), and I'm struggling to understand why Variant1 2 should be included while Variant4 1 shouldn't. – JigglyNaga Feb 7 at 12:26
1

As proposed in comment:

awk '
NR == FNR       {RGMIN[++IX] = $3 - 5000                        
                 RGMAX[IX]   = $3 + 5000  
                 CHROM[IX]   = $2              
                 next
                }
FNR == 1        {print                                         
                 next
                }
                {PR = 1                                         
                 for (i=2; i<=IX; i++) if ($2 == CHROM[i]) PR = PR * ($3 < RGMIN[i] || $3 > RGMAX[i])                                              
                }
PR          
' file2 file1 
Variant   Chromsome   Chromosome Position  
Variant1      2             14000     
Variant3      8             37000          
0

Depending on how heavy file2 is in relation to file1 then you could also map file2 and then just test if the location in file1 is not already mapped.....

 awk 'FNR==NR{for (i=$3-5000; i<=$3+5000; i++){v[$2"."i]=1}}
      FNR!=NR{if (! v[$2"."$3]) print $0}
      ' file2 file1 

Slower than @RudiC for this data set but if there are large numbers of multiple overlapping ranges in file2 it may become more worthwhile.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.