1
$\begingroup$

I am trying to follow the 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.

According to figure 2 enter image description here

Before using the program itself, it is highly recommended to filter my raw read using bwa and then samtools

Unfortunately I have no idea where to start - I have browsed official bwa and samtools sites, I have installed them using the following commands:

sudo apt-get update -y
sudo apt-get install -y bwa

I have indexed my reference genome with the following command (as I read it is required before the analysis)

bwa index sequence.fasta

Then, I have installed samtools

cd ..
wget https://github.com/samtools/samtools/releases/download/1.9/samtools-1.9.tar.bz2
tar -vxjf samtools-1.9.tar.bz2
cd samtools-1.9
make
sudo apt install samtools

I wonder what will be the further steps to get these "UNMATCHED READS at the end of step i", since Di-tector program takes virus_reference genome (I have it) and these unmatched reads.

First of all I want to align my .fastq files vs GCF_000001405.39_GRCh38.p13_genomic.fna - what will be the output file? If I got it right - I need this output file to filter with samtools, am I correct? If I am not mistaken, the output after samtools will be the .sam file that I will have to put into a DI-tector - correct?

Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

If you want to keep track of the unmapped reads, I recommend you to use Bowtie2 instead of bwa. Start by building an index:

bowtie2-build genome.fa genome.fa

Then map your reads, using arguments to save unmapped reads in a specific file like this:

bowtie2 -x genome.fa -1 reads_1.fq -2 reads_2.fq -S bowtie2_mapping.sam --un-conc-gz unmappedreads_%.fq.gz

It will write the alignment as well as unmappedreads_1.fq.gz and unmappedreads_2.fq.gz that you may use in DI-tector.

Improve this answer
$\endgroup$
10
  • 1
    $\begingroup$ My thoughts exactly. Personally I don't know why a paper in 2018 would want to use BWA over Bowtie2. Any ideas? $\endgroup$
    – M__ ♦
    Commented Jun 19, 2020 at 11:59
  • 1
    $\begingroup$ after putting a glance at the article, it is no really clear how they switch from bwa mem to bwa samse to align reads. Maybe it has something to see with the flags in the output .sam file, required for DI-tector proper running. But it is only a supposition... $\endgroup$
    – thomas duge de bernonville
    Commented Jun 19, 2020 at 12:30
  • $\begingroup$ Thank you. But what about samtools step? I mean, my read has to be aligned with bwa (or bowtie2 - great I will try that) FIRSTLY with the human genome, then SECONDLY with virus genome and only after these two steps it has to be filtered through SAMTOOLS and only after that I have to have the unmapped reads that I can put into a DI-tector. $\endgroup$
    – Dmitrii Trubetskoy
    Commented Jun 19, 2020 at 16:16
  • $\begingroup$ I think it depends on the format of input data required for DI-tector. Bowtie2 will return .fastq files. In the article, they used samtools with -Bs -f 4 options probably to return unaligned reads into a .sam with specific CIGAR labels. $\endgroup$
    – thomas duge de bernonville
    Commented Jun 19, 2020 at 19:59
  • $\begingroup$ @thomas duge de bernonville OK, I did bowtie2-build GCF_000001405.39_GRCh38.p13_genomic.fna human_genome.fa first, and the output is 4 files human_genome.fa.1.bt2, human_genome.fa.2.bt2, human_genome.fa.3.bt2, human_genome.fa.4.bt2. I wonder how can I map my reads using these 4 files while in the command you have proposed you have only genome.fa Anyways I have already tried this bowtie2 -x human_genome.fa /mnt/e/nastya/SLX066-02/B-dVMV-RIG-1/B-dVMV-RIG-1_ACAGTG_L008_R1_001.rc.fastq.gz -S bowtie2_mapping.sam --un-conc-gz unmappedreads_%.fq.gz and it gave me an (ERR) bowtie2 exited $\endgroup$
    – Dmitrii Trubetskoy
    Commented Jun 19, 2020 at 20:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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