Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-03-31T16:45:58.123Z
Answers.Text: Jared_mamrot, Thank you so much.
Answers.UserDisplayName: user14714429
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2023-03-27T22:28:51.427Z
Answers.LastActivityDate: 2023-03-27T22:28:51.427Z
Answers.UserDisplayName: jared_mamrot
Body: <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>
  
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-03-23T14:42:57.007Z
Comments.Text: 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.
Comments.UserDisplayName: Community
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-03-23T16:21:58.813Z
Comments.Text: 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?
Comments.UserDisplayName: Maximilian Press
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-03-23T18:33:34.597Z
Comments.Text: I agree with @MaximilianPress given the approach used. The thing about Biopython is this already done - so it's easy
Comments.UserDisplayName: M__
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-03-23T19:48:43.740Z
Comments.Text: 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.
Comments.UserDisplayName: gringer
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-03-24T12:33:41.323Z
Comments.Text: Will you always have just one alignment in your file or can you have muiltiple? Does our solution need to deal with both cases?
Comments.UserDisplayName: terdon
Comments.ContentLicense: CC BY-SA 4.0
Comments.CreationDate: 2023-03-24T12:56:22.727Z
Comments.Text: 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?
Comments.UserDisplayName: terdon
Tags: <blast><linux><shell>
Title: Extract specific nucleotide base position for query and subject from text file
ThreadText: # 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 @username_2 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.  
  
  ---
  
  