Hello everyone,
To explain, I have a SAG (single amplified genome) from a bacterium of interest, I have assembly reads into contigs (with SPAdes) and then I have extract the 16s RNA sequence.
I have also several metagenomics datas (filtered reads, ~100pb, ~200 000 000 reads) from differents lakes and the objectif is to calculate the relative abundance of my bacterium in these metagenomics datas.
For that, I thought about aligning the 16s RNA against the metagenomics reads with this simplified script:
for meta in metagenome: do blastn -evalue 1e-7 -db metadb -query ARN_16s -perc_identity 97 -outfmt 6 -out resultats/$meta/resultats.blasttab (I keep the hits only if the alignment is greater than 40pb)
The problem is that I have only few hits in the first lake, but no hits in the four others .. (I'm supposed to have results in the other lakes)
Did you thinks the parameters that I utilise are good ? (40pb for the alignement, evalue or perc_identity .. ?) Is my method good ?
Thank's in advance
I recommend to add the code formatting around your code snippets for more clarity
the 200M reads are from one or multiple experiments, whole genome shotgun or 16S amplicon seq?