To map reads against the assembled draft genome, we will use Bowtie2. Bowtie2 defaults to finding global alignments (no soft-clipping of terminal bases) and it allows for gaps in the alignment, thereby increasing mapping accurracy (Schlotterer et al. (2014) Nature Reviews Genetics). Although the aligner TMAP was specifically compiled for Ion Torrent data, in this tutorial we want use an aligner that can also be used with Illumina or SoliD data - as you are likely to work with sequencing data that were obtained from these platforms in the future. An alternative aligner that is currently widely used is BWA.

What you need for mapping are two fasta files, one containing the quality-trimmed reads and the other containing the draft genome. Create a new folder named Mapping (with mkdir) and copy the two fasta files into it (with cp). For example:

mkdir Mapping

cp GenomeAssembly/IonTorrentDeNovoAssembly_d_results/\
IonTorrentDeNovoAssembly_out.unpadded.fasta \
Mapping/

cp RAWREADS_trimmed_trimmed.fq \
Mapping/

cd Mapping

In the Mapping folder, we first create an index for the reference genome using the following command:

bowtie2-build -f IonTorrentDeNovoAssembly_out.unpadded.fasta GENOMENAME

You can change GENOMENAME to whatever you like.

Now, to align the trimmed reads against the genome, use the following command:

nohup bowtie2 -p 1 \
-q \
--phred33 \
--sensitive \
--no-unal \
--al-gz FILE_Aligned.fq.gz \
--un-gz FILE_Unaligned.fq.gz \
--met-file MetricsFile.txt \
-x GENOMENAME \
-U RAWREADS_trimmed_trimmed.fq \
-S MAPPEDREADS.sam >Logfile.log &

Here an overview of the meaning of the used options:

-p 1

Causes Bowtie 2 to use a single thread. Depending on the number of users and libraries we will probably increase this.

-q

Informs the program that the reads to map are saved in fastq files.

--phred33

Sets the quality encoding of the fastq files to “Phred+33”.

--sensitive

sets several options at once regarding the seeding and other adjustments.

--no-unal

Suppress SAM records for reads that do not align.

--al-gz

Write unpaired reads that align at least once to to the specified file.

--un-gz

Write unpaired reads that failed to align to the specified file.

--met-file

Write bowtie2 metrics to Metricsfile.txt.

-x

Specifies the name of the genome.

-U

Specify the unpaired reads to align (can contain a comma-separated list of several fq files).

-S

Specify the sam file to which the alignment shall be saved.

You can’t set the exact number of mismatches in the seed, but you can adjust the mismatch penalty.

The program should run no longer than 10-20 mins. The resulting output file will be in the SAM format. For a detailed description of this format, see here.

To remove unmapped reads, reads below a mapping quality of 20, and reads that were not aligned uniquely (reads that were mapped to >1 places in the genome), use the python script Bowtie2Filtering.py:

Bowtie2Filtering.py -mq -u -a -s MAPPEDREADS.sam

Your filtered reads will be saved in MAPPEDREADSfiltered.sam

Alternatively, you can use samtools to filter out reads with a mapping quality <20:

samtools view -Sh -q 20 -o MAPPEDREADS_QualityAbove20.sam MAPPEDREADS.sam

Options:

-S

Input is in the sam format

-h

Include the samfile header in the output

-q

Skip alignments with a mapping quality below 20

A widely used tool to identify sequence variants is samtools mpileup. See here for an overview of its options. The tool samtools mpileup defaults to creating a pileup file, which summarizes aligned base calls in text format (see here for a detailed characterization of a pileup file http://samtools.sourceforge.net/pileup.shtml). If you call samtools mpileup with the -u or -g option, instead, the output format is a vcf or bcf (compressed binary version of vcf) file; vcf stands for ‘variant call format’. Its format specifications are described here and summarized in Fig. 3.

The first step for calling SNPs from your aligned and deduplicated reads is:

samtools mpileup -g \
-f \
IonTorrentDeNovoAssembly_out\
.unpadded.fasta \
-q 20 \
-Q 20 \
-t DP \
-t SP \
MAPPEDREADS_dedup.bam  > MAPPEDREADS_dedup.bcf

The chosen options are described on this page. By setting the -t SP and -t DP tags, samtools mpileup provides:

-t SP

per-sample Phred-scaled strand bias P-value

-t DP

per sample read depth

To call SNPs from the bcf file, we use bcftools:

bcftools call -vm -V indels MAPPEDREADS_dedup.bcf >  MAPPEDREADS_variants.vcf

Options:

-v

Output variant sites only

-V indels

Skip indels

-m

model for multiallelic and rare-variant calling

To count how many SNPs were found, use the following command:

grep -v -c '^#' MAPPEDREADS_variants.vcf

The option -v in combination with ^# excludes all header lines that start with (^) the #-sign. With the -c option, grep counts the lines instead of writing them out.

To filter out SNPs that are low quality or covered by low depth, we can use the vcfutils.pl varFilter that comes with samtools:

vcfutils.pl varFilter -d 5 -w 3 -Q 20  MAPPEDREADS_variants.vcf > MAPPEDREADS_variants_filtered.vcf

Options used:

-d 5

minimum read depth of 5

-w 3

SNP within 3 bp around a gap to be filtered. This may be an alternative solution to re-alignment around indels

-Q 20

minimum mapping quality of 20

Another useful option can be:

-1 0.0001

min P-value for strand bias (given the PV4-tag in the vcf file). We obtained the PV4-tag by setting the -t SP tag in samtools mpileup. This option filters out the SNPs that have a strong strand-bias: SNPs that are supported by one strand and not the other.

Count how many SNPs are left after filtering

grep -v -c '^#' MAPPEDREADS_variants_filtered.vcf

The SNPs can be visualized with IGV. For this, we first need to compress and index the vcf files:

bgzip -c \
MAPPEDREADS_variants_filtered.vcf \
> MAPPEDREADS_variants_filtered.vcf.gz

tabix \
-p vcf \
MAPPEDREADS_variants_filtered.vcf.gz

Open IGV and load the indexed bam file and the indexed vcf file.