3
$\begingroup$

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.

Improve this question
$\endgroup$
6
  • $\begingroup$ 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. $\endgroup$
    – Community Bot
    Commented Mar 23, 2023 at 14:42
  • 3
    $\begingroup$ Have you tried looking for ways to get alternative output formats from blast that would be easier to parse? $\endgroup$
    – Maximilian Press
    Commented Mar 23, 2023 at 16:21
  • 1
    $\begingroup$ I agree with @MaximilianPress given the approach used. The thing about Biopython is this already done - so it's easy $\endgroup$
    – M__ ♦
    Commented Mar 23, 2023 at 18:33
  • 1
    $\begingroup$ 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. $\endgroup$
    – gringer ♦
    Commented Mar 23, 2023 at 19:48
  • $\begingroup$ Will you always have just one alignment in your file or can you have muiltiple? Does our solution need to deal with both cases? $\endgroup$
    – terdon ♦
    Commented Mar 24, 2023 at 12:33

1 Answer 1

Reset to default
2
$\begingroup$

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

Improve this answer
$\endgroup$
1
  • $\begingroup$ Jared_mamrot, Thank you so much. $\endgroup$
    – user14714429
    Commented Mar 31, 2023 at 16:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.