Answers.Body: <p>If you want to keep track of the unmapped reads, I recommend you to use <a href="http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml" rel="nofollow noreferrer">Bowtie2</a> instead of bwa.
  Start by building an index:</p>
  <pre><code>bowtie2-build genome.fa genome.fa
  </code></pre>
  <p>Then map your reads, using arguments to save unmapped reads in a specific file like this:</p>
  <pre><code>bowtie2 -x genome.fa -1 reads_1.fq -2 reads_2.fq -S bowtie2_mapping.sam --un-conc-gz unmappedreads_%.fq.gz
  </code></pre>
  <p>It will write the alignment as well as unmappedreads_1.fq.gz and unmappedreads_2.fq.gz that you may use in DI-tector.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-06-19T11:59:10.800Z
Answers.Text: My thoughts exactly. Personally I don't know why a paper in 2018 would want to use BWA over Bowtie2. Any ideas?
Answers.UserDisplayName: M__
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-06-19T12:30:13.193Z
Answers.Text: 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...
Answers.UserDisplayName: thomas duge de bernonville
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-06-19T16:16:53.087Z
Answers.Text: 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.
Answers.UserDisplayName: Dmitrii Trubetskoy
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-06-19T19:59:10.137Z
Answers.Text: 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.
Answers.UserDisplayName: thomas duge de bernonville
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-06-19T20:25:22.247Z
Answers.Text: @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
Answers.UserDisplayName: Dmitrii Trubetskoy
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-06-20T12:17:34.827Z
Answers.Text: try with `bowtie2 -x human_genome.fa -U /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`, adding the -U argument to indicate the single end read file
Answers.UserDisplayName: thomas duge de bernonville
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-06-20T18:29:46.683Z
Answers.Text: `18035544 reads; of these:
    18035544 (100.00%) were unpaired; of these:
      17979950 (99.69%) aligned 0 times
      55594 (0.31%) aligned exactly 1 time
      0 (0.00%) aligned >1 times
  0.31% overall alignment rate` this is what I have now. I also have unmappedreads_1.fq and unmappedreads_2.fq archives of 1Kb each.. Is that ok? Why are there 2 files? Which one can I use with DI-tector? I also got a `bowtie2_mapping.sam` file which is about 3Gb. Thanks!
Answers.UserDisplayName: Dmitrii Trubetskoy
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-06-21T11:45:02.903Z
Answers.Text: This is because you asked for it `--un-conc-gz unmappedreads_%.fq.gz`. But I don't know if your input data is indeed paired. What is the size of each unmapped*fq?  What input is required for DI-tector?
Answers.UserDisplayName: thomas duge de bernonville
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-06-21T20:14:13.357Z
Answers.Text: Thanks again, and sorry for taking your time. There are two unmapped reads and the size is 1Kb each. In fact, they are empty. DI-tector requires unmatched reads in .fastq format + Virus_ref genome (I believe it's an indexed .fasta viral genome)
Answers.UserDisplayName: Dmitrii Trubetskoy
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-06-22T06:50:15.260Z
Answers.Text: ok I think it is my fault, you do not have paired-end reads but single end. So the correct option should be `--un-gz`, instead of `--un-conc-gz`, as described in the [bowtie2](http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml) manual.
Answers.UserDisplayName: thomas duge de bernonville
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-06-19T04:24:31.920Z
Answers.LastActivityDate: 2020-06-19T04:24:31.920Z
Answers.UserDisplayName: thomas duge de bernonville
Body: <p>I am trying to follow the Beauclair et al. paper (free version here <a href="https://rnajournal.cshlp.org/content/24/10/1285.long" rel="nofollow noreferrer">https://rnajournal.cshlp.org/content/24/10/1285.long</a>) for identifying defective genomes using their DI-tector program.</p>
  <p>According to figure 2 <a href="https://i.sstatic.net/8e1SZ.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/8e1SZ.jpg" alt="enter image description here" /></a></p>
  <p>Before using the program itself, it is highly recommended to filter my raw read using <strong>bwa</strong> and then <strong>samtools</strong></p>
  <p>Unfortunately I have no idea where to start - I have browsed official bwa and samtools sites, I have installed them using the following commands:</p>
  <pre><code>sudo apt-get update -y
  sudo apt-get install -y bwa
  </code></pre>
  <p>I have indexed my reference genome with the following command (as I read it is required before the analysis)</p>
  <pre><code>bwa index sequence.fasta
  </code></pre>
  <p>Then, I have installed samtools</p>
  <pre><code>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
  </code></pre>
  <p>I wonder what will be the further steps to get these &quot;UNMATCHED READS at the end of step i&quot;, since Di-tector program takes virus_reference genome (I have it) and these unmatched reads.</p>
  <p>First of all I want to align my .fastq files vs <code>GCF_000001405.39_GRCh38.p13_genomic.fna</code> - 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?</p>
  
Tags: <samtools><bwa>
Title: How can I get unmatched reads for defective genomes analysis using bwa and samtools?
ThreadText: # How can I get unmatched reads for defective genomes analysis using bwa and samtools?
  
  (Asked by: username_0 on 2020-06-18)
  
  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](https://i.sstatic.net/8e1SZ.jpg)](https://i.sstatic.net/8e1SZ.jpg)
  
  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?
  
  **Tags:** samtools, bwa
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2020-06-19. Score: 1
  
  If you want to keep track of the unmapped reads, I recommend you to use [Bowtie2](http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml) 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.
  
  > **Comments:**
  > **username_2**: My thoughts exactly. Personally I don't know why a paper in 2018 would want to use BWA over Bowtie2. Any ideas?  
  > **username_1**: 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...  
  > **username_0**: 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.  
  > **username_1**: 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.  
  > **username_0**: @username_1 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  
  > **username_1**: try with `bowtie2 -x human_genome.fa -U /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`, adding the -U argument to indicate the single end read file  
  > **username_0**: `18035544 reads; of these:   18035544 (100.00%) were unpaired; of these:     17979950 (99.69%) aligned 0 times     55594 (0.31%) aligned exactly 1 time     0 (0.00%) aligned >1 times 0.31% overall alignment rate` this is what I have now. I also have unmappedreads_1.fq and unmappedreads_2.fq archives of 1Kb each.. Is that ok? Why are there 2 files? Which one can I use with DI-tector? I also got a `bowtie2_mapping.sam` file which is about 3Gb. Thanks!  
  > **username_1**: This is because you asked for it `--un-conc-gz unmappedreads_%.fq.gz`. But I don't know if your input data is indeed paired. What is the size of each unmapped*fq?  What input is required for DI-tector?  
  > **username_0**: Thanks again, and sorry for taking your time. There are two unmapped reads and the size is 1Kb each. In fact, they are empty. DI-tector requires unmatched reads in .fastq format + Virus_ref genome (I believe it's an indexed .fasta viral genome)  
  > **username_1**: ok I think it is my fault, you do not have paired-end reads but single end. So the correct option should be `--un-gz`, instead of `--un-conc-gz`, as described in the [bowtie2](http://bowtie-bio.sourceforge.net/bowtie2/manual.shtml) manual.  
  
  ---
  
  