extract multiple sequences from a fasta file with a list of numbers in header

Question

[barvazduck]

Hi all,

I got a large fasta file that I need to extract multiple sequences from the header for each sequence is composed of several parameters such as Chromosome, Genetic distance etc, and also has an ID number in the end. I need to extract sequences based on the ID number.

I tried using 'grep':

grep -A 2 -wFf LIST.txt IN.fa > OUT.fa

But this matches also non-specific numbers. For example if I search for ID: 79695 I also get sequences with IDs such as 1379695 and 7969522.

Any idea how to solve this or other solutions?

Thanks

[PoGibas]

Can you post example of your header? Is it ID: 79695 or ID:79695 or just 79695?

[barvazduck]

>chr7B, genetic location: 72.6 cM, bin num: 3629, phase id: 1, scaf14, marker id: 1079695

[Devon Ryan]

This is probably simplest in either biopython or bioperl. They'll handle the presence of multi-line entries transparently and allow you to just use a regex to get out the ID (just store the target IDs in a hash/dict and query that).

Answer 1

[robert.davey]

As long as I'm understanding your question accurately, the simplest way would be to combine grep and sed. If your LIST.txt contains IDs, and you are specifically looking for "marker id: "+<exact ID match>, then the following should work:

grep -A 1 -wFf <( sed -r 's/^/marker id: /' LIST.txt ) IN.fa > OUT.fa

Edit: apologies, missed the -w flag off in my retype. The previous answer failed when the start of the ID field contained a search ID.

[barvazduck]

Thanks,

But this still wont get the exact ID number but also all the nested ones

[robert.davey]

EDIT: Fixed missing -w flag

Seems to work properly for me:

test.fa:

>chr7B, genetic location: 72.6 cM, bin num: 3629, phase id: 1, scaf14, marker id: 1079695
ATGCTAGCTAGCTAGCCGGCTAGCTACTATCGGCTATCGTACGTAGC
>chr7B, genetic location: 72.6 cM, bin num: 7969, phase id: 1, scaf14, marker id: 1079696
GTATCTGGCATCTTACTGACGGCGATCGATGCGCGCTAGCTAGCTAT

When LIST.txt:

davey:~/temp$ cat LIST.txt
7969

davey:~/temp$ grep -A 1 -wFf <( sed -r 's/^/marker id: /' LIST.txt ) test.fa
davey:~/temp$

I get no output. With LIST.txt:

davey:~/temp$ cat LIST.txt
1079696

I get:

davey:~/temp$ grep -A 1 -wFf <( sed -r 's/^/marker id: /' LIST.txt ) test.fa
>chr7B, genetic location: 72.6 cM, bin num: 3629, phase id: 1, scaf14, marker id: 1079696
GTATCTGGCATCTTACTGACGGCGATCGATGCGCGCTAGCTAGCTAT

What bash/grep versions are you using?

[barvazduck]

bash version 4.3.11(1)-release

grep 2.16

Anyway, it works partially because if you use 107969 then you get

>chr7B, genetic location: 72.6 cM, bin num: 3629, phase id: 1, scaf14, marker id: 1079695
ATGCTAGCTAGCTAGCCGGCTAGCTACTATCGGCTATCGTACGTAGC
>chr7B, genetic location: 72.6 cM, bin num: 7969, phase id: 1, scaf14, marker id: 1079696
GTATCTGGCATCTTACTGACGGCGATCGATGCGCGCTAGCTAGCTAT

Thanks for the efforts.

[robert.davey]

I've fixed my previous post. Sorry about the omission of the -w option.

[barvazduck]

Yeah I tried that, seems like something is wrong with my LIST.txt

[robert.davey]

Do you have non-standard (i.e. Windows) carriage returns in your file? You can use dos2unix to get rid of them.

It does indeed fail with the following:

1079696^M
1079695^M

If you open up your LIST.txt in vi, the status bar will help tell you if you have a DOS file or not: "testlist.txt" [dos] 2L, 18C

[barvazduck]

Actually it was a very stupid mistake, I just didn't use -w when I checked the result file using another grep line

Answer 2

[Pierre Lindenbaum]

Sort your list of ID LIST.txt

Linearize the fasta, extract the id after marker_id, sort on the first column, join with the LIST, convert back to fasta.

awk '/^>/ {i=index($0,"marker id:"); printf("%s%s\t%s\t",(N==0?"":"\n"),substr($0,i+11),$0);++N;next;} {printf("%s",$0);} END{printf("\n");}' input.fa|\
sort -t '    '  -k1,1 |\
join -t '    ' -1 1 -2 1 LIST.txt - |\
cut -f 2,3 |\
tr "\t" "\n"

Answer 3

[Bara'a]

If you are aiming at retrieving certain sequences by ID for further manipulation, you can use the Bio.SeqIO.index() function of Bio.SeqIO module provided in Biopython .

This functions allows you to index moderately large sequence files without consuming much memory as it stores where each record is within the file and parse it on demand given identifiers .

A short Biopython code will do the work :

from Bio import SeqIO

def get_ID (identifier):
    parts= identifier.split (",")

    assert len (parts)==6 and parts[0]=="chr" and parts[1]=="genetic location:" and parts[2]=="bin num:" and parts[3]=="phase id:" and parts[4]=="scaf" and parts[5]=="marker id:"

    return parts [5]

file_dictionary=SeqIO.index(" your_file.fasta ", " fasta ", key_function=get_ID)
file_dictionary.keys()

handle=open(" selected.fasta" , "w")

for ID in [ "id1", "id2", "id3", ...  ] :
      handle.write( file_dictionary(ID).seq )

handle.close()

You can refer to Biopython tutorial section 5.4.2 for more details and: http://biopython.org/DIST/docs/tutorial/Tutorial.html

Also, you can find other indexing methods for extremely large sequence files in section 5.4.

Hope you find this useful :)

EDITED to suit your file's header format.

[robert.davey]

I'd hate to write out the IDs by hand in your for loop :) Consider reading them in by file:

with open('LIST.txt','r') as ids:
  for ID in ids:
    ID = ID.rstrip()
    handle.write( file_dictionary(ID).seq )

[Bara'a]

I wrote that command assuming that she wants to retrieve very few sequences, so why bother reading them from a file?!

Anyways, All roads lead to Rome !! :D

[robert.davey]

Quite so! :)