Extract sequences from a list of ID

Question

[doinelpierrot]

Hello everybody,

Here i my issue, I have a Trinity.fasta file like this

>TRINITY_DN5631_c0_g2_i1 len=947 path=[0:0-946]
TACAACTTGAACATCAACAATGGTTGCGCAGCTATTGCCATCCGCGACGTTCGAGGACTGCGTGCGAA
>TRINITY_DN62279_c1_g1_i1 len=298 path=[0:0-297]
TATTACCATTATTATTATTATCATATTTATGTTCATTATTATCATTATCATAATCATTATCATCTTGATA
 ...

And I also have a list of id in a id.txt file :

TRINITY_DN16359_c0_g1_i4
TRINITY_DN62279_c1_g1_i1
...

I am trying to extract for my fasta file the sequences that have their id in my txt.file. I am using seqkit for that, but with no success :

 seqkit grep -n -f id.txt Trinity.fasta -o result.fa

Does anyone know how to fix it ?

[GenoMax]

This has been asked multiple times. Here is one thread: Extract fasta sequnces not matching a list Answers include extracting both matching and non-matching sequences. If you want to keep using seqkit try: C: Extract specific sequence from FASTA file

[doinelpierrot]

Yes I have checked biostars first that's how I found the seqkit command I am trying to use. My goal is not to remove sequence from my fasta but to have a file with all the sequences that are in my list

[GenoMax]

Did you check the links I included above? Second link has seqkit command, if you want to stay with that program. -nrif option mentioned there should cover your use case, I would think.

[doinelpierrot]

seqkit grep -nrif ids.txt test.fa totally did the job thank you very much @genomax !

[doinelpierrot]

finally it does not work so great since as shenwei says, " it may produce some unwanted results. For example, seq_1 matches seq_10 with -nri"

[Fatima]

If your sequences are only one line you can use this command:

cat id.txt | while read line ; do grep -A 1 ${line} file.fasta >> outputfile.fasta ; done

If your sequences are multi-lines you can convert them to a one-liner fasta file first and then use the above command:

awk '/^>/ {printf("\n%s\n",$0);next; } { printf("%s",$0);}  END {printf("\n");}' <  Trinity.fasta >> file.fasta

[doinelpierrot]

Thank for your answer I have tried this but there is 6835 sequences in outputfile.fasta whereas I have a liste of 6750 ID. Is it possible to fix that ?

[Fatima]

I usually use this for accession numbers that are unique, and ${line} works well for them, but maybe you can try "${line} ", because this one requires a space after ID, and 10 won't be picked when looking for 1.

Alternatively, you can use grep --perl-regex "${line}\s", but I haven't tried mixing this with -A1 or -A 1 option. You can try!

[shenwei356]

Option -w should be added to grep for exact matching.

[Neel]

Hi, i want to extract seq from list of id but Seqkit gives this below error and empty output file generated

[WARN] symbol ">" detected, it should not be a part of the sequence ID/name: >CP077971.1_4943 P_aeruginosa_ZPPH33resistance-nodulation-cell division (RND) antibiotic efflux pumpType B NfxB

Thank you

[GenoMax]

You need to remove > from ID's as the warning tells you to do.

[Neel]

Thank you....It has been done.

Answer 1

[shenwei356]

Referring to the usage, you should not switch on -n , just use seqkit grep -f id.txt seqs.fa.

-n, --by-name               match by full name instead of just id
-i, --ignore-case           ignore case
-r, --use-regexp            patterns are regular expression

For -nrif, it's partly matching full FASTA header by regular expression and case-ignored, with patterns from file. This may produce some unwanted results. For example, seq_1 matches seq_10 with -nri.

[doinelpierrot]

Hmm unwanted results, it's exactly what I got since I ended up with 20 mysterious sequences more than expected (out of 7000) wit -nrif.

I have tried seqkit grep -f id.txt Trinity.fasta > contaminants.fasta but unsuccesfully ...

[shenwei356]

Hmm... I've tested with your sample data, everything is OK.

You may check the input data (check extra whitespace or tab at front/end of IDs, file encoding, even manually copy and search in the file with text editor.), and try with additional options, e.g., switch on -i to ignore cases.

[doinelpierrot]

My bad it did the job, thanks again !

[Ishak]

thanks

seqkit grep -r -f id.txt Cucumis_melo.Melonv4.cds.all.fa -o result.fa

it works!

Answer 2

[Renesh]

You can also try Python package bioinfokit (v1.0.2 or later) for extracting the sequences from the FASTA file (check how to install https://reneshbedre.github.io/blog/howtoinstall.html#how-to-install)

from bioinfokit.analys import fasta

fasta.extract_seq(file='Trinity.fasta', id='id.txt')

Check detailed usage at https://reneshbedre.github.io/blog/howtoinstall.html#extract-sequences-from-fasta-file

Answer 3

[hxj5]

Could also use seqtk or AWK from an old thread Tool: Retrieve a subset of FASTA from large Illumina multi-FASTA file:
seqtk subseq large.fa id.txt
or
awk 'BEGIN{while((getline<"id.txt")>0)l[">"$1]=1}/^>/{f=l[$1]}f' large.fa

Answer 4

[Md]

you can use samtools command :

samtools faidx Trinity.fasta TRINITY_DN16359_c0_g1_i4 TRINITY_DN62279_c1_g1_i1