1
$\begingroup$

I recently prefetched 157 SRA files from an NCBI BioProject using the SRAtoolkit. I then used the toolkit to download those files in FASTA format.

Each individual FASTA file looks something like this:

(example for SRR5678544 file)

>SRR5678544.1 HWI-ST1146_88:5:1112:6472:81473 length=194
CCATGCGGGGTATCGTATGCTTCCTTCTGCACTACCCTTTAGCTGTTCTATATGCTGCCACTCCTCAATTGGATTAGTCTCATCCTTCAATGCTATCACAAGAAGTAGAGACACAAATGCAAGAGGAGCATATAAATTACAAAACACCATCACTGAGGGCCCTAAAGCGGTTCCCACGAAAAAAAGGAGAGTAG   
>SRR5678544.2 HWI-ST1146_88:5:1113:13218:62635 length=194
 TGGATGGTGGTTGGGAGTGGTAAGGTTGAATGAGACACGGTAACGAGTGGGGAGGTAGGGTAATGGAGGGTAAGTTGAGAGACAGGTTCGTCAGGGGACACACCACACCACACCCACACCCACACACCCACACACCACACCCACACACCACACCCACACCCACACCCCCACACACCACACCCCCACCCCACCCA  
>SRR5678544.3 HWI-ST1146_88:5:2206:11224:22269 length=194
 GTCTCTTAACTTACCCTCCATTACCCTACCTCCCGACTCGTTACCCTGTCTCATTCAACCATACCACTCCCGACCACCATCCATCCCTCTACTTACTACAGTATGGTGAGTGGGACATGGTGGATGGTAGGGTAAGTGGCAGTGGAGTTGGATATGGGCAATTGGAGGGTAACGGTTATGGTGGACGGGGGGTG

And so on and so forth. (Ignore the " and ' marks I had to input those to format it right on here)

Question 1: Can someone explain why this file has so many sequences (>)? Every FASTA file I have worked with has one line with ">" at the top and a header sentence followed by a single big paragraph of the genome sequence.

Question 2: If I wanted to combine all the sequences in a single FASTA file to look like a "normal" FASTA file - just a single paragraph - how would I do that?

Reasoning: I'm asking simply because I need to make a BLAST database out of the 157 sequences and the "makeblastdb" command only takes in a single file as the input. I was going to try and combine each individual FASTA file into single paragraphs and then make a huge FASTA file that includes all 157 sequences separated by ">".

$\endgroup$
6
  • 1
    $\begingroup$ This is a fairly straight forward question, generally, someone here should answer this, but your biological objectives are quite important. Firstly, do you have a reference genome, or is that going to be the purpose of the Blast? What are you looking for in the Blast, is it some sort of mapping? $\endgroup$
    – M__
    Commented Sep 6, 2022 at 20:59
  • 1
    $\begingroup$ I've edited the sequences into a coding block, rather than raw text. When I did this the taxa id (Genbank code and data) were on the same line as the sequence data. There was a carriage return this line and the next one starting with >. Please check whether my editing is what was in the original file, i.e. I inserted a carriage return between the taxa id, e.g. >taxa id and the sequence, e.g. GTCTCTTA. Thats quite important $\endgroup$
    – M__
    Commented Sep 6, 2022 at 21:01
  • $\begingroup$ crossposted biostars.org/p/9537532 $\endgroup$
    – Pallie
    Commented Sep 7, 2022 at 9:31
  • 1
    $\begingroup$ The question appears to be a confusion regarding short read data. The point fastqc is normal for assemblies rather than raw fasta is a good point, ie. its not cool to make an assembly directly from the data set in the question. $\endgroup$
    – M__
    Commented Sep 7, 2022 at 15:06
  • $\begingroup$ I was confused about the data but now I understand it's short unassembled reads. My next question is then what now? Do I assemble them? Or do I align them? I have a reference genome I could use but the purpose of the BLAST database is to eventually do some phylogenetic analysis, analyzing closely related species, looking at specific genes in those species, etc.. There's a lot I want to be able to do with the data I am collecting I am just confused on how to process them before those analysis. $\endgroup$
    – rimo
    Commented Sep 7, 2022 at 16:07

1 Answer 1

2
$\begingroup$

Briefly the best strategy I think is to download fastQ files and process this via an assembler. I wouldn't enter raw short read Illumina data into a local Blast database. What I'm pretty certain of is thats a bad move.

Which assembler I personally think is a separate question on Bioinformatics SE. It will depend on the genome, for large eukaryote genomes BWA-mem2 seems popular, but its not my thing and there's other stuff upstream and downstream (someone might dispute that). For assembling small genomes, Spades is popular, but preferences are changing and there's also upstream QC. The other question is whether you have a reference genome to assist assembly.


The update is this a yeast genome with a reference. This is complicated because its a small genome but will comprise heterozygosity. If it was me - i.e. someone way outside the yeast community - I'd use Spades and provide a reference genome to perform an assembly. However, I'm not sure thats at all cool, because you will will reduce the eukaryotic genome to a haplotype. This could get complicated because what you would be left with may not be biologically real, e.g. two alleles at different parts of the same chromosome on different strands may result in a in silico artefact. I think its a separate question personally.

$\endgroup$
4
  • $\begingroup$ I do have a reference genome (S288C: ncbi.nlm.nih.gov/nuccore/NC_001224.1?report=fasta) that I was planning on using. However, I would like to implement de novo assemblies as well just to have that in case it's needed. I have some experience with SPAdes but I assume the tool used doesn't really matter too much? $\endgroup$
    – rimo
    Commented Sep 7, 2022 at 16:49
  • 1
    $\begingroup$ My other question is why would you assemble over aligning? These are all yeast genomes so they are pretty ambiguous in terms of their ploidy/heterozygosity/etc., so I want to make sure I'm doing the right thing for my data. This question is purely curiosity driven and I really appreciate your help so far. $\endgroup$
    – rimo
    Commented Sep 7, 2022 at 16:51
  • $\begingroup$ Hi @rimo on Bioinfo SE this is a separate question. For me its a difficult question because its a eukaryote but a small genome. The issue is (I presume) there'll be heterozygosity in the data and I don't know how the yeast community deal with this (they might ignore it). Most small genomes (my stuff) are haplotypes ( obviously I would ignore it). Most higher eukaryotes use VCF, but this is yeast. It gets complicated because you can't work out the linkage so decomposing a eukaryote into two haplotypes isn't possible. Personally I would assemble via a reference genome. I've edited the reply. $\endgroup$
    – M__
    Commented Sep 7, 2022 at 17:57
  • $\begingroup$ That's exactly what I was thinking. I just posted a separate question to the Bioinfo SE about this exact thing. Thank you so much for all your help with everything else. I'm going to try and assemble via a reference genome to start and then go from there. $\endgroup$
    – rimo
    Commented Sep 7, 2022 at 18:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.