Short read alignment#
In the last session you learned about aligning sequences with blast. But have you ever tried to align 1 million reads with blast? Blast is a nice tool and an even nicer web interface to quickly query search for hits in big databases. But its unsuitable for alignment of large number of reads. For that there are other tools such as bwa for dna-vs-dna alignment or star for rna-vs-dna alignment. We will align the reads we generated in the previous step against a database and inspect results.
1 $ ls /science/teaching/ecoli_wgs/database
2 ecoli.fasta ecoli.fasta.amb ecoli.fasta.ann ecoli.fasta.bwt ecoli.fasta.pac ecoli.fasta.sa
The database can be found in the database folder. It consists of a fasta file and the tool specific database files.
We want to align reads against this database using the tool bwa.
1$ cd /science/teaching/ecoli_wgs/backmapping
2$ ml BWA/0.7.17-foss-2018b
3$ ml SAMtools/1.9-foss-2018b
4$ bwa mem /science/teaching/ecoli_wgs/database/ecoli.fasta /science/teaching/ecoli_wgs/qc/SRR11669351.noadapters.nophix.qualitytrimmed_1.fastq.gz /science/teaching/ecoli_wgs/qc/SRR11669351.noadapters.nophix.qualitytrimmed_2.fastq.gz > /science/teaching/ecoli_wgs/backmapping/SRR11669351.sam
bwa takes as an input a database and either one (single mode) or two read files (paired mode) and produces a [SAM](https://samtools.github.io/hts-specs/SAMv1.pdf) file as output. The SAM file is an alignment format that contains information about the input read, the reference it aligns to and the position where it aligned to. SAM (and the binary version BAM) files can be inspected with SAMtools. E.g. we want to see a couple alignments that bwa produced:
1 $ samtools view SRR11669351.sam | head -n 20
2
3 SRR11669351.1 99 NODE_18_length_88469_cov_139.157967 50064 60 150M = 50298 384 CGCGTAATGTGCAGAAACGCTTTGCCGAATGGGGCATAAAACAAGGGCTGAATAATCACGTTCATCCGCATAAATTACGTCACTCGTTCGCCACGCATATGCTGGAGTCGAGCGGCGATCTTCGTGGTGTGCAGGAGCTGCTGGGTCATG FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:150 MC:Z:150M AS:i:150 XS:i:0
4 SRR11669351.1 147 NODE_18_length_88469_cov_139.157967 50298 60 150M = 50064 -384 AACGGGGGAAATAATGCGTTTTTACCGGCCTTTGGGGCGCATCTCGGCGCTCACCTTTGACCTGGATGATACCCTTTACGATAACCGTCCGGTGATTTTGCGCACCGAGCGAGAGGCGCTTACCTTTGTGCAAAATTATCATCCGGCGCT :FFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF,FFFFFF:FFFFFFFFFFFF,:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:150 MC:Z:150M AS:i:150 XS:i:0
5 SRR11669351.2 99 NODE_18_length_88469_cov_139.157967 66191 60 150M = 66440 399 TGAATCACACCCTCGGTTTCCCTCGCGTTGGCCTGCGTCGCGAGCTGAAAAAAGCGCAAGAAAGTTATTGGGCGGGGAACTCCACGCGTGAAGAACTGCTGGCGGTAGGGCGTGAATTGCGTGCTCGTCACTGGGATCAACAAAAGCAAG FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:150 MC:Z:150M AS:i:150 XS:i:0
6 SRR11669351.2 147 NODE_18_length_88469_cov_139.157967 66440 60 150M = 66191 -399 AAGATGGTTCGGTAGATATCGACACCCTGTTCCGTATTGGTCGTGGACGTGCGCCGACTGGCGAACCTGCGGCGGCAGCGGAAATGACCAAATGGTTTAACACCAACTATCACTACATGGTGCCGGAGTTCGTTAAAGGCCAACAGTTCA FF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:150 MC:Z:150M AS:i:150 XS:i:0
7 SRR11669351.3 99 NODE_3_length_283843_cov_128.512507 8282 60 150M = 8401 269 AGTCACGTAGCGGCGGCAGATTGAGCGTCAACACGTTCAGACGATAATAGAGATCTTCACGGAACATGCCTTTTTGCACCAGTTCGACCAGATTCTTCTGCGTAGCGCAAATCACCCGCACATCGACATGCACCTCATGGTCTTCGCCAA FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:150 MC:Z:150M AS:i:150 XS:i:0
8 SRR11669351.3 147 NODE_3_length_283843_cov_128.512507 8401 60 150M = 8282 -269 ACATCGACATGCACCTCATGGTCTTCGCCAACCCGACGGAAAGTGCCATCATTAAGGAAACGCAGTAATTTCGCCTGCATCCGTGGTGACATTTCCCCTATTTCATCCAACAGCACCGAACCACCGTTCGCCTGCTCAAAGAATCCTTTC FFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:150 MC:Z:150M AS:i:150 XS:i:0
9 SRR11669351.5 83 NODE_4_length_268175_cov_132.594797 1897 60 150M = 1570 -477 AAGTATGTGGAAAACTGGCTGTTGTGGGTGATTATTAACGTGATTAGCGTCGTTATTTTTGCACTTCAGGGCGTTTACGCCATGTCTCTGGAGTACATCATCCTGACCTTTATTGCGCTCAACGGCAGCCGGATGTGGATCAACAGCGCA FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:150 MC:Z:150M AS:i:150 XS:i:0
10 SRR11669351.5 163 NODE_4_length_268175_cov_132.594797 1570 60 150M = 1897 477 CTGCTATTACAGGTGTTTTTCTTTGCCGCGAATATTTACGGTTGGTATGCGTGGTCGCGACAAACCAGTCAGAACGAGGCGGAGTTGAAAATTCGCTGGTTGCCATTGCCGAAGGCACTCAGCTGGTTGGCGGTTTGCGTTGTTTCGATT FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF,FFF,FF NM:i:0 MD:Z:150 MC:Z:150M AS:i:150 XS:i:0
11 SRR11669351.6 99 NODE_5_length_224476_cov_123.202020 159349 60 150M = 159534 335 AAGGGCTACAACTCATTCAGTATGTACTGGAACCGGCAAGTAGAATCGATTGGCTGCATACTGCTGGTATCAACACGCACGACCATATCGTGTAACTGTTTGTCCGTTATAGGGTGTCGATTTGGCGTGACCGAGACTTCCGCCTTGTCA FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:150 MC:Z:150M AS:i:150 XS:i:0
12 SRR11669351.6 147 NODE_5_length_224476_cov_123.202020 159534 60 150M = 159349 -335 CCAGCGGCAGTAACCACCTGTTTCATTCATCTGCGTACTGTCGCACAACTTCCCATCTTTCATCAAATACGTTGAAAGCGTTTTCTCAACCACACCGCCCGTCGACTGCAGGGTGAGCAGTTTTGACTTACTGACATAATTCTCATTCCA FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF::FFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF NM:i:0 MD:Z:150 MC:Z:150M AS:i:150 XS:i:0
13 SRR11669351.3011 77 * 0 0 * * 0 0 ACACCTTGCTTCATAGTAACCAGTTGAGATGAAGCACGTCGTTAGAACGTTGTTGGACACCATGTTGAATTATTCCCCCATCGGTTGTGAAGAACTGTGCTACATTCAGGCTTACCCATTGAACTCAGTATATATATTTTTTTCCTTCCT FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFF,,FFF:F:FFF:FFFFFFFFFFFFFFF:F AS:i:0 XS:i:0
14 SRR11669351.3011 141 * 0 0 * * 0 0 CAAGAATGGTATCCTGCCAGACAAAAGACAGGAAGGAAAAAAATATATATACTGAGTTCAATGGGTAAGCCTGAATGTAGCACAGTTCTTCACAACCGATGGGGGAATAATTCAACATGGTGTCCAACAACGTTCTAACGACGTGCTTCA FFFFFFFFFFFFFFFFFFFF:FFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:,FFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF AS:i:0 XS:i:0
The SAM file can be read the following. Insert SRR11669351.1 aligns fulllength (150M + 150M) against the reference NODE_18_length_88469_cov_139.157967. The alignment of the R1 sequence starts at position 50064 `. The alignment of the R2 sequence starts at position `50298.
Exercises#
Go through the SAM specification and try to understand the fields in a SAM file. Ask questions if you have problems to understand.
Insert SRR11669351.4 is missing in the SAM file. Where did it go?
How do you find unaligned inserts/reads?
Download a random e.coli genome from NCBI, build an index (bwa index ecoli.fasta) and align. Explain differences, especially the NM:i: field and in the CIGAR string (column 6).