<p>Usually, when you run your BLAST alignment you specify the output format (e.g. as described in &quot;outfmt&quot; in <a href="https://www.ncbi.nlm.nih.gov/books/NBK279684/table/appendices.T.options_common_to_all_blast/" rel="nofollow noreferrer">https://www.ncbi.nlm.nih.gov/books/NBK279684/table/appendices.T.options_common_to_all_blast/</a>). If, for some reason, you can't re-run the BLAST command (e.g. if this output file was given to you, or if you can't access your query or reference sequence/s) you <em>can</em> parse this 'readable' output format but it's not straightforward.</p>
<p>Here is a potential solution using awk:</p>
<pre><code>cat positions.txt
129
249
366
372
392

cat input.txt
&gt; CD630_08510_19
Query  121  TGGCTTGGGTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  121  TGGCTTGGTTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180

Query  241  GCAATTGTGCTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  241  GCAATTGTACTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300

Query  361  CCATACTTTATAGGTATGTTTATTGGTATGGTTGGATAG  399
            ||||| ||||| ||||||||||||||||||| |||||||
Sbjct  361  CCATATTTTAT-GGTATGTTTATTGGTATGGCTGGATAG  398
&gt; CD630_08510_20
Query  121  TGGCTTGGATTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  121  TGGCTTGGGTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180

Query  241  GCAATTGTGCTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  241  GCAATTGTACTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300

Query  361  CCATACTTTATAGGTATGTTTATTGGTATGGTTGGATAG  399
            ||||| ||||| ||||||||||||||||||| |||||||
Sbjct  361  CCATATTTTAT-GGTATGTTTATTGGTATGGCTGGATAG  398

awk 'BEGIN {
    OFS = &quot;\t&quot;
}

NR == FNR {
    a[$0]
    next
}

<span class="math-container">$1 == "&gt;" {
b[++c] = $</span>2
}

<span class="math-container">$1 != "&gt;" {
for (i in a) {
    if ($1 == "Query" &amp;&amp; i &gt;= $2 &amp;&amp; i &lt;= $4) {
        d[++e] = b[c] OFS i OFS substr($3, (i - $2) + 1, 1)
    } else if ($1 == "Sbjct" &amp;&amp; i &gt;= $2 &amp;&amp; i &lt;= $4) {
        f[++g] = substr($3, (i - $2) + 1, 1)
        h[c] = h[c] substr($3, (i - $</span>2) + 1, 1)
        }
    }
}

END {
    printf &quot;%s\n&quot;, &quot;Seq_name\tbase_position\tquery_base\tsubject_base\tall_subject_base&quot;
    for (j = 1; j &lt;= e; j++) {
        printf &quot;%s\t%s\t%s\n&quot;, d[j], f[j], (j % length(a) == 1 ? h[int(j / length(a)) + 1] : &quot;&quot;)
    }
}' positions.txt input.txt
Seq_name    base_position   query_base  subject_base    all_subject_base
CD630_08510_19  129 G   T   TA-CT
CD630_08510_19  249 G   A
CD630_08510_19  372 A   -
CD630_08510_19  392 T   C
CD630_08510_19  366 C   T
CD630_08510_20  129 A   G   GA-CT
CD630_08510_20  249 G   A
CD630_08510_20  372 A   -
CD630_08510_20  392 T   C
CD630_08510_20  366 C   T

</code></pre>
<p>If you have more positions:</p>
<pre><code>cat positions.txt
129
249
366
372
373
374
375
376
392

awk 'BEGIN {
    OFS = &quot;\t&quot;
}

NR == FNR {
    a[$0]
    next
}

<span class="math-container">$1 == "&gt;" {
b[++c] = $</span>2
}

<span class="math-container">$1 != "&gt;" {
for (i in a) {
    if ($1 == "Query" &amp;&amp; i &gt;= $2 &amp;&amp; i &lt;= $4) {
        d[++e] = b[c] OFS i OFS substr($3, (i - $2) + 1, 1)
    } else if ($1 == "Sbjct" &amp;&amp; i &gt;= $2 &amp;&amp; i &lt;= $4) {
        f[++g] = substr($3, (i - $2) + 1, 1)
        h[c] = h[c] substr($3, (i - $</span>2) + 1, 1)
        }
    }
}

END {
    printf &quot;%s\n&quot;, &quot;Seq_name\tbase_position\tquery_base\tsubject_base\tall_subject_base&quot;
    for (j = 1; j &lt;= e; j++) {
        printf &quot;%s\t%s\t%s\n&quot;, d[j], f[j], (j % length(a) == 1 ? h[int(j / length(a)) + 1] : &quot;&quot;)
    }
}' positions.txt input.txt
Seq_name    base_position   query_base  subject_base    all_subject_base
CD630_08510_19  129 G   T   TA-GGTACT
CD630_08510_19  249 G   A
CD630_08510_19  372 A   -
CD630_08510_19  373 G   G
CD630_08510_19  374 G   G
CD630_08510_19  375 T   T
CD630_08510_19  376 A   A
CD630_08510_19  392 T   C
CD630_08510_19  366 C   T
CD630_08510_20  129 A   G   GA-GGTACT
CD630_08510_20  249 G   A
CD630_08510_20  372 A   -
CD630_08510_20  373 G   G
CD630_08510_20  374 G   G
CD630_08510_20  375 T   T
CD630_08510_20  376 A   A
CD630_08510_20  392 T   C
CD630_08510_20  366 C   T
</code></pre>
<p>NB this won't work if you only have a single position in positions.txt, but if that is a potential use-case you can modify the script to suit, e.g. replace <code>(j % length(a) == 1 ? h[int(j / length(a)) + 1] : &quot;&quot;)</code> with <code>f[j]</code>.</p>

CC BY-SA 4.0
2023-03-31T16:45:58.123Z
Jared_mamrot, Thank you so much.
user14714429
CC BY-SA 4.0
2023-03-27T22:28:51.427Z
2023-03-27T22:28:51.427Z
jared_mamrot
<p>I have an input file looks like this:</p>
<pre><code>&gt; CD630_08510_19
Query  121  TGGCTTGGGTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  121  TGGCTTGGTTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180

Query  241  GCAATTGTGCTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  241  GCAATTGTACTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300

Query  361  CCATACTTTATAGGTATGTTTATTGGTATGGTTGGATAG  399
            ||||| ||||| ||||||||||||||||||| |||||||
Sbjct  361  CCATATTTTAT-GGTATGTTTATTGGTATGGCTGGATAG  398
</code></pre>
<p>The input file is the result file after running the command line blast.
It's a txt file generated after blast alignment and in txt file header format are like this.</p>
<p>I am trying to extract the specific nucleotide base for position 129, 249, 366, 372 and 392 from query and subject aligned sequence. python and perl solution will be fine as well. thank you.</p>
<p>For example: I need my output to have format below like this:</p>
<pre><code>Seq_name        base_position       query_base  subject_base    all_subject_base
CD630_08510_19  129                 G           T               TAT-C
CD630_08510_19  249                 G           A
CD630_08510_19  366                 C           T
CD630_08510_19  372                 A           -
CD630_08510_19  392                 T           C

CD630_08510_20  129                 G           G               GATCC
CD630_08510_20  249                 G           A
CD630_08510_20  366                 C           T
CD630_08510_20  372                 A           C
CD630_08510_20  392                 T           C
</code></pre>
<p>all_subject_base: are the bases from position 129, 249, 366, 372 and 392 of subject_base column. I want to keep this in one column thats why all_subject_base.</p>
<p>I tried to extract this with grep but have not been able to. I am new in bioinformatics so any suggestion would be very helpful.</p>

CC BY-SA 4.0
2023-03-23T14:42:57.007Z
Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.
Community
CC BY-SA 4.0
2023-03-23T16:21:58.813Z
Have you tried looking for ways to get [alternative output formats](https://www.ncbi.nlm.nih.gov/books/NBK569862/) from blast that would be easier to parse?
Maximilian Press
CC BY-SA 4.0
2023-03-23T18:33:34.597Z
I agree with @MaximilianPress given the approach used. The thing about Biopython is this already done - so it's easy
M__
CC BY-SA 4.0
2023-03-23T19:48:43.740Z
FYI, the BLAST text format output is for reading, not parsing with a computer program, and is subject to change without notice. I don't recommend using it, especially if you want to do similar things in the future.
gringer
CC BY-SA 4.0
2023-03-24T12:33:41.323Z
Will you always have just one alignment in your file or can you have muiltiple? Does our solution need to deal with both cases?
terdon
CC BY-SA 4.0
2023-03-24T12:56:22.727Z
And why aren't your numbers in sequence? Why do we go from 180 as the last position of the first line to 241 as the first position of the second? Why isn't that 181? Is that actually what your data look like?
terdon
<blast><linux><shell>
Extract specific nucleotide base position for query and subject from text file
# Extract specific nucleotide base position for query and subject from text file

(Asked by: username_0 on 2023-03-23)

I have an input file looks like this:
> CD630_08510_19
Query  121  TGGCTTGGGTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  121  TGGCTTGGTTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180

Query  241  GCAATTGTGCTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  241  GCAATTGTACTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300

Query  361  CCATACTTTATAGGTATGTTTATTGGTATGGTTGGATAG  399
            ||||| ||||| ||||||||||||||||||| |||||||
Sbjct  361  CCATATTTTAT-GGTATGTTTATTGGTATGGCTGGATAG  398

The input file is the result file after running the command line blast.
It's a txt file generated after blast alignment and in txt file header format are like this.
I am trying to extract the specific nucleotide base for position 129, 249, 366, 372 and 392 from query and subject aligned sequence. python and perl solution will be fine as well. thank you.
For example: I need my output to have format below like this:
Seq_name        base_position       query_base  subject_base    all_subject_base
CD630_08510_19  129                 G           T               TAT-C
CD630_08510_19  249                 G           A
CD630_08510_19  366                 C           T
CD630_08510_19  372                 A           -
CD630_08510_19  392                 T           C

CD630_08510_20  129                 G           G               GATCC
CD630_08510_20  249                 G           A
CD630_08510_20  366                 C           T
CD630_08510_20  372                 A           C
CD630_08510_20  392                 T           C

all_subject_base: are the bases from position 129, 249, 366, 372 and 392 of subject_base column. I want to keep this in one column thats why all_subject_base.
I tried to extract this with grep but have not been able to. I am new in bioinformatics so any suggestion would be very helpful.

**Tags:** blast, linux, shell
> **Comments:**
> **username_1**: Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.  
> **username_2**: Have you tried looking for ways to get [alternative output formats](https://www.ncbi.nlm.nih.gov/books/NBK569862/) from blast that would be easier to parse?  
> **username_3**: I agree with @MaximilianPress given the approach used. The thing about Biopython is this already done - so it's easy  
> **username_4**: FYI, the BLAST text format output is for reading, not parsing with a computer program, and is subject to change without notice. I don't recommend using it, especially if you want to do similar things in the future.  
> **username_5**: Will you always have just one alignment in your file or can you have muiltiple? Does our solution need to deal with both cases?  
> **username_5**: And why aren't your numbers in sequence? Why do we go from 180 as the last position of the first line to 241 as the first position of the second? Why isn't that 181? Is that actually what your data look like?  

---

## 1 answers

### Answer by username_6 on 2023-03-27. Score: 2

Usually, when you run your BLAST alignment you specify the output format (e.g. as described in "outfmt" in https://www.ncbi.nlm.nih.gov/books/NBK279684/table/appendices.T.options_common_to_all_blast/). If, for some reason, you can't re-run the BLAST command (e.g. if this output file was given to you, or if you can't access your query or reference sequence/s) you can parse this 'readable' output format but it's not straightforward.
Here is a potential solution using awk:
cat positions.txt
129
249
366
372
392

cat input.txt
> CD630_08510_19
Query  121  TGGCTTGGGTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  121  TGGCTTGGTTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180

Query  241  GCAATTGTGCTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  241  GCAATTGTACTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300

Query  361  CCATACTTTATAGGTATGTTTATTGGTATGGTTGGATAG  399
            ||||| ||||| ||||||||||||||||||| |||||||
Sbjct  361  CCATATTTTAT-GGTATGTTTATTGGTATGGCTGGATAG  398
> CD630_08510_20
Query  121  TGGCTTGGATTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  121  TGGCTTGGGTTTATATTTGATACTTTAGGAACATTTACAATGAGTAACATAGCCAACAGC  180

Query  241  GCAATTGTGCTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300
            |||||||| |||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  241  GCAATTGTACTTATGCTTTTTCATGCAAGTTGGGCTACTTTTGTCCTATACAAAGATGAT  300

Query  361  CCATACTTTATAGGTATGTTTATTGGTATGGTTGGATAG  399
            ||||| ||||| ||||||||||||||||||| |||||||
Sbjct  361  CCATATTTTAT-GGTATGTTTATTGGTATGGCTGGATAG  398

awk 'BEGIN {
    OFS = "\t"
}

NR == FNR {
    a[$0]
    next
}

$1 == ">" {
b[++c] = $2
}

$1 != ">" {
for (i in a) {
    if ($1 == "Query" && i >= $2 && i <= $4) {
        d[++e] = b[c] OFS i OFS substr($3, (i - $2) + 1, 1)
    } else if ($1 == "Sbjct" && i >= $2 && i <= $4) {
        f[++g] = substr($3, (i - $2) + 1, 1)
        h[c] = h[c] substr($3, (i - $2) + 1, 1)
        }
    }
}

END {
    printf "%s\n", "Seq_name\tbase_position\tquery_base\tsubject_base\tall_subject_base"
    for (j = 1; j <= e; j++) {
        printf "%s\t%s\t%s\n", d[j], f[j], (j % length(a) == 1 ? h[int(j / length(a)) + 1] : "")
    }
}' positions.txt input.txt
Seq_name    base_position   query_base  subject_base    all_subject_base
CD630_08510_19  129 G   T   TA-CT
CD630_08510_19  249 G   A
CD630_08510_19  372 A   -
CD630_08510_19  392 T   C
CD630_08510_19  366 C   T
CD630_08510_20  129 A   G   GA-CT
CD630_08510_20  249 G   A
CD630_08510_20  372 A   -
CD630_08510_20  392 T   C
CD630_08510_20  366 C   T


If you have more positions:
cat positions.txt
129
249
366
372
373
374
375
376
392

awk 'BEGIN {
    OFS = "\t"
}

NR == FNR {
    a[$0]
    next
}

$1 == ">" {
b[++c] = $2
}

$1 != ">" {
for (i in a) {
    if ($1 == "Query" && i >= $2 && i <= $4) {
        d[++e] = b[c] OFS i OFS substr($3, (i - $2) + 1, 1)
    } else if ($1 == "Sbjct" && i >= $2 && i <= $4) {
        f[++g] = substr($3, (i - $2) + 1, 1)
        h[c] = h[c] substr($3, (i - $2) + 1, 1)
        }
    }
}

END {
    printf "%s\n", "Seq_name\tbase_position\tquery_base\tsubject_base\tall_subject_base"
    for (j = 1; j <= e; j++) {
        printf "%s\t%s\t%s\n", d[j], f[j], (j % length(a) == 1 ? h[int(j / length(a)) + 1] : "")
    }
}' positions.txt input.txt
Seq_name    base_position   query_base  subject_base    all_subject_base
CD630_08510_19  129 G   T   TA-GGTACT
CD630_08510_19  249 G   A
CD630_08510_19  372 A   -
CD630_08510_19  373 G   G
CD630_08510_19  374 G   G
CD630_08510_19  375 T   T
CD630_08510_19  376 A   A
CD630_08510_19  392 T   C
CD630_08510_19  366 C   T
CD630_08510_20  129 A   G   GA-GGTACT
CD630_08510_20  249 G   A
CD630_08510_20  372 A   -
CD630_08510_20  373 G   G
CD630_08510_20  374 G   G
CD630_08510_20  375 T   T
CD630_08510_20  376 A   A
CD630_08510_20  392 T   C
CD630_08510_20  366 C   T

NB this won't work if you only have a single position in positions.txt, but if that is a potential use-case you can modify the script to suit, e.g. replace (j % length(a) == 1 ? h[int(j / length(a)) + 1] : "") with f[j].

> **Comments:**
> **username_0**: Jared_mamrot, Thank you so much.  

---

