Since you're not working with a model organism, there will be many reads that do not have a BLAST hit but still are probably derived from your bivalve. I do not suggest using BLAST to filter your reads for that reason.

I would try a few things:

• Assemble your transcriptomes then map the transcripts to the draft genome.
• Map the reads to the genome, then assemble the transcriptomes.
• blast the two transcriptomes against each other and see if any of the transcripts were not hit. You can investigate the reasons why afterward.

Most importantly for non-model organisms and highly heterozygous species - do not assemble a transcriptome using data from more than one individual. The number of allelic differences in highly heterozygous species will make the transcript graphs a jumbled mess and you'll end up with a bunch of partial and erroneous transcripts. You can merge transcriptomes from multiple individuals later with tools like orthoFinder.

First, to answer your question about mapping to a low-quality reference:

1. Mapping

2. Reasons to not BLAST your raw reads

There are a few reasons that I do not suggest starting with BLAST, but other options may better suit your needs.

• BLASTing 2 billion reads will take several weeks to months even on a machine with many cores (90+).
• Since you're not working with a model organism, there will be many reads that do not have a BLAST hit but still are probably derived from your bivalve. I do not suggest using BLAST to filter your reads for that reason.
• Similarly, many contaminating reads will not hit to anything in NCBI databases, so you will miss many of the contaminants as well. Leaving them in for now and removing with the strategy I describe in #4, then reassembling, may work better.

3. Using multiple individuals of highly heterozygous species in assembly.

Most importantly for non-model organisms and highly heterozygous species - do not assemble a transcriptome using data from more than one individual. The number of allelic differences in highly heterozygous species will make the transcript graphs a jumbled mess and you'll end up with a bunch of partial and erroneous transcripts. You can merge transcriptomes from multiple individuals later with tools like orthoFinder.

4. A strategy to filter transcripts based on mapping and BLAST

Given #1, #2, and #3 above, I propose the following strategy to filter your reads that may prevent removing true bivalve reads and will guarantee removal of reads from known contaminating species.

I would try this:

• Assemble the transcriptome of each individual separately.
  • This avoids misassemblies due to heterozygous sites and over-filtering by Trinity or other assemblers.
• Concatenate all of the transcriptomes into one file. From this point on we will use the information from all individuals to help filter.
  • Use somthing like cat txome_individual1.fasta txome_individual2.fasta > all_individuals.fasta.
• Map all the transcripts to the reference genome s using a spliced aligner like STAR or HISAT2.
• Get a bam file of all the transcripts that didn't map to the reference genome. These may be from contaminants.
  • Use samtools view -f 4 -b mapped_transcripts.bam | samtools sort - > output.sorted.bam to only keep the transcripts that did not mapped to the reference. -f 4 is the magic flag here that keeps things that didn't map.
• Make a list of all the transcripts that did not map to the reference...
  • samtools view output.sorted.bam | cut -f1 | sort | uniq > unmapped_transcripts.txt
• ...and get only the transcripts from your fasta file that are potential contaminants.
  • seqtk subseq all_transcripts.fa unmapped_transcripts.txt > unmapped.fasta
• BLAST all of the unmapped.fasta and identify which ones are from contaminating species. blobtools works well for this if there are too many to keep track of.
• Use seqtk again to only get the known contaminating seqs from BLAST or blobtools. Let's say you call this known_contaminants.fasta
• Map all reads from all individuals to the known_contaminants.fasta.
• Get all the reads that map to known contaminating transcripts, known_contaminants.fasta using:
  • samtools view -F 4 all_reads_to_contams.bam | cut -f1 | sort | uniq > contaminating_reads.txt
• Now remove all of the contaminating reads from your individual fastq files using bbmap filterbyname.sh.
• Reassemble the filtered fastq files.

This strategy has the following benefits:

• Does not require weeks/months of compute time to BLAST 2 billion reads.
• Ensures that all reads from the reference genome are left untouched.
• Ensures that all reads from known contaminants are removed.
• Leaves in reads with no hits to BLAST, which could be from unassembled parts of the reference genome.
• Avoids potentially filtering out low-abundance contaminants by not inputting all 2 billion reads to the transcriptome assembler, but assembling on a per-individual basis.
• Leverages contaminants from different individuals to filter reads from all individuals.