Extracting subsequence from FASTA using Samtools faidx not working

Question

[hmg]

Hello,

I have a FASTA file, and I want to extract a subsequence from it.

>first
AGCATGCATCGTACGATCGTACGATCATCATC
GTACGATCGATCGATCGATCGATCGTGTGCTA
GACGATCGATCGATCGA
>second
GTACGATGCATCGATCGATCGATCGATCGATC
GTGATCGTGTCGTCGTCGTCGATAAGCATCGT
AGTACGATCGATCGTA

I know I can do this with samtools:

samtools faidx seq.fa first:1-5

output is:

>first:1-5
AGCAT

However, my actual sequence names are much more complicated, and anything other than one word doesn't seem to work.

If seq2.fa is:

>first seq
AGCATGCATCGTACGATCGTACGATCATCATC
GTACGATCGATCGATCGATCGATCGTGTGCTA
GACGATCGATCGATCGA
>second seq 
GTACGATGCATCGATCGATCGATCGATCGATC
GTGATCGTGTCGTCGTCGTCGATAAGCATCGT
AGTACGATCGATCGTA

Then the output of:

samtools faidx seq2.fa first:1-5

>first:1-5
seqAG

​```​ samtools faidx seq2.fa first seq:1-5

first seqAGCATGCATCGTACGATCGTACGATCATCATCGTACGATCGATCGATCGATCGATCG TGTGCTAGACGATCGATCGAT seq:1-5 ```

samtools faidx seq2.fa 'first seq':1-5

>first seq:1-5

Also, if my sequence description contains a number (like seq3.fa):

>seq123
AGCATGCATCGTACGATCGTACGATCATCATC
GTACGATCGATCGATCGATCGATCGTGTGCTA
GACGATCGATCGATCGA
>seq456
GTACGATGCATCGATCGATCGATCGATCGATC
GTGATCGTGTCGTCGTCGTCGATAAGCATCGT
AGTACGATCGATCGTA

samtools faidx seq3.fa seq123:1-5

>seq123:1-5

Is there a way to use samtools to extract a subsequence when the descriptor is more complicated than just a word? I've tried with samtools 0.1.18 and 0.1.19.

Thanks

Answer 1

[Pierre Lindenbaum]

faidx stops reading the name after the 1st blank: https://github.com/lh3/samtools/blob/master/faidx.c#L89

while ((ret = razf_read(rz, &c, 1)) != 0 && !isspace(c)) {

you could fix this by converting the spaces to an underline

tr " " "_" < in.fa > out.fa && samtools faidx out.fa

[hmg]

I tried that, but no luck:

>first_seq
AGCATGCATCGTACGATCGTACGATCATCATC
GTACGATCGATCGATCGATCGATCGTGTGCTA
GACGATCGATCGATCGA
>second_seq 
GTACGATGCATCGATCGATCGATCGATCGATC
GTGATCGTGTCGTCGTCGTCGATAAGCATCGT
AGTACGATCGATCGTA

samtools faidx seq4.fa first_seq:1-5

>first_seq:1-5

[Pierre Lindenbaum]

Please show me the output of grep first_seq seq4.*

[hmg]

Just:

​>first_seq

[Pierre Lindenbaum]

If you really ran grep first_seq seq4.*, then you didn't correctly indexed the fasta sequence.

[SES]

I think the blank lines are the issue, and samtools should have warned when the index creation was attempted. The output of samtools faidx seq4.fa may be informative.

[hmg]

Yep--I didn't index it at all. Silly error by me. Thanks for the help.