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:

  1. (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
  2. Indexed human genome with bwa

    bwa index GCF_000001405.39_GRCh38.p13_genomic.fa
  3. 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
  4. Transform .sam --> .bam

    samtools view -S -b whole.sam -o whole.bam
  5. Separated unmapped reads (as it is recommended in Materials and Methods using -f4)

    samtools view -f4 whole.bam > sample.unmapped.sam
  6. 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?


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
  • $\begingroup$ Thank you! What about step 0 - converting from .fna to .fa? $\endgroup$ – Dmitrii Trubetskoy Jun 21 '20 at 11:52
  • 3
    $\begingroup$ Unnecessary since this is only a suffix, the content of the file is not changed. You could have simply renamed the file with mv. I am 99.99% sure bwa does not care about the suffix of the fasta. $\endgroup$ – ATpoint Jun 21 '20 at 12:45
  • $\begingroup$ Thanks! I also wonder what are these 5 additional files after indexing, e.g. *.fa.amb, *.fa.ann, *.fa.bwt, *.fa.pac and *.fa.sa ? $\endgroup$ – Dmitrii Trubetskoy Jun 21 '20 at 15:32
  • 1
    $\begingroup$ Internal files that bwa needs for alignment (suffix arrays, burrows-wheeler transformed genomes...), nothing that you as a user have to worry about. Most are binary anyway. $\endgroup$ – ATpoint Jun 21 '20 at 16:12
  • 1
    $\begingroup$ Hi @DmitriiTrubetskoy please remember to upvote ""ATpoint" for their answer and mark the answer as "accepted" if this correctly answers your question. $\endgroup$ – M__ Jun 21 '20 at 19:02

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.