Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-09-23T20:24:21.760Z
Answers.LastActivityDate: 2020-09-23T23:11:17.050Z
Answers.LastEditDate: 2020-09-23T23:11:17.050Z
Answers.UserDisplayName: Alex Reynolds
Answers.Body: <p>I couldn't get it to work with varied: <code>pcl2bed</code>, <code>samtools</code>, <code>bedtools</code> commands. I think this biopython should work:</p>
  <pre><code>import gzip
  from Bio import SeqIO
  n=5
  
  exact_matches = {}
  blat_qresult = SearchIO.read('blat_output.pslx', 'blat-psl', pslx=True)
  for hit in blat_qresult:
      for f in hit.fragments:
          if f.hit_span == f.query_span:
              chrom = hit.id
              if chrom not in exact_matches:
                  exact_matches[chrom] = []
              exact_matches[chrom].append(f.hit_range)
  
  
  with gzip.open('/genomes/mm10.fa.gz', &quot;rt&quot;) as handle:
      for record in SeqIO.parse(handle, &quot;fasta&quot;):
          if record.id in exact_matches:
              for subseq in exact_matches[record.id]:
                  print(record.seq[subseq[0]-n:subseq[1]+n])
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-09-03T22:28:34.873Z
Answers.LastActivityDate: 2020-09-03T22:28:34.873Z
Answers.UserDisplayName: 0x90
Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-09-03T17:55:04.303Z
Comments.Text: https://biopython.org/DIST/docs/api/Bio.SearchIO.BlatIO-module.html
Comments.UserDisplayName: 0x90
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-09-03T23:00:30.780Z
Comments.Text: Creat a bed file with `chr2 95 118` inside, and then use it with `bedtools getfasta` to extract your region
Comments.UserDisplayName: user3479780
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2020-09-03T23:01:44.633Z
Comments.Text: @user3479780 I tried that. Add some issue with indexing my reference genome. But this should work too. You are right.
Comments.UserDisplayName: 0x90
Tags: <python><bedtools><data-preprocessing><blat>
Title: How to get the genomic sequences from a blat result?
ThreadText: # How to get the genomic sequences from a blat result?
  
  (Asked by: username_0 on 2020-09-03)
  
  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
  
  ```
  
  **Tags:** python, bedtools, data-preprocessing, blat
  
  > **Comments:**
  > **username_0**: https://biopython.org/DIST/docs/api/Bio.SearchIO.BlatIO-module.html  
  > **username_1**: Creat a bed file with `chr2 95 118` inside, and then use it with `bedtools getfasta` to extract your region  
  > **username_0**: @username_1 I tried that. Add some issue with indexing my reference genome. But this should work too. You are right.  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_2 on 2020-09-23. Score: 1
  
  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](https://bioinformatics.stackexchange.com/a/5374/776) 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>
  
  
  ---
  
  ### Answer by username_0 on 2020-09-03. Score: 0
  
  I couldn't get it to work with varied: `pcl2bed`, `samtools`, `bedtools` commands. I think this biopython should work:
  
  ```
  import gzip
  from Bio import SeqIO
  n=5
  
  exact_matches = {}
  blat_qresult = SearchIO.read('blat_output.pslx', 'blat-psl', pslx=True)
  for hit in blat_qresult:
      for f in hit.fragments:
          if f.hit_span == f.query_span:
              chrom = hit.id
              if chrom not in exact_matches:
                  exact_matches[chrom] = []
              exact_matches[chrom].append(f.hit_range)
  
  with gzip.open('/genomes/mm10.fa.gz', "rt") as handle:
      for record in SeqIO.parse(handle, "fasta"):
          if record.id in exact_matches:
              for subseq in exact_matches[record.id]:
                  print(record.seq[subseq[0]-n:subseq[1]+n])
  
  ```
  
  
  ---
  
  