Merging bowtie output

Question

[Sbrillo]

Hi,

I'm working on a pipeline that needs illumina and pacbio reads as input.

The first step consists of the removal of mt reads from the illumina reads.

To remove the reads i used this command

bowtie2 -p 30 -x /MtRef_bowtie2 -1 SRR1508169_1.fastq.gz -2 SRR1508169_2.fastq.gz --un-conc output1/f.fastq 2> f_bowtie2.log

I have two files as output (f.1.fastq and f.2.fastq) as expected. However the pipeline require a single fastq file as an input.

DO you know how I can obtain a single file output from bowtie2? Do you think is possible to just merge f.1 and f.2 files just using cat command?

Best

[ATpoint]

Hello, I am not quite sure I understand. So you want to remove reads of mt origin from this paired-end dataset, is that correct? For this you map to the mt sequence, guess that makes sense, but then what is the output you need? Is it a single fastq file, or again two fastq files, and if a single one, into what kind of tool or command does this go then, and why does it need a single file? You could make a so-called interleaved fastq file, which means that R1 and R2 are in alternating order, whereas a naive cat (while technically possible) most likely is not what you want as the two reads belong to a pair and order must be preserved for them to be meaningful. Can you elaborate?

Edit: I guess it makes more sense to map against the entire genome incl the mt sequence and then extract reads not mapping to mt. If you only align to mt then homologs and similar sequences with true origin being the genome might be incorrectly removed in your approach.

[Sbrillo]

Thanks for replying!

Yes, i need to remove mt from paired-end fastq file and i need a fastq file (not bam) as input for the next steps. I need a single file because i need to map the illumina reads from males and females samples to the pacbio reads in order to create chromosomal bins and the pipeline is designed to work with one single input from males and females (not 4 in in the case i used R1 and R2 fastq files).

You suggestion of mapping the illumina reads to the whole genome makes definitely sense and i will try.

Do you think that i should use a specific tool like pandaseq to interveal the reads? Or i have another strategy in my mind.

Use pandaseq BEFORE the step of mt reads removal and run bowtie2 with this parameters:

bowtie2 -p 30 -x /MtRef_bowtie2 --intervealed f.merged.fastq --un-conc output1/f.fastq 2> f_bowtie2.log

DO you think is gonna work?

[GenoMax]

i need to remove mt from paired-end fastq file and i need a fastq file (not bam) as input for the next steps.

You can use bbduk.sh from BBMap suite in filter mode. You will use mitochondrial genome to filter the reads against. A guide for BBDuk is available here.

If you want to be extra careful you could use bbsplit.sh with your genome plus mitochondrial genome to bin the reads. Something like:

bbsplit.sh in1=reads1.fq in2=reads2.fq ref=genome.fa,mito.fa basename=out_%.fq outu1=clean1.fq outu2=clean2.fq

[Mark]

Use the samtools merge to merge the bam files. Do not sure cat command.