What tool to use to map reads (samtools, sequencher,...)?

Question

I am a beginner in Bioinformatics. I have 4 files:

2 fastq files (aln1.fastq and aln2.fastq)

2 bam files (aln.bam and aln.bam.bai)

I know that:

• the raw file is aligned to hg19 human genome.
• The sequenced file is pair-ended from Miseq sequencing platform.
• This file is the result of an amplicon sequencing design.

I also have sequence information about the primers (4 forward and 4 reverse) and 2 adaptors.

And I have to answer these two questions:

1. What is the constitution of each read, for example (adaptor + primer + amplified region) ?
2. Which gene region did those reads mapped to ?

My question is what kind of tool do I have to use to respond to these questions? samtools? sequencher? FastQC? Can anybody help me just to start because I am totally lost.

Thanks

Answer 1

Samtools is not an aligner; it is used to analyse alignments. FastQC is for analysing the read qualities and average composition. You can get to know the adaptor sequence by running FastQC.

The next step is adaptor trimming. There are many tools that do that. Trimmomatic is currently popular. You can find the adaptor composition of each read by subtracting the read length post trimming from the read length pre-trimming (this would be same for all reads).

Then you have to align your trimmed reads to a reference, using an aligner. Again, there are many aligners available; Bowtie, STAR and now a new one called HISAT. STAR and HISAT are faster than Bowtie.

Since you already have the alignments (aln.bam) to hg19, you do not need to perform alignment. You want to know, to which genes do these reads map to. This information is not available unless you have the genome annotations. The annotations are available in the form of GTF/GFF files. You can obtain the association between genomic locations and the reads using the BEDOPS toolkit as mentioned in this Biostars post.

I don't have much experience with this tool. What I would do is convert the BAM to SAM (using samtools) or to BED (using bedtools:bamtobed) and parse the columns of the SAM/BED file describing the read locations with the GTF file. See this link for the details on the SAM format. You can use any scripting language like awk or perl to parse.

Samtools does have an option to filter the reads according to regions specified in the BED format but it will not automatically annotate them.

samtools view -hL regions.bed aln.bam > aln_filtered.sam

Answer 2

The aln.bam file is likely a binary and compressed version of a file in SAM format, indicating where reads have aligned. I guess it has been generated by some alignment program using the two fastq files.

For a start, you might use a program to visualize where the reads have mapped.

For instance, you may use IGV and load your aln.bam file to see where in hg19 they aligned. IGV will also use the aln.bam.bai file to help it find the aligned reads in the aln.bam file.

As far as I know, sequencher is used to clean and assemble Sanger sequencing data, but my experience with that program is almost 10 years old.

Samtools is a command-line program that can help you extract information from the aln.bam file. You may use it to count how many alignments are present in a given genomic region, for instance, using option -c of samtools view. See the manual for more details: http://www.htslib.org/doc/samtools.html