1
$\begingroup$

First of all, sorry if this question has been posted previously (I could not found a solution accessing the previous Q&A).

I have a fasta file as follow:

>LNIV02000036.519060.520603 Pseudomonas aeruginosa
GAACTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAA...

>LDLH01000045.361.1876 Pseudomonas aeruginosa
AGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACACATG...

>LNMQ01000087.5825.7394 Enterococcus faecium
TTTTTATGAGAGTTTGATCCTGGCTCAGGACGAACGCTGGCGGCGTGC...

...

I would like to filter sequences or extract entries for a specific species, let's say Pseudomonas (I have others in mind...). I tried a few things but I think the main problem is my skills associated with bioinformatics. I would like to create a file content Pseudomonas entries (~9k entries in the current database). Any suggestion or tip?

Thank you for your time.

Best,
Fabricio

Improve this question
$\endgroup$

3 Answers 3

Reset to default
4
$\begingroup$

You can do this using seqkit as follows:

seqkit grep -r -n -p '.*Pseudomonas.*' temp.fa

To explain a little, seqkit grep will allow you to search FASTA/Q files by sequence name or sequence itself. In this instance:

  • -r tells that the pattern is a regular expression
  • -n to match by full name instead of just id
  • -p to specify the regular expression pattern to search

There are a bunch of other tools that come with seqkit that can prove indispensable for basic sequence analyses.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Hi @vkkodali, thank you for your time posting it. I did not know about seqkit (living and learning). I'm testing it out and looks like a great tool. Thanks again. $\endgroup$
    – Fabricio
    Feb 19, 2021 at 13:08
  • $\begingroup$ Hi again @vkkodali, another question (sorry, I'm a novice in the world of bioinformatics). Is there a way to retrieve the whole sequence header or ID using seqkit? I filtered the sequences that belong to Pseudomonas and the fasta file contains 38K entries of Pseudomonas. I would like to extract the sequence ID to count how many species of Pseudomonas are in this file (38K entries). illustration: >AB001446.1.1538 Pseudomonas amygdali pv. mori AACTGAAGAGTT... >AB001440.1.1538 Pseudomonas coronafaciens pv. atropurpurea AACTGAAGAGTTT... > ... Thank you for your time and attention. $\endgroup$
    – Fabricio
    Feb 19, 2021 at 13:46
  • $\begingroup$ So you just want the AB001446.1.1538 portion of the header? So many ways to do that... Easiest would be something like grep '^>' in.fa | cut -f1 -d ' ' | sed 's/>//'. Another way is to do seqkit faidx in.fa and that creates a in.fa.fai which is essentially a tab-delimited file with the first column containing the identifiers. $\endgroup$
    – vkkodali
    Feb 19, 2021 at 17:01
  • $\begingroup$ Hi, can be the whole header or portion of it. The portion is the genus and species name (e.g., Pseudomonas putid). I would like to extract the species name from the filtered Psedumonas file (e.g., Pseudomonas amygdali pv. mori, Pseudomonas coronafaciens pv. atropurpurea etc.). I would like to know how many different species of Pseudomonas are there. Thanks again for your time. Fabricio $\endgroup$
    – Fabricio
    Feb 19, 2021 at 20:22
  • 1
    $\begingroup$ This can be done using a combination of grep and sed but I'd be making assumptions about how the organism names in the headers for all your 9k sequences are formatted. Perhaps you should describe the entire issue with an example of what you are starting with and what you would like to see in the end and I can update my answer? $\endgroup$
    – vkkodali
    Feb 19, 2021 at 20:58
2
$\begingroup$

Using Biopython:

from Bio import SeqIO

for record in SeqIO.parse("input.fa", "fasta"):
    if record.description.split()[-2:] == ["Pseudomonas", "aeruginosa"]:
        print(record.format("fasta"))
Improve this answer
$\endgroup$
3
  • $\begingroup$ Hi Chris, right after the post I realize that the issue was not properly described. So, this database has ~9k entries of Pseudomonas species, i.e., hundreds of P. species. The command above asks for the name of each Pseudomonas entry, correct ("Pseudomonas", "X")? Is there a way to retrieve the Pseudomonas entries without specifying all Pseudomonas names (I do not have a list name)? Thanks. $\endgroup$
    – Fabricio
    Feb 17, 2021 at 19:47
  • $\begingroup$ @Fabricio Hi Fabricio, do you understand the above code? If so you can easily modify the if line to search for entries with "Pseudomonas" $\endgroup$
    – Chris_Rands
    Feb 17, 2021 at 19:50
  • $\begingroup$ @Fabricio Add this at the top of the code, after the other import statement: import re . Change the if expression to this: if re.findall(r'pseudomonas', record.description, re.I): , in order to search in the fasta header for the word pseudomonas, in case-insensitive manner. $\endgroup$
    – Timur Shtatland
    Feb 17, 2021 at 21:19
0
$\begingroup$

"describing the entire issue".

I have a database (silva 16S formatted by DADA2 - fasta format) that contains ~300k entries for bacteria (16S rRNA gene). e.g.,

AB001441.1.1538 Pseudomonas syringae pv. broussonetiae AACTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGG... KR054096.1.1224 Marinococcus halophilus GAAGGTAAAGGCTTCCCCAAGGGCGACGATGGGTAGCCGACTTGA... AB681405.1.1478 Paenibacillus barengoltzii GACGAACGCTGGCGGCGTGCCTAATACATGCAAGTCGAGCGGAGT... ...

I would like to retrieve sequences for some genus, e.g., Pseudomonas. After running the seqkit command (mentioned by you above) I ended up with ~9 of Pseudomonas entries in a fasta file (example below).

AB001441.1.1538 Pseudomonas syringae pv. broussonetiae AACTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACA... AB001442.1.1538 Pseudomonas amygdali pv. eriobotryae AACTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCCTAACA... ...

This pool of 9k sequences of Pseudomonas has several species of Pseudomonas (e.g., P. putida, P. syringe Et Cetera). Since some species of Pseudomonas have multiple sequences in the database (e.g., P. putida has 1.8k entries) I would like to know how many different species of Pseudomonas are there (in the 9k Pseudomonas entries), by names (i.e., in the filtered Psedumonas fasta file). At the end of the day, I would like to say that in this database I have 9k entries of Pseudomonas which represents 'X' number of Pseudomonas species.

What I'm thinking... Retrieve the full ID (AB001442.1.1538 Pseudomonas amygdali pv. eriobotryaeall) or if possible the portion of the ID (headers) that contains the species names (e.g. Psudomonas amygdali pv. eriobotryaeall) and with this, I can list the species names and count them. Does this explanation help to understand what I'm trying to do?

Let me know if you need more info.

Thank you for your time. Fabricio

Improve this answer
$\endgroup$

Your Answer

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

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