Aligning FASTQs to FASTA reference

Question

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

Answer 1

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.

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" 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.