Answers.Body: <blockquote>
  <p>Can I expect that after running the command <code>bowtie2 -x ref_index -U seq.fastq.gz -S seq.sam</code>, my <code>seq.sam</code> will only contain the DNA that I'm interested in (i.e., Human mtDNA)?</p>
  </blockquote>
  <p>Not quite as that depends on what you had in your input. If you had a lot of non-human mitochondrial sequences, for example, those might still be there if they align well enough against your reference. Essentially, what you can assume here is that the SAM file will have those input sequences that aligned well against your mt DNA reference. Whether that will be a clean set of human mitochondrial sequences depends entirely on what you had as input going in.</p>
  <p>However, as a general rule, yes, aligning against a specific target should leave you with the reads that align to that target which in your case means human mitochondrial sequences.</p>
  <hr />
  <p>And here is a brief explanation of what each command does:</p>
  <ol>
  <li><p><code>bowtie2 -p $THREADS -x host_DB -1 $FW_READS -2 $RV_READS --un-conc-gz SAMPLE_host_removed &gt; $STUDY-data.sam </code></p>
  <p>This will align the input fastq sequences against the reference to produce a <a href="https://samtools.github.io/hts-specs/SAMv1.pdf" rel="nofollow noreferrer">SAM file</a>.</p>
  </li>
  <li><p><code>samtools view -bS $STUDY-data.sam &gt; $STUDY-data.bam </code></p>
  <p>This simply converts the SAM file to a BAM file. Essentially a compressed SAM.</p>
  </li>
  <li><p><code>samtools view -b -f 12 -F 256 $STUDY-data.bam &gt; $STUDY-data-no-human.bam</code></p>
  <p>Here, you are filtering your bam. The <code>-b</code> flag just means the output will be BAM, and the <code>-f</code> and <code>-F</code> control what alignments are kept or discarded respectively depending on the SAM flag. This... is a bit complicated, you need to read the SAM specification to understand fully, but briefly, each alignment in a SAM/BAM file has a binary flag whose value gives you information about the alignment. You can use the helpful &quot;<a href="https://broadinstitute.github.io/picard/explain-flags.html" rel="nofollow noreferrer">Explain SAM flags</a>&quot; page from the Broad Institute to understand what each number means.</p>
  <p>12 means &quot;read unmapped and mate unmapped&quot; so <code>-f 12</code> means &quot;keep only alignments where the bits given are set&quot;. this will match read pairs where neither pair has been mapped. Next, <code>-F</code> is the inverse option and means &quot;remove alignments where the bits given are set&quot;. 256 means &quot;this is not a primary alignment&quot;. Taken together then, this command would only keep unmapped reads, excluding reads that are unmapped as a secondary alignment. In other words, it will result in a bam file with only those reads that couldn't be mapped at all. Presumably, in the example you found this from, this is to collect the non-human reads, i.e. those that don't match the reference genome.</p>
  </li>
  <li><p><code>samtools sort -n -o $STUDY-sorted-data.bam $STUDY-data-no-human.bam 1&gt;&gt; $REPORTS_DIR/report_stdout.txt</code></p>
  <p>This just sorts the bam file. It doesn't remove anything.</p>
  </li>
  <li><p><code>samtools fastq -1 $STUDY-reads-r1.fq -2 $STUDY-reads-r2.fq $STUDY-sorted-data.bam 1&gt;&gt; $REPORTS_DIR/report_stdout.txt</code></p>
  <p>And this just extracts fastq files from the BAM alignment.</p>
  </li>
  </ol>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-04-09T15:56:04.120Z
Answers.LastActivityDate: 2024-04-09T17:47:21.437Z
Answers.LastEditDate: 2024-04-09T17:47:21.437Z
Answers.UserDisplayName: terdon
Body: <p>I'd like to align some FASTQs to an average mtDNA FASTA file that I have downloaded so I can have the human mtDNA isolated from those FASTQs. For that, I used bowtie2. Can I expect that after running the command <code>bowtie2 -x ref_index -U seq.fastq.gz -S seq.sam</code>, my <code>seq.sam</code> will only contain the DNA that I'm interested in (i.e., Human mtDNA)?</p>
  <p>I have this example below from a class in which is supposed that the human DNA is discarded (Opposite of what I want). But I don't understand the roles of the line 3. Is it necessary to filter out undesired DNA?</p>
  <ol>
  <li><p><code>bowtie2 -p $THREADS -x host_DB -1 $FW_READS -2 $RV_READS --un-conc-gz SAMPLE_host_removed &gt; $STUDY-data.sam </code></p>
  </li>
  <li><p><code>samtools view -bS $STUDY-data.sam &gt; $STUDY-data.bam </code></p>
  </li>
  <li><p><code>samtools view -b -f 12 -F 256 $STUDY-data.bam &gt; $STUDY-data-no-human.bam</code></p>
  </li>
  <li><p><code>samtools sort -n -o $STUDY-sorted-data.bam $STUDY-data-no-human.bam 1&gt;&gt; $REPORTS_DIR/report_stdout.txt</code></p>
  </li>
  <li><p><code>samtools fastq -1 $STUDY-reads-r1.fq -2 $STUDY-reads-r2.fq $STUDY-sorted-data.bam 1&gt;&gt; $REPORTS_DIR/report_stdout.txt</code></p>
  </li>
  </ol>
  
Tags: <phylogenetics><samtools><fastq><bowtie2>
Title: Aligning FASTQs to FASTA reference
ThreadText: # Aligning FASTQs to FASTA reference
  
  (Asked by: username_0 on 2024-04-09)
  
  I'd like to align some FASTQs to an average mtDNA FASTA file that I have downloaded so I can have the human mtDNA isolated from those FASTQs. For that, I used bowtie2. Can I expect that after running the command `bowtie2 -x ref_index -U seq.fastq.gz -S seq.sam`, my `seq.sam` will only contain the DNA that I'm interested in (i.e., Human mtDNA)?
  
  I have this example below from a class in which is supposed that the human DNA is discarded (Opposite of what I want). But I don't understand the roles of the line 3. Is it necessary to filter out undesired DNA?
  
  1. `bowtie2 -p $THREADS -x host_DB -1 $FW_READS -2 $RV_READS --un-conc-gz SAMPLE_host_removed > $STUDY-data.sam`
  2. `samtools view -bS $STUDY-data.sam > $STUDY-data.bam`
  3. `samtools view -b -f 12 -F 256 $STUDY-data.bam > $STUDY-data-no-human.bam`
  4. `samtools sort -n -o $STUDY-sorted-data.bam $STUDY-data-no-human.bam 1>> $REPORTS_DIR/report_stdout.txt`
  5. `samtools fastq -1 $STUDY-reads-r1.fq -2 $STUDY-reads-r2.fq $STUDY-sorted-data.bam 1>> $REPORTS_DIR/report_stdout.txt`
  
  **Tags:** phylogenetics, samtools, fastq, bowtie2
  
  
  ---
  
  ## 1 answers
  
  ### Answer by username_1 on 2024-04-09. Score: 5 | Accepted answer
  
  > Can I expect that after running the command `bowtie2 -x ref_index -U seq.fastq.gz -S seq.sam`, my `seq.sam` will only contain the DNA that I'm interested in (i.e., Human mtDNA)?
  
  Not quite as that depends on what you had in your input. If you had a lot of non-human mitochondrial sequences, for example, those might still be there if they align well enough against your reference. Essentially, what you can assume here is that the SAM file will have those input sequences that aligned well against your mt DNA reference. Whether that will be a clean set of human mitochondrial sequences depends entirely on what you had as input going in.
  
  However, as a general rule, yes, aligning against a specific target should leave you with the reads that align to that target which in your case means human mitochondrial sequences.
  
  ---
  
  And here is a brief explanation of what each command does:
  
  1. `bowtie2 -p $THREADS -x host_DB -1 $FW_READS -2 $RV_READS --un-conc-gz SAMPLE_host_removed > $STUDY-data.sam`
  
     This will align the input fastq sequences against the reference to produce a [SAM file](https://samtools.github.io/hts-specs/SAMv1.pdf).
  2. `samtools view -bS $STUDY-data.sam > $STUDY-data.bam`
  
     This simply converts the SAM file to a BAM file. Essentially a compressed SAM.
  3. `samtools view -b -f 12 -F 256 $STUDY-data.bam > $STUDY-data-no-human.bam`
  
     Here, you are filtering your bam. The `-b` flag just means the output will be BAM, and the `-f` and `-F` control what alignments are kept or discarded respectively depending on the SAM flag. This... is a bit complicated, you need to read the SAM specification to understand fully, but briefly, each alignment in a SAM/BAM file has a binary flag whose value gives you information about the alignment. You can use the helpful "[Explain SAM flags](https://broadinstitute.github.io/picard/explain-flags.html)" page from the Broad Institute to understand what each number means.
  
     12 means "read unmapped and mate unmapped" so `-f 12` means "keep only alignments where the bits given are set". this will match read pairs where neither pair has been mapped. Next, `-F` is the inverse option and means "remove alignments where the bits given are set". 256 means "this is not a primary alignment". Taken together then, this command would only keep unmapped reads, excluding reads that are unmapped as a secondary alignment. In other words, it will result in a bam file with only those reads that couldn't be mapped at all. Presumably, in the example you found this from, this is to collect the non-human reads, i.e. those that don't match the reference genome.
  4. `samtools sort -n -o $STUDY-sorted-data.bam $STUDY-data-no-human.bam 1>> $REPORTS_DIR/report_stdout.txt`
  
     This just sorts the bam file. It doesn't remove anything.
  5. `samtools fastq -1 $STUDY-reads-r1.fq -2 $STUDY-reads-r2.fq $STUDY-sorted-data.bam 1>> $REPORTS_DIR/report_stdout.txt`
  
     And this just extracts fastq files from the BAM alignment.
  
  
  ---
  
  