How to extract unmatched reads using bwa and samtools?

Question

I have a single read (NOT paired) that I need to pass through the workflow described in Beauclair et al. paper (free version here https://rnajournal.cshlp.org/content/24/10/1285.long) for identifying defective genomes using their DI-tector program.

Here in Materials and methods the procedure is described as the following:

The first step of the workflow consists of an alignment of the reads against the host genome (Fig. 2i). This step aims at discarding reads that map to the host genome and may partially map to the viral genome after segmentation, and at reducing the working file size. For example, MV and rMV-ΔV data sets were generated from total RNA samples of infected cells and contained mostly reads mapping the human genome (≈99%). This step uses a combination of bwa mem and samtools view with the parameters –bS –f4. An additional step consists of an alignment of the reads against the viral genome of interest, in order to exclude reads perfectly mapping the viral genome. Therefore, only unmapped reads are further analyzed. Of note, clipped reads (i.e., CIGAR motif contain S or H) are also conserved. Some of these reads may map to viral genome recombination junctions that are present in DI genomes.

I have already been suggested to use bowtie2 instead of bwa but, firstly, the output is not clear for me and secondly, I would like to test the official protocol.

Since the article suggests to use bwa and samtools on this very first step, that is what I have done so far:

0. (Optional) not sure if this is important but as suggested by someone I have transformed .fna --> .fa

   cp GCF_000001405.39_GRCh38.p13_genomic.fna GCF_000001405.39_GRCh38.p13_genomic.fa
   

1. Indexed human genome with bwa

   bwa index GCF_000001405.39_GRCh38.p13_genomic.fa
   

2. aligned my only read with indexed genome from step 1

   bwa mem GCF_000001405.39_GRCh38.p13_genomic.fa /mnt/e/nastya/SLX066-02/B-dVMV-RIG-1/B-dVMV-RIG-1_ACAGTG_L008_R1_001.rc.fastq.gz > whole.sam
   

3. Transform .sam --> .bam

   samtools view -S -b whole.sam -o whole.bam
   

4. Separated unmapped reads (as it is recommended in Materials and Methods using -f4)

   samtools view -f4 whole.bam > sample.unmapped.sam
   

5. Converted unmapped reads into .fastq format (since this is the format used by the software later)

   samtools fastq sample.unmapped.sam > unmatched.fastq
   

I wonder if these steps are correct?

Answer 1

Yes, correct but overly many steps. There is no need to convert between bam and sam. samtools can read from stdin and handles both sam and bam and samtools fastq can interpret flags, therefore one can shorten this to:

bwa mem (...options) | samtools view -o out.bam
samtools fastq -f 4 out.bam > unmatched.fastq