Module 3 -- Mapping and Genome Rearrangement LAB #1 OBJECTIVE: Map FASTQ files to the human reference genome and explore the tools available for working with alignments. TOOLS USED: We will be using: * bwa (read mapper) * samtools (utilities to view/modify alignment files) * Picard (utilities to view/modify alignment files) * IGV (genome browser) Instructions: In the following exercise, lines starting with a # symbol are comments. The lines without # symbols are the commands that you should run in the terminal. ################### # Start of Exercise ################### # # Step 0. Create a directory to store your working data # Run the following commands by copying and pasting or typing them into the terminal: # cd ~/workspace mkdir module3 cd module3 # # Step 1. Map reads to the reference genome using bwa # # bwa mem is the latest mapping algorithm in the bwa software package # You can read about how it works here: http://arxiv.org/abs/1303.3997 # This command may take a few minutes to run bwa mem -t 4 ~/CourseData/CG_data/Module3/human_g1k_v37.fasta ~/CourseData/CG_data/Module3/reads_1.fastq ~/CourseData/CG_data/Module3/reads_2.fastq > sample.sam # If this command is taking a long time, you can just copy the results to your working directory # to proceed with the step 2 of the exercise: cp ~/CourseData/CG_data/Module3/sample.sam . # # Step 2. Explore the alignment data # # SAM is a plain-text format that can be viewed from the command line # Use the head command to look at your alignments. # You will see the SAM header containing the reference information # followed by a few alignments. head -100 sample.sam # In this SAM file, the reads are ordered by their position in the original FASTQ file. # We usually want to order the alignments by their position on the reference genome. # We will use Picard's SortSam tool to do this. The output will be in BAM format # which takes up much less disk space than SAM. java -jar /usr/local/picard/picard.jar SortSam I=sample.sam O=sample.sorted.bam SO=coordinate CREATE_INDEX=true # The samtools view command can be used to convert BAM to SAM (and vice-versa) samtools view sample.sorted.bam | head -100 # 'samtools view' can also extract the alignments for regions of the genome # this command will give you all the alignments for a 10kb region of chromosome 20 samtools view sample.sorted.bam 20:26,000,000-26,010,000 # As a tab-delimited file, the SAM format is easy to manipulate with common # unix tools like grep, awk, cut, sort. For example, this command # uses cut to extract just the reference coordinates from the alignments samtools view sample.sorted.bam 20:26,000,000-26,010,000 | cut -f3-4 # We can samtools flagstat command to get a summary of the alignments samtools flagstat sample.sorted.bam # We can use 'samtools idxstats' to count the number of reads mapped to each chromosome samtools idxstats sample.sorted.bam # This command will count the number of reads on chromosome 20 samtools idxstats sample.sorted.bam | awk '$1 == "20"' # # Step 3. Examining alignments # # You can view all of the aligned reads for a particular reference base using mpileup. # The following command will show the read bases and their quality scores at a heterozygous SNP. # The "." and "," symbols indicate bases that match the reference. # There are 18 reads that show a "G" base at this position. # The individual's genotype at this position is likely A/G. samtools mpileup -f ~/CourseData/CG_data/Module3/human_g1k_v37.fasta -r 20:32,001,292-32,001,292 sample.sorted.bam # Load the data into IGV by performing the following: # Open IGV and change the genome from hg19 to 'human_g1k_v37' # Choose 'Load from URL' from the file menu # Type: http://cbw#.entrydns.org/module3/sample.sorted.bam where # is your student ID # Navigate to 20:32,001,292 # Notice that the alignments have high mapping quality. # Now we will view the alignments for a different position: samtools mpileup -f ~/CourseData/CG_data/Module3/human_g1k_v37.fasta -r 20:25,997,273-25,997,273 sample.sorted.bam # In this case, 11 reads show a T base at this position. # Look at the alignments in IGV by navigating to 20:25,997,273. # The reads colored white have mapping quality 0. # This means the alignment is ambiguous and should not be trusted. # It is unclear whether the T->C alignments are true SNPs. # This is the end of lab 1. You can use the remaining time to explore the alignments # and ask questions if you notice anything unusual or interesting.