Can you assemble with merged paired end reads and unmatched reads as "single ended" reads?

Question

[O.rka]

I used BBMap to isolate r1 and r2 reads while providng unmatched reads as well. It generated a single file. It is wise to use this in my assembler (SPAdes) as single-ended?

I ended up using the outm=./bbmap_output/mapped.fq from below as input into SPAdes as single ended with the -s flag. My mapped.fq file has not only paired r1 and r2 reads but also the singletons.

bbwrap.sh -Xmx40g in=./reads/r1.fq,$SINGLETONS in2=./reads/r2.fq,null outm=./bbmap_output/mapped.fq outu=./bbmap_output/unmapped.fq ref=./reference/assembly.fa out=./bbmap_output/output.sam lengthtag=t idtag=t covstats=./bbmap_output/output.covstats.txt rpkm=./bbmap_output/output.rpkm.txt threads=$N_JOBS usemodulo append

For SPAdes, my command was the following:

python spades-3.9.0/bin/spades.py -t $N_JOBS -s ./bbmap_output/mapped.fq -o ./spades_output/

[GenoMax]

I used BBMap to isolate r1 and r2 reads. It generated a single file

Which tool from BBMap did you use? reformat.sh? If yes then you must have made a interleaved reads file.

You can look at SPAdes manual section 3.1 on how to specify interleaved reads.

--12 <file_name>
    File with interlaced forward and reverse paired-end reads.

[O.rka]

@genomax thank you. I've added my command above in the original description. Are they still interlaced if it includes both paired r1/r2 and the unpaired reads?

[GenoMax]

You used bbmap.sh to align your data to assembly.fa. You obtained a SAM format file output.sam from that alignment, which is no longer in the fastq format but is in SAM format. You seem to have captured reads that did not align in unmapped.fq and that is probably in the interleaved format since you started with PE reads but provided only single output file name.

If you wanted to assemble the data then there is no need to align it first. You should start with your scanned/trimmed original data files and then go into SPAdes directly after that.

Note: If this question is a related to/follow-up on Extracted mapped shotgun metagenomic reads to reference genome. SPAdes or metaSPAdes for de-novo assembly? then you are on the right track. You can retrieve reads that mapped to your assembly by doing

reformat.sh in=output.sam out1=R1.fq out2.R2.fq

You can then use these reads in your SPAdes assembly.

[O.rka]

Apologies, I forgot to add the bbwrap.sh command that I ended up using. I've updated the question with the correct command and my input into spades. Can reformat.sh create singletons? If I was able to regenerate the cleaned r1.fq, r2.fq, and singletons.fq fastq files, would assembling these with spades.py -1 r1.fq -2 r2.fq -s singletons.fq result in a better assembly than spades.py -s mapped.fq where mapped.fq has everything merged into one?

[GenoMax]

Unless you have a really high number of singletons, I suggest that you ignore them for now and try an assembly with properly paired reads.

[pdhrati02]

GenoMax Hi, sorry to reply on such an old post but I have a similar question. The kneaddata output has unmatched_1.fq and unmatched_2.fq which are reads whose mates are lost but they themselves passed both trimmomatic and bowtie2 step. In this case would at what step would have reads without mates be an issue in downstream processing? Thanks in advance

[GenoMax]

You will need to provide some context. What exactly are you trying to do?

[pdhrati02]

My apologies, I wish to use these sequences for taxonomic and functional profiling (plan to use metaphlan and humann) and also to assemble using megahit/metaspades. I wish to know if using those unmatched reads will affect these steps?

[GenoMax]

I am not sure which steps you are referring to. If programs you are planning to use can use the singleton reads then you would be able to use them. If they require paired-end data then these reads are not going to be useful.

Answer 1

[O.rka]

An adaptation from @genomax. I used BBMap for the paired separately and the de-interlaced the reads.

 bbmap.sh in=./reads/r1.fq in2=./reads/r2.fq ref=./reference/assembly.fa outm=$R1_R2 lengthtag=t idtag=t threads=$N_JOBS usemodulo
 reformat.sh in=$R1_R2 out1=$R1 out2=$R2
 bbmap.sh in=$SINGLETONS ref=./reference/assembly.fa outm=./bbmap_output/singletons.mapped.fq lengthtag=t idtag=t threads=$N_JOBS usemodulo
 python spades.py -t $N_JOBS -1 $R1 -2 $R2 -s $SINGLETONS -o ./spades_output/