<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>

CC BY-SA 4.0
2020-09-23T20:24:21.760Z
2020-09-23T23:11:17.050Z
2020-09-23T23:11:17.050Z
Alex Reynolds
<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>

CC BY-SA 4.0
2020-09-03T22:28:34.873Z
2020-09-03T22:28:34.873Z
0x90
<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>

CC BY-SA 4.0
2020-09-03T17:55:04.303Z
https://biopython.org/DIST/docs/api/Bio.SearchIO.BlatIO-module.html
0x90
CC BY-SA 4.0
2020-09-03T23:00:30.780Z
Creat a bed file with `chr2 95 118` inside, and then use it with `bedtools getfasta` to extract your region
user3479780
CC BY-SA 4.0
2020-09-03T23:01:44.633Z
@user3479780 I tried that. Add some issue with indexing my reference genome. But this should work too. You are right.
0x90
<python><bedtools><data-preprocessing><blat>
How to get the genomic sequences from a blat result?
# 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**: @user3479780 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 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])


---

