match two files and check if position of chrom is within exon region

Question

I have two files, file1 and file2. I want to match file1 chrom with file2 chrom, and if there is a match then I want to check file2 if the position of the match chrom1 is within the exon region that is in file1. i want to do this with the whole file and print out rows that match and the position is within the exonic regions.

For instance, file1 has a chrom 1 and POS 736689 which match chrom 1 in file2, now I take POS 736689 from file1 and check if it is within the exon region 29553,30563,30975, 30039,30667,31097,... in file2.

If it is within the region it will print out the file1 row that matches the pattern.

The exon regions can be read as:

29553 to 30039
30563 to 30667
30975 to 31097

file1 look like this:

CHROM   POS REF ALT FATHER  MOTHER  DAUGHTER1   DAUGHTER2   DAUGHTER3   SON1    SON2    INFO
1   736689  T   A   0/0 1/0 0|0 0|1 0|1 0|0 0|1 AN=2184;AC=6    
1   758663  CGCC    TTAG    0/0 0/1 0|0 0|1 0|1 0|0 0|1 AC=2;AN=2184    
1   762005  G   C   1/0 0/0 1|0 1|0 1|0 1|0 0|0 AN=2184;AC=13   
1   771950  C   T   0/0 0/1 0|1 0|0 0|0 0|1 0|0 AC=110;AN=2184  
1   772172  C   T   0/1 0/0 0|0 0|0 0|0 0|0 1|0 AN=2184;AC=8    
1   790590  C   T   0/1 0/1 0|0 0/1 0/1 0|0 1|1 AN=2184;AC=84   

file2 look like this

chrom   Start   End
chr1    29553,30563,30975,  30039,30667,31097,
chr1    53048,54829,    53067,54936,
chr1    69090,  70008,
chr1    160445,161313,  160690,161525,
chr1    324755,327551,  326514,328453,
chr1    367639, 368634,
chr1    459655,461153,461750,   459992,461288,461954,
chr1    523008,529614,  523081,530148,

Answer 1

Despite an awk solution is favored by the OP, here is an R solution. (The task, especially dealing with exon coordinates in file2, is rather complex for a straightforward awk solution, at least for me).

res <- data.frame() # empty df for storing results

for (i in 1:nrow(file1)) {

  chr <- file1[i,1] # get chr
  pos <- file1[i,2] # get position of interest
  print("pos")
  print(pos)

  if(chr %in% file2$chrom) {

    file2_subset <- file2[file2$chrom == chr,] # subset file2

    for(j in 1:nrow(file2_subset)) {

      # create exon start - exon end pairs

      exon_start <- file2_subset[j,2]
      exon_start <- as.numeric(unlist(strsplit(exon_start,",")))
      print(exon_start)

      exon_end <- file2_subset[j,3]
      exon_end <- as.numeric(unlist(strsplit(exon_end,",")))
      print(exon_end)

      # if position of interest lie within exons
      # need to be careful if start/end position are counted as exon

      if(any(pos >= exon_start & pos <= exon_end)) {
        print(file1[i,])
        res <- rbind(res, file1[i,]) # write line i from file1 if satisfies above conds
      }
    }
  }
}

> res
  CHROM    POS REF ALT FATHER MOTHER DAUGHTER1 DAUGHTER2 DAUGHTER3 SON1 SON2          INFO
6  chr1 790590   C   T    0/1    0/1       0|0       0/1       0/1  0|0  1|1 AN=2184;AC=84

Answer 2

BEDTools already has a function to take a vcf and bedfile as input and look for intersections. Most people would get the exon coordinates into bed format, and use that.

Answer 3

Make a sorted BED file called exons.bed and run BEDOPS bedmap to map variants to exons:

$ bedmap --echo --echo-map --delim '\t' --skip-unmapped exons.bed <(vcf2bed < variants.vcf) > answer.bed

To make a sorted BED file from your exons, I'd suggest using sort-bed.

I would also suggest using awk to make sure that the chromosome naming scheme is consistent between exons and variants, e.g.:

$ vcf2bed < variants.vcf | awk '{ print 'chr$0' }' > variants.bed

Then use that result, instead, for mapping:

$ bedmap --echo --echo-map --delim '\t' --skip-unmapped exons.bed variants.bed > answer.bed

If chromosome names are not consistent, mapping or other set operations with BEDOPS or other toolkits will not work.