Hi All,
I have:
A) a fasta file containing a few thousand short sequences (40bp each)
B) an 8Gb fastq.gz file containing millions of unassembled illumina paired end reads (100bp each)
I would like to blast each of the sequences in A against those in B and count the hits of each.
Does anyone know a good strategy for this? To use blastn I would need to convert the fastq.gz to fasta first in order to create a blast database, right? Is there any better way to do this, and any tools that might be computationally fast?
hope someone can advise!
cheers
David
You can create a blast database of your 8Gb fastq.gz file ensuring first they are assembled before creating the database. Then you can run blast of your short sequences (A) based on your created database (B). You can do that using CLC software.
You go for standalone blast or try blaststation for similarity search of millions of read against few thousand query sequence
You can create a blast database of your 8Gb fastq.gz file ensuring first they are assembled before creating the database. Then you can run blast of your short sequences (A) based on your created database (B). You can do that using CLC software.
As you can see in my answer, i mentioned that data B (8Gb fastq.gz) should be assembled first. Don't know what point you are trying to make?. #YiXu
One question - are the short reads from a single organism or metagenomic sample?
If single organism, then you could possibly approach it the other way around - align the FASTQ file to the FASTA file using bowtie2 in 'local' mode, rather than BLAST. Although you are doing it the other way round, it might hopefully answer your question.
Alternatively, if the FASTQ file contains sequences from several different organisms, then you could try one of the metagenomics pipelines such as qiime/MG-RAST first, to identify which organism each read came from.
Thanks for the input from everyone.
@Jonathan David Moore, the fastq reads are from one organism.
@Lishuang Shen, Your suggestion (also Jonathan's) to align the FASTQ reads to the FASTA file is really interesting.
A bit more detail about what i am trying to figure out: The 8GB of FASTQ reads are a whole genome (a plant genome) that is unassembled. I would like to know how many times each of the 40bp FASTA sequences map to that genome but would prefer not to have to assemble it!
Combination of "grep" and "REGEX" may help since blast excludes low complexity region even you turn the filter off.
I think following could be the solution for your problem
1. Merge/Concatenate
cat file_1.fq 2.fq > lib.fq
2. Conversion of fastq to fasta files
seqtk seq -a lib.fq > lib.fa
3. Making Database
formatdb -i lib.fa -p F
4. Blast
nohup blastall -p blastn -i ref.fa -d lib.fa -e 0.00001 -b 100000000 -v 100000000 -o out.txt &
5. To show reads
grep '^E00602:' out.txt |cut -f1 >read_ID_list
6. To extract final fasta file
nohup read_ID_list lib.fa 1 >read.fa &
then use seqman pro for contig assembly.
- Badges
- Science topic
- Similar topics
- Bioinformatics
- Bioinformatics and Computational Biology