Additional question regarding BBtools

[steepale]

Hello,
This thread helped me with an issue I'm dealing with only I've become stuck and seek some advice.

I am trying to use BBMap or BBTools to fix what might be either an incorrectly interleaved file or a file that has reads formatted as single-end but really contains a collection of paired and single end reads.

I have interleaved some paired-end RNAseq Illumina reads in order to run them through a program called 'sortmerna' to remove rRNA from my reads. Here is a testformat.sh on the interleaved input file.

/mnt/home/steepale/Apps/bbmap/testformat.sh ./data/017798-1_1_AGTGAG_L001_interleaved_001_100K.fastq
sanger	fastq	raw	interleaved	125bp

If BBTools can provide this ability (filter out rRNA), I would appreciate any advice. The output file seems to either be incorrectly interleaved or are just in single end format. Here is a testformat.sh of the output file I am interested in.

/mnt/home/steepale/Apps/bbmap/testformat.sh ./data/017798-1_1_AGTGAG_L001_norRNA_001_100K.fastq
-1	fastq	raw	single-ended

I would ultimately like to the reads from this output file and separate them by their forward and reverse reads into two files; essentially I want to map them with tophat.

Is anyone familiar with BBTools and how to fix such an issue? I've hit a roadblock.

[GenoMax]

I have moved your question to a new thread to make it visible.

Solution for your problem can come from BBTools themselves. There is a program called bbsplit.sh that you can use with your original data (R1/R2 files). Provide this program with rDNA repeat sequence (provided one is available for your genome). BBsplit can then separate any reads that align to this in one file where as the rest will go to a different file. BBMap is also splice aware so you could use it to align your RNAseq data (should perform better than tophat).

You can use reformat.sh from BBMap to see if the reads are correctly interleaved

$ reformat.sh in=reads.fastq verifypairing

then optionally de-interleave them

$ reformat.sh in=reads.fastq out1=r1.fastq out2=r2.fastq

[steepale]

Here's some additional advice which might add to the convo.

Also, thanks GenoMax, I've located Chicken-specific rDNA clusters and am lifting them over to the correct genome build.

Assuming the reads stayed in the same order but sortmerna just removed some of them, you can use repair.sh like this:

repair.sh in=017798-1_1_AGTGAG_L001_norRNA_001_100K.fastq out1=r1.fq out2=r2.fq outs=single.fq fint
"
If the reads were reordered you'd need the "repair" flag instead of "fint" but they probably were not. The "repair" flag will always work, it just uses more memory than "fint".

However, you can avoid this problem in the first place if you use BBDuk for kmer-matching to remove rRNAs, if you have the ribosomal sequence, since BBDuk will keep pairs together:

bbduk.sh in=interleaved.fq out1=filtered1.fq out2=filtered2.fq outm1=rrna1.fq out2=rrna2.fq ref=ribosomes.fa k=31

You can also use a bulk set of ribosomal sequences like Silva, but using the species' specific ribosomal sequences is much more precise."

[GenoMax]

I would not bother doing any liftover since the sequence of rDNA is unlikely to change between builds.

Use one copy of the full rDNA repeat (don't bother with multiple copies since those are just tandem repeats in most organisms) with whichever tool (bbduk or bbsplit) that you choose to use.