Assume we have a `query.fa` file that contains sequences and we run:

`blat -stepSize=5 -repMatch=2253 -minScore=20 -minIdentity=0 -out=pslx /genomes/mm10.fa.qz query.fa output.pslx`

the output `output.pslx` file looks like this:


    match   mis-    rep.    N's     Q gap   Q gap   T gap   T gap   strand  Q               Q       Q       Q       T               T       T       T       block       blockSizes      qStarts  tStarts
            match   match           count   bases   count   bases           name            size    start   end     name            size    start   end     count
    ---------------------------------------------------------------------------------------------------------------------------------------------------------------
    20      0       0       0       0       0       0       0       +       seq     20      0       20      chr9    124595110       44046930        44046950   20,      0,      44046930,       aaaagtatcagtgtgtatag,   aaaagtatcagtgtgtatag,
    20      0       0       0       0       0       0       0       +       seq     20      0       20      chr9    124595110       44046930        44046950   20,      0,      44046930,       aaaagtatcagtgtgtatag,   aaaagtatcagtgtgtatag,


What would be a reasonable way to get the genomic contexts (5bp upsteam and 5bp downstream) for each aligned sequence. 

For example, assume that `blat` found that the seq: `AAATTGGGGAAAA` aligns to `chr2:100-113`, so the question is how to get `chr2:95-118` easily without reinventing the wheel. 

-------------------------

I couldn't make it work with bedtools, because my genome's index file is corrupted, but this should work for others who have successfully used `bwa` or `samtools` to index their reference genome:


    blat -stepSize=5 -repMatch=2253 -minScore=20 -minIdentity=0 -out=pslx /genomes/mm10.fa.qz query.fa output.pslx
    awk 'NR>5 {print $14 "\t" $16-10"\t" $17+10}' output.pslx > regions.bed
    bedtools getfasta -fi /genomes/mm10.fa.gz -bed regions.bed