makeblastdb reports ambiguous nucleotides (how do I debug a fasta file?)
1
0
Entering edit mode
9.5 years ago
arronar ▴ 280

Hello. I had a fastq file from an illumina machine and convert it in fasta. Then i used dos2unix command to fix any problem with new lines.

dos2unix merged_zea.fasta

after that I used the makeblastdb command to create database to use it with tblastx.

makeblastdb -dbtype nucl -in merged_zea.fasta -input_type fasta -out zeaDB -max_file_sz 2GB

The problem is that makeblastdb returns me this error:

Error: (1431.1) FASTA-Reader: Warning: FASTA-Reader: First data line in seq is about 100% ambiguous nucleotides (shouldn't be over 40%)

I searched int net and everyone is talking about 2 possibilities.

  1. Some mistake in the fasta file
  2. There are enough NNNN in sequencies.

My sequencies has no NNNN and here is how they look like in fasta file

>M02381-6-000000000-AA1W0-1-1101-15670-1350 1-N-0-1
GAAGCAGTGTTGACGTAGTTCTGAGCCATGGCCATCATGTACTGAACATCGAGGTTAGCTTCAACACAAGAGTTGGGGTTGAAAGAGCAAGCACTGGGGTTGTTGGGTCCGATGACCTTGTCAACCTTGTCGGAAGCGGTTCCCATAGCAG
>M02381-6-000000000-AA1W0-1-1101-15485-1351 1-N-0-1
CCAGAAGTCAAGTCGAACAAAACGCCAATTGAGGATGGCGACGATCCCGACCCGCCAAGTAAGAATGCCACGCCAGAGGAGAAAGCTTCGGCCATCCGCGCAATACAAGACAACTCGGATGAGATGTTCGACGTTATTG
>M02381-6-000000000-AA1W0-1-1101-15538-1354 1-N-0-1
GGGAAGGACTGACTGCTATGAACGCGCTGTTTTATTCACGGTAGCGTACCGGCGCGCACCATTGCTTCACGTGTTGGGTTAGTGAAGAGGGCTCATCGAACGACCCAATTAAGGGCCCCTCTC
>M02381-6-000000000-AA1W0-1-1101-16188-1354 1-N-0-1
GGACAGTTGGGGGTGCTAGTATTTAACGGCCAGAGGTGAAATTCTTGGATTCGTTAAAGACTAACTAATGCGAAAGCATTCACCAAGGATGTCTTCTTTAATCAAGAACGAAAGTTGGGGGATCGAAGAGGATCAGATACCCTCGTAGT
>M02381-6-000000000-AA1W0-1-1101-15230-1355 1-N-0-1
GTGCACACACACGAATGCTCAACACCGCCCGGCCCTATTTGGGTTGAGGTAGAGTGTGCCTGTTTGGACCCGAAAGGTGGTGAACTATGCCTGAACAAGGTGAAGCCAGGGGAAACCCCGGTGGAGGCCTGTAG
>M02381-6-000000000-AA1W0-1-1101-16511-1355 1-N-0-1
ATACTGGAGTTCCCACAGCAATAGATAGATCATGGTCATAGTCATCATAATTATCTTTGCTTCCTGACCTGTTCCCAGAGGAGTCGTAGTC
>M02381-6-000000000-AA1W0-1-1101-14969-1356 1-N-0-1
GGTTTAAAGGGTCCGTAGGCGGTTTTATAAGTCAGTGGTGAAAGTTTGCGGCTTAACCGTAAAGTTGCCATTGATACTGTAGAACTTGAATAATTGTGAAGTGGTTAGAATAAGTAGTGTAGCGGTGAAATGCATAG
>M02381-6-000000000-AA1W0-1-1101-15132-1356 1-N-0-1
CAATTCCGTTGAATTCGAGGAAGGGATAATCCAAGCAATCAGACATTATTTCTTTTTGTTTGAGGTGTCCAAGGCTCCATTCACACTCCACCAGTTATTCCCGTTTTTGAAGGGTACCGACGTTCCCGATCAAGAGATACCCTCCACACT
>M02381-6-000000000-AA1W0-1-1101-15192-1356 1-N-0-1
CGTGAATCGTCTGGCGGCAGACGGCGTTGCCCTCGCGGGCAATGCGCACGATCGCGTAGTTGGCGGTGCAGCCGTCGGCCTCGCCCTGCTCGGTGGAGTTCACCTTCGGGTCCGAAAGGAGGGCGTCGAGGTCCGGCGAGTTGCACTCGG
>M02381-6-000000000-AA1W0-1-1101-16071-1357 1-N-0-1
GCTATGTATGTTGCTATTCAAGCTGTCCTTTCTCTTTATGCATCTGGACGTACAACTGGTATTGTTTTAGATACTGGTGATGGTGTTTCTCACACAGTCCCAATTTATGAAGGATATGCAC
>M02381-6-000000000-AA1W0-1-1101-15463-1357 1-N-0-1
GCTGTGAGCCATGTTGCGGTCAACCTTGCAACACTAGTGACGGATTATATTGCTGTCTTTCTTTCTTCCCTGGATGTCTTTGTGCTTACCCACAATTCTATGCTTCCACTCAAGAACAACAATGTGCTTGGGTCAACCACTGTGTCCCAA
>M02381-6-000000000-AA1W0-1-1101-16160-1357 1-N-0-1
TGTTTGAATTTCTTCACTTTGACATTCAGAGCACTGGGCAGAAATCACATTGTGTTAAAATCTTTTCATGACCATCACAATGCTTTGTTTTAATTAAACAGTCGGATTCCCCTGGTCAGTAACAGTTCTAAATTAGCTGTTCATTGTATA
>M02381-6-000000000-AA1W0-1-1101-15586-1357 1-N-0-1
GGATCATACTCTACCTCGCGCAAACTAGAGATAGGCACAAGACTATGGAAGTTTCGTTGGAAAGCGGCGCCAGAGAGGGTGAAAGCCCCGTGGATTAGATTTGTGTAGCGTGTGAGTTGGGGGTGGCCCCGAGCGAGTCGTGTTGTTTGGG
>M02381-6-000000000-AA1W0-1-1101-16490-1359 1-N-0-1
TTGCATACCCGATTGCCTTCTTAATGATCTCAACCCAGTTGGTCCTCTCCTCCTCTTTGAGAAGATAGAATGTTCTGGTCTTCGAAGGGAAGATCAGGGTGAACGGATAGAGAGTGTTGGAC
>M02381-6-000000000-AA1W0-1-1101-15381-1359 1-N-0-1
AGATACCTTTATGCCGGCTGCTGACTACTTATGCAGGATTGGCATTATGACTGCTTTGACTAGAAAACTCTGGAAAAGAAAAAAAATTGTAAGGCTTGGGCGCTTAGTGTCTCTT

So what produces this error? Can somebody give me a hint?

Thanks.

sequence next-gen • 2.9k views
ADD COMMENT
0
Entering edit mode

answer number 1 above is the very likely answer.

that your file is incorrect even though it may look correct to you - dos2unix won't fix all line ending problems - take a single line and see what happens. The file that you list above works fine on my system.

ADD REPLY
0
Entering edit mode

I am taking a single line and it works. Maybe there is a specific line produces this error. But how can I "debug" a fasta file?

ADD REPLY
1
Entering edit mode
9.5 years ago

do a binary search - take upper half and lower half of your fasta file, then repeat on those

actually do also a:

egrep "N+" --color=always yourfile | head

just to make sure you don't actually have sequences with lots of Ns

PS: I changed the title to add how to "debug" a fasta file - that is a funny and appropriate way to put it

ADD COMMENT
0
Entering edit mode

I run the egrep command without the | head, and actually it seems that there are some NNNNs in some sequencies

>M02381-6-000000000-AA1W0-1-1102-17197-13281 1-N-0-1
>M02381-6-000000000-AA1W0-1-1102-29267-13281 1-N-0-1
>M02381-6-000000000-AA1W0-1-1102-6311-13281 1-N-0-1
>M02381-6-000000000-AA1W0-1-1102-26477-13281 1-N-0-1
>M02381-6-000000000-AA1W0-1-1102-10131-13281 1-N-0-1
>M02381-6-000000000-AA1W0-1-1102-28945-13281 1-N-0-1
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
>M02381-6-000000000-AA1W0-1-1102-26404-13281 1-N-0-1
>M02381-6-000000000-AA1W0-1-1102-24512-13281 1-N-0-1
>M02381-6-000000000-AA1W0-1-1102-21563-13281 1-N-0-1
>M02381-6-000000000-AA1W0-1-1102-20297-13281 1-N-0-1
>M02381-6-000000000-AA1W0-1-1102-13548-13281 1-N-0-1
>M02381-6-000000000-AA1W0-1-1102-3763-13281 1-N-0-1
>M02381-6-000000000-AA1W0-1-1102-9024-13281 1-N-0-1
>M02381-6-000000000-AA1W0-1-1102-10869-13281 1-N-0-1

I've made a little mistake in a command before, by telling that there are no NNNNs in sequences.

So now what? Should I replace them with blanks?

ADD REPLY
1
Entering edit mode

blank won't really do better, you need to remove those sequences.

If all sequences are on one line you could do something like this to remove all sequences with more than 5 consecutive Ns:

egrep -B 1 "N{5,}" yourfile | grep -v '^--'  > newfile.fa

if the sequences are on multiple lines then you would need to use an approach that properly parses your sequences, say bioawk

cat yourfile | bioawk -c fastx ' $seq ~ /NNN/ { print ">" $name "\n" $seq } '
ADD REPLY
0
Entering edit mode

Are you sure that the command with bioawk works properly? I run it and nothing changed in the file.

ADD REPLY
1
Entering edit mode

actually never just run stuff you see on the internet - I just typed that in on a whim who knows if it is right - I think it is at least 90% there but I do have many other things on my mind as well

looks like I am selecting the pattern instead of opposing it with !~ - read a bit on how awk works then build your own little script and test every step of it. It should work but simple things can be tricky as well.

ADD REPLY

Login before adding your answer.

Traffic: 3256 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6