Read alignment - advanced

Learning outcomes

After having completed this chapter you will be able to:

• Use samtools to mark duplicates from an alignment file
• Use samtools to add readgroups to an alignment file
• Use a for loop in bash to perform the same operation on a range of files
• Use samtools in a pipe to efficiently do multiple operations on an alignment file in a single command

Material

• samtools documentation

Exercises

1. Adding readgroups

During several steps of variant calling gatk uses read group information. For each read, this gives information on the sequencing platform, the library, the lane and of course the sample. Have a look at the description of the different levels of read group information gatk uses here.

Exercise: The documentation mentions several read group fields that are used by gatk. Have a look at the fastq header. Does that give you the information that is required? Do we have that information for our sample? Can you specify it for our sample?

zcat mother_R1.fastq.gz | head

@EAS139:136:FC706VJ:2:2104:15343:197393 1:Y:18:ATCACG

Exercise: Have a look at the documentation of AddOrReplaceReadGroups. Specify the required arguments, and run the command. Do this from a script called B05_add_readgroups.sh (in ~/workdir/scripts/B-mother_only).

#!/usr/bin/env bash

cd ~/workdir/results

gatk AddOrReplaceReadGroups \
--INPUT alignments/mother.bam \
--OUTPUT alignments/mother.rg.bam \
--RGLB lib1 \
--RGPU H0164.2.ALXX140820 \
--RGPL ILLUMINA \
--RGSM mother \
--RGID H0164.2

Exercise: Compare the header and first alignments of mother.bam and mother.rg.bam. Notice any differences?

samtools view -H <alignment.bam>

samtools view <alignment.bam> | head

@RG     ID:H0164.2      LB:lib1 PL:ILLUMINA     SM:mother       PU:H0164.2.ALXX140820

2. Mark duplicates

Now that we have specified read groups, we can mark the duplicates with gatk MarkDuplicates.

Exercise: Have a look at the documentation, and run gatk MarkDuplicates with the three required arguments. Do this from a script called B06_mark_duplicates.sh (in ~/workdir/scripts/B-mother_only).

#!/usr/bin/env bash

cd ~/workdir/results

gatk MarkDuplicates \
--INPUT alignments/mother.rg.bam \
--OUTPUT alignments/mother.rg.md.bam \
--METRICS_FILE alignments/marked_dup_metrics_mother.txt 

Exercise: Run samtools flagstat on the alignment file with marked duplicates, and write the output to a file called mother.rg.md.bam.flagstat. Create a script called B07_get_alignment_stats_after_md.sh (in ~/workdir/scripts/B-mother_only). How many reads were marked as duplicate?

#!/usr/bin/env bash

cd ~/workdir/results/alignments

samtools flagstat mother.rg.md.bam > mother.rg.md.bam.flagstat

133477 + 0 in total (QC-passed reads + QC-failed reads)
0 + 0 secondary
317 + 0 supplementary
17329 + 0 duplicates
132892 + 0 mapped (99.56% : N/A)
133160 + 0 paired in sequencing
66580 + 0 read1
66580 + 0 read2
131470 + 0 properly paired (98.73% : N/A)
131990 + 0 with itself and mate mapped

585 + 0 singletons (0.44% : N/A)
0 + 0 with mate mapped to a different chr
0 + 0 with mate mapped to a different chr (mapQ>=5)

3. Indexing

To look up specific alignments, it is convenient to have your alignment file indexed. An indexing can be compared to a kind of ‘phonebook’ of your sequence alignment file. Indexing can be done with samtools as well, but it first needs to be sorted on coordinate (i.e. the alignment location). You can do it like this:

samtools index <bam file>

Exercise: Create a script called B08_index_alignment.sh (in ~/workdir/scripts/B-mother_only) to perform the alignment.

#!/usr/bin/env bash

cd ~/workdir/results/alignments/

samtools index mother.rg.md.bam

4. Recap - mother only

Now we have performed now the following steps on the sample mother:

• Read alignment
• Adding readgroups
• Marking of duplicates
• Indexing

Your scripts directory should look like this:

├── A-prepare_references
│   ├── A01_download_course_data.sh
│   └── A02_create_bwa_index.sh
├── B-mother_only
│   ├── B01_alignment.sh
│   ├── B02_get_alignment_statistics.sh
│   ├── B03_sort_alignment.sh
│   ├── B04_compress_alignment.sh
│   ├── B05_add_readgroups.sh
│   ├── B06_mark_duplicates.sh
│   ├── B07_get_alignment_stats_after_md.sh
│   └── B08_index_alignment.sh
└── C-all_samples

5. Apply it on all three samples with pipes and loops

We now apply these steps to all three samples. In order to do that, we combine the alignment with the sorting and compression in one command. We can do that with piping the output of bwa to samtools sort and samtools view, like this:

SAMPLE="mother"

bwa mem data/reference/Homo_sapiens.GRCh38.dna.chromosome.20.fa \
data/fastq/"$SAMPLE"_R1.fastq.gz \
data/fastq/"$SAMPLE"_R2.fastq.gz \
| samtools sort \
| samtools view -bh > results/alignments/"$SAMPLE".bam

Exercise: Make a directory in the scripts directory C-all_samples (so ~/workdir/scripts/C-all_samples). In here, create a script called C01_alignment_sorting_compression.sh. Within that script use the above snippet to make a loop that performs the alignment, sorting and compression for all three samples (i.e. mother, father and son).

scripts
├── A-prepare_references
│   ├── A01_download_course_data.sh
│   └── A02_create_bwa_index.sh
├── B-mother_only
│   ├── B01_alignment.sh
│   ├── B02_get_alignment_statistics.sh
│   ├── B03_sort_alignment.sh
│   ├── B04_compress_alignment.sh
│   ├── B05_add_readgroups.sh
│   ├── B06_mark_duplicates.sh
│   ├── B07_get_alignment_stats_after_md.sh
│   └── B08_index_alignment.sh
└── C-all_samples
    └── C01_alignment_sorting_compression.sh

#!/usr/bin/env bash

cd ~/workdir

for SAMPLE in mother father son
do
    bwa mem data/reference/Homo_sapiens.GRCh38.dna.chromosome.20.fa \
    data/fastq/"$SAMPLE"_R1.fastq.gz \
    data/fastq/"$SAMPLE"_R2.fastq.gz \
    | samtools sort \
    | samtools view -bh > results/alignments/"$SAMPLE".bam
done

Now we continue with adding the readgroups. For each sample, we have to add specific information to the different readgroup fields. We can do that by looping over a tab delimited file with sample-specific information in each row. Let’s create that tab-delimited file.

Exercise Generate a tab-delimited file called sample_rg_fields.txt and store it in ~/workdir/results/. In this file, each line should represent a sample (mother, father and son), and you specify the SM, LB, PU and ID fields. E.g., the first line (for ‘mother’) would look like:

mother  lib1    H0164.2.ALXX140820  H0164.2

mother  lib1    H0164.2.ALXX140820  H0164.2
father  lib2    H0164.3.ALXX140820  H0164.3
son lib3    H0164.6.ALXX140820  H0164.6

Exercise Generate a script called C02_add_readgroups.sh (in ~/workdir/scripts/C-all_samples) to loop over the tab-delimited file (have a look at the last exercise in Setup), and add the correct readgroups to the bam file of each sample with gatk AddOrReplaceReadGroups.

cd ~/workdir/results

cat sample_rg_fields.txt | while read SAMPLE LB PU ID
do
    echo $SAMPLE $LB $PU $ID
done

#!/usr/bin/env bash

cd ~/workdir/results

cat sample_rg_fields.txt | while read SAMPLE LB PU ID
do
    gatk AddOrReplaceReadGroups \
    --INPUT alignments/"$SAMPLE".bam \
    --OUTPUT alignments/"$SAMPLE".rg.bam \
    --RGLB "$LB" \
    --RGPU "$PU" \
    --RGPL ILLUMINA \
    --RGSM "$SAMPLE" \
    --RGID "$ID"
done 

As final step, we will mark the duplicates and perform the indexing for the three samples.

Exercise: Generate two scripts called C03_mark_duplicates.sh and C04_index_alignments.sh, in which you loop over the sample names and perform the respective calculations. You can use B06_mark_duplicates.sh and B08_index_alignment.sh as a template.

#!/usr/bin/env bash

cd ~/workdir/results

for SAMPLE in mother father son
do
    gatk MarkDuplicates \
    --INPUT alignments/"$SAMPLE".rg.bam \
    --OUTPUT alignments/"$SAMPLE".rg.md.bam \
    --METRICS_FILE alignments/marked_dup_metrics_"$SAMPLE".txt 
done

#!/usr/bin/env bash

cd ~/workdir/results

for SAMPLE in mother father son
do
    samtools index alignments/"$SAMPLE".rg.md.bam
done  

6. Recap - all samples

Now, we have performed on all three samples:

• Alignment
• Sorting
• Compression
• Adding readgroups
• Marking of duplicates
• Indexing

Your scripts directory should look like this:

scripts
├── A-prepare_references
│   ├── A01_download_course_data.sh
│   └── A02_create_bwa_index.sh
├── B-mother_only
│   ├── B01_alignment.sh
│   ├── B02_get_alignment_statistics.sh
│   ├── B03_sort_alignment.sh
│   ├── B04_compress_alignment.sh
│   ├── B05_add_readgroups.sh
│   ├── B06_mark_duplicates.sh
│   ├── B07_get_alignment_stats_after_md.sh
│   └── B08_index_alignment.sh
└── C-all_samples
    ├── C01_alignment_sorting_compression.sh
    ├── C02_add_readgroups.sh
    ├── C03_mark_duplicates.sh
    └── C04_index_alignment.sh

Questions & Answers