bbmap split paired-end reads back into separated fastq files?

Question

[c.e.chong]

Hi,

I am trying to remove human reads from my metagenomic sequences (skin samples).

I have used the command recommended in the bbmap documentation:

bbmap.sh minid=0.95 maxindel=3 bwr=0.16 bw=12 quickmatch fast minhits=2 path=/path/to/hg19masked/ qtrim=rl trimq=10 untrim -ea -Xmx100g in=reads.fq outu=clean.fq outm=human.fq threads=12

My input is Illumina HiSeq 4000 Paired-end, 2x150 bp reads, so I want to split paired-end reads back into separated fastq files .._r1 .._r2.

Is there a way to do this using BBTools?

Thanks in advance!

[GenoMax]

As shown your command is only using one read file (in=). So unless that file is interleaved there is no chance that out*= files will be interleaved. @ATPoint shows what an interleaved file should look like. Note: Your headers may look different but Read 1 and 2 from one clusters should be one after the other in the file.

[c.e.chong]

Sorry that was a mistake, I have actually used in=r1.fq in2= r2.fq!

Thanks!

Answer 1

[GenoMax]

Use reformat.sh from BBMap suite. This assumes your clean.fq is interleaved.

reformat.sh in=clean.fq out1=clean_R1.fq out2=clean_R2.fq

[c.e.chong]

Thank you very much for your help!

How would I find out if the clean.fq file is interleaved?

[ATpoint]

Check if two adjacent read pairs always have the same read name. This is an example of an interleaved file from a HiSeq3000:

@HWI-ST-J00104_BSF_0515:8:1101:10003:10036#DMSO_rep2_S43794/1
GGAAGAAAGACAGTCTGTGGCCCTGCCTGGGGACCTACACTGTCTGCTGTAACAGGCTTTCCCTTCATCTCAAGAG
+
AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ
@HWI-ST-J00104_BSF_0515:8:1101:10003:10036#DMSO_rep2_S43794/2
GGTTGAAAAGTGGCCTCGGTGTTCCAGACCTACCTGGTTCACTTAGCTTTTTCCTCCTTTCCTTCCTTTTATACTC
+
AAFFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJJJJJJJJJJJJJJFJJJJJJJFJJJJJJJJJ
@HWI-ST-J00104_BSF_0515:8:1101:10003:10141#DMSO_rep2_S43794/1
GTGATTCTCTTCAAGGCCACTCCTAAACAACTGTTGAATCCTTGATCCAGGCACCAAGCCCAGAGGTCCCTATCAG
+
AAAFFJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ<FJJJJJJJJJJJJJFJJJJJJJJJJJ
@HWI-ST-J00104_BSF_0515:8:1101:10003:10141#DMSO_rep2_S43794/2
GTCTAGGAGGAAACCAGGCTGTTAACTCCTCAGCAAAGACAAAGGAGGAGTTGGGCGGGGGCAAGACTAACTCCTA
+
AAFFFJJJJJFJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJJAJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ

[robert.murphy]

Sorry to necro this but is there a difference in using reformat.sh vs bbsplit.sh to do this operation?

[GenoMax]

reformat.sh is simply separating R1/R2 reads that are interleaved as shown in @ATPoint's example above. bbsplit.sh on the other hand uses reference(s) genome(s) to bin the reads depending on how they align to those genomes.

[robert.murphy]

can reformat.sh also be used to separate R1/R2 reads that are merged ?

[GenoMax]

No program can do that if the reads are already merged. Use following alternate method.

You could get the fastq headers from the merged files and then use them to retrieve reads from original R1/R2 files using filterbyname.sh program from BBMap suite.

[bioknown]

Again sorry to necro this post, but is the correct way to work with paired end data to use in1= in2= and then have a single outu and outm , then run reformat.sh to get the paired unmapped reads when trying to remove human contaminants?

[GenoMax]

You can specify two output files as in outu1= outu2= etc. If you provide only one output file then the data will be written interleaved. Not many programs understand interleaved data and thus you will need to separate them into two files using reformat.sh.

[Kermit]

It's validating that reformat.sh is the approach that the library author recommends here: Tool to separate human and mouse rna seq reads

[nasreenbano91]

Can reformat.sh also be used to merge the replicate file?

like-

reformat.sh in1=file1_R1.fastq in2=file1_R2.fastq in1=file2_R1.fastq in2=file2_R2.fastq out=merged_R1.fastq out2=merged_R2.fastq

where file1_R1.fastq and file1_R2.fastq are the paired-end reads from the first replicate, and file2_R1.fastq and file2_R2.fastq are the paired-end reads from the second replicate?

Is it the correct way to merge the replicate?

[Ram]

How is this question an answer to the original question? I'm moving it to a comment this time, please be more mindful in the future.

[GenoMax]

If you simply want to merge "technical" sequencing replicates then you can do the following

cat file1_R1.fastq file1_R1.fastq > merged_R1.fastq
cat file1_R2.fastq file2_R2.fastq > merged_R2.fastq