text: user734: 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
  
  
  user776: Via BEDOPS convert2bed (psl2bed) and bedops operations:
  $ psl2bed < hits.psl | bedops --range 5 --everything - > answer.bed
  
  The file answer.bed will contain target intervals from the PSL (BLAT) input, padded up- and downstream by five bases.
  This BED file can be run through samtools faidx or similar to get sequence data.
  References:
  
  https://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/psl2bed.html
  
  https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#range-range
  
  
  
original_text: 734: <p>Assume we have a <code>query.fa</code> file that contains sequences and we run:</p>
  <p><code>blat -stepSize=5 -repMatch=2253 -minScore=20 -minIdentity=0 -out=pslx /genomes/mm10.fa.qz query.fa output.pslx</code></p>
  <p>the output <code>output.pslx</code> file looks like this:</p>
  <pre><code>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,
  </code></pre>
  <p>What would be a reasonable way to get the genomic contexts (5bp upsteam and 5bp downstream) for each aligned sequence.</p>
  <p>For example, assume that <code>blat</code> found that the seq: <code>AAATTGGGGAAAA</code> aligns to <code>chr2:100-113</code>, so the question is how to get <code>chr2:95-118</code> easily without reinventing the wheel.</p>
  <hr />
  <p>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 <code>bwa</code> or <code>samtools</code> to index their reference genome:</p>
  <pre><code>blat -stepSize=5 -repMatch=2253 -minScore=20 -minIdentity=0 -out=pslx /genomes/mm10.fa.qz query.fa output.pslx
  awk 'NR&gt;5 {print <span class="math-container">$14 "\t" $</span>16-10&quot;\t&quot; $17+10}' output.pslx &gt; regions.bed
  bedtools getfasta -fi /genomes/mm10.fa.gz -bed regions.bed
  </code></pre>
   776: <p>Via BEDOPS <code>convert2bed</code> (<code>psl2bed</code>) and <code>bedops</code> operations:</p>
  <pre><code>$ psl2bed &lt; hits.psl | bedops --range 5 --everything - &gt; answer.bed
  </code></pre>
  <p>The file <code>answer.bed</code> will contain target intervals from the PSL (BLAT) input, padded up- and downstream by five bases.</p>
  <p>This BED file can be run through <code>samtools faidx</code> or <a href="https://bioinformatics.stackexchange.com/a/5374/776">similar</a> to get sequence data.</p>
  <p>References:</p>
  <ul>
  <li><p><a href="https://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/psl2bed.html" rel="nofollow noreferrer">https://bedops.readthedocs.io/en/latest/content/reference/file-management/conversion/psl2bed.html</a></p>
  </li>
  <li><p><a href="https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#range-range" rel="nofollow noreferrer">https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html#range-range</a></p>
  </li>
  </ul>
  