The Journey of each DNA Fragment from Organism to Sequencing Read

In Section 7.2.5 we discussed one solution for downloading gzipped FASTQ files from a sequencing center. This turns out to be the last step in a very long process of acquiring sequencing data from a collection of organisms that you are studying. We could say that the DNA of each organism that you are studying took an extraordinary journey to get to you. This journey included, 1) you (or a colleague) sampling tissue or blood from that organism, 2) DNA being extracted from that tissue or blood, 3) that DNA getting Illumina adapters (with sample-specific barcodes) attached to it and then getting PCR amplified some number of times in the library preparation stage, 4) the products of the library preparations are sequenced, sometimes on multiple flow cells and often across multiple lanes, during the sequencing stage, 5) the sequencer processes the read data into FASTQ files, 6) finally, you download the data.

Anyone who has ever been involved in research knows that at every stage of the above process there are opportunities for errors or biases to occur that could affect the final data (and inferences) produced. The effects of some of the errors and biases that occur during the library prep stage and the sequencing stage can be eliminated or lessened. Doing so requires that we keep careful track of the route that each read took during those those steps in its journey. This is the role of read groups, which we discuss below. For now, however, we will offer that the main errors/effects that can be controlled for by careful accounting of reads into read groups are: 1) controlling for effects on the reported base quality scores of the particular flowcell or the particular lane on a flowcell that a read was sequenced on; 2) The identification of sequencing reads which are PCR duplicates of another read in the data set. Thus, most of the focus on defining read groups is concerned with identifying features of a read’s journey that can be used to ameliorate those two types of errors/biases.

The process of controlling for batch effects of flowcells and lanes on base quality scores is called, unsurprisingly, Base Quality Score recalibration. This might be discussed in a later section. It is a very important process when one is searching for de novo mutations, or otherwise attempting to very accurately identify variants that occur at very low frequency. Doing base quality score recalibration is considered a “best practice” when using GATK (Section XX), but I have my doubts as to whether it will make appreciable differences in analyses where the importance of rare variants is downweighted. Nonetheless, whether you are going to do base quality score recalibration or not, you should certainly prepare your data so that you can do it should you decide to. In order to do it, you must know which flowcell and upon which lane each read was sequenced. That information can be found in the FASTQ file as part of the name for each sequence, however, methods for doing base quality score recalibration do not typically use the names of the reads to access that information. Rather it must be stored in the read group information.

A PCR duplicate is a template sequence that happens to be a copy (made by PCR) of another sequence that has also been sequenced. They are considered “errors” in bioinformatics, because many of the models for inferring genotypes from next generation sequencing data assume a process of independent sampling of gene copies (i.e. the maternal and paternal copies of a choromosome) within a diploid individual, and that model and the data are used to estimate how likely an individual is a homozygote versus a heterozygote at each position. For example, if Individual X produces one read covering position Y—a position known to be a variable SNP with alleles A and G in the species—and the read carries an A at position Y, then you would conclude that Individual X likely has genotype AA or AG at position Y, and is less likely to have genotype GG, because our observed read with an A at that position would have only occurred because of a sequencing error in that case. If we saw 10 reads from individual X, and they all had the A base at position Y, however, we could be quite confident that the genotype is AA rather than AG, because if the indvidual had genotype AG, then each read should have a 50% chance of carrying an A and a 50% chance of carrying a G, and it is very unlikely, in that case, that all 10 would carry the A allele. However, imagine how your conclusions would change if someone were to tell you that 95% of the reads were PCR duplicates of another one of the reads. In that case, even if the individual has genotype AG, there is high probability that only a single, true template carrying an A came from the individual, and the remaining 9 reads that carry an A are just PCR copies of that first original template.

PCR duplicates are copies of the same template fragments. Therefore in paired-end data if a pair of reads maps to exactly the same location as another pair of reads, it is likely that one of them is a PCR duplicate of the other. Figure 19.1 provides a simplified schematic (with merely single-end data) describing a situation with high and low numbers of PCR duplicates. The PCR duplicate fragments are denoted with gray-colored flanking sequence. In the figure, note that the length of the fragments gives information about which fragments are duplicated and which are not—if two fragments of the same length are identified, then one of them is considered a PCR duplicate. It is worth noting that it isn’t (to my knowledge) really known which fragment is the duplicate and which is the “original,” but, in order to correct the situtation, one (or more of them) will be identified as duplicates, and one of them will be identified as an “original.”

 An example of PCR duplicates. Note that in paired end data, duplicates are identified by the lengths of both reads as well as the insert length. This figure represents a simplified, toy, situation, showing just single-end data.

PCR duplicates occur during the library prep stage, and also are only possible for fragments originating from the same individual sample. In other words, if you found two different paired-end reads that mapped to the same location in the genome, but they were from different individuals (i.e. had different sample-specific DNA barcodes), you would not suspect that they were PCR duplicates. Likewise, DNA fragments from two different library preps, even if they came from the same individual—in the case that DNA from that individual had been prepared in two separate libraries—would not be considered to be PCR duplicates. (We should pause here to note what is considered a “library prep.” Basically, if a batch of samples were all processed together in the same day with reagents from the same kit, that is considered a single “library prep.”) Accordingly, when keeping track of read groups for the purposes of identifying PCR duplicates, the most important information is which individuals the DNA came from and which library preparation it was prepared in.

What we can conclude from this section is that, for downstream analyses, we need a way to quickly (and efficiently, in terms of data storage space) identify, for each read the:

• Sample it originated from
• The library in which it was prepared for sequencing.
• The flowcell that it was sequenced on.
• The lane on that flowcell that it was sequenced on

This information can be strored in the alignment file using read group information as described below.

Read Groups

The SAM file header can contain lines like this one:

@RG ID:HTYYCBBXX.1.CH_plate1_A01    SM:CH_plate1_A01    LB:Lib-1    PU:HTYYCBBXX.1.TCAATCCG+TTCGCAGT    PL:ILLUMINA

The purpose of this line is to identify a read group that is named via the ID tag. In the above case the read group is named HTYYCBBXX.1.CH_plate1_A01. The attributes of this read group are stored using the remaining tags in the @RG line. Happily, those attributes mirror exactly what we might need to know to identify PCR duplicates and to perform base quality score recalibration; namely:

• SM : give the name of the sample, CH_plate1_A01, the read is from,
• LB : gives the name of the library prep, Lib-1, the read is from,
• PU : gives the name of the flowcell, HTYYCBBXX, and the lane, 1, the read is from.

It is important to note that this line in the SAM header is just a header line that gives information about a particular read group that is found in the file. It is not, by itself giving you information about any particular read found in the alignment file. Read group information about a particular read in the alignment file is given on the alignment line in a column (way off near the end of the line typically) tagged with RG:Z. That RG:Z column in an alignment row simply holds the ID of the read group that the read belongs to. In this way, all the information that applies to the read—i.e., the sample it came from, (SM), the library prep it was in (LB), and information about the flowcell and lane (PU)—can be assigned to the read, without writing all of that information on every line in the alignment file. Accordingly, it should be obvious why the read group IDs must be unique across all the different read groups (@RG tags) in a SAM/BAM file: the ID is used to identify all the other information that goes along with that read, and is used to differentiate those different read groups!

You can download the first 100 lines (without the @SQ lines) of a SAM file that has 8 different read groups in it here. Let’s do that and have a look at it.

The way that read group information—particularly the platform unit, or PU information—is recorded is dictated by which software will be used for downstream processing. We describe the conventions that are useful for passing the aligned data to the Genome Analysis Toolkit (GATK), an extensive and well-supported piece of software developed by the Broad Institute. Currently, the GATK conventions for read group tags are these:

• ID : any string that uniquely identifies a read group. The user must make sure that these uniquely identify the read groups. If you have multiple samples, libraries, flowcells, and lanes, and easy way to ensure this is to just bung those all together in a format like {Sample}.{Library}.{Flowcell}.{Lane}, where the parts in curly braces should be replaced by the actual names of those items. As long as there is a PU field, however, you can make the ID anything you want, as long as it is unique for each read group.
• SM : any string (no spaces please!) that serves as an identifier for the individual the sample was taken from. This is used to assign genotypes to indvididuals, so even if you have different “samples” from the same individual (i.e. some blood and some tissue, prepped in different libraries; or an original extraction prepared in an earlier library and a later extraction prepared in a library for individuals that needed more sequening coverage after an initial sequencing run) you should still assign them all the same SM tag because they are the same individual, and you will’ want all of the DNA sequences from them to apply to calling genotypes for them.
• LB : any string which gives the name of the library prep these reads were prepared in. If you prepped 192 individuals in one day, they are all part of the same library.
• PU : GATK expects this information to be given in a specific format. You can provide it as {Flowcell}.{Lane} (for example, HTYYCBBXX.1), but with later versions of the base quality score recalibration methods in GATK, it appears that it might be useful to also add information about the particular individual/library the read came from. It is suggested that this be done by including the barcode of reads, as {Flowcell}.{Lane}.{Barcode}, like HTYYCBBXX.1.TCAATCCG+TTCGCAGT.
• PL : It is also important for base calibration to indicate what type of sequencer was used. Most of the time, currently, that will be ILLUMINA.

Usually, when you get data back from the sequencing center, it will be in a series of different FASTQ files, each one containing reads from a single DNA sample well (on a plate that you sent to them), that was sequenced on a certain lane on a given flow cell. The files will often be named something like: DPCh_plate1_C07_S31_L7_R2.fq.gz where the first part gives a sample identifier, the part after the L gives the lane, and the part after the R tells whether the file holds read 1 or read 2 of paired end reads. You can get the flowcell name by looking at the sequence names inside the file (See Section 17.2.1). You will map the reads in each file separately, and when you do so, you will attach the read group information to them (see below). Accordingly, it is a great idea to maintain a spreadsheet that links the different files to the attributes and the read group information of the reads inside them, like this:

index   file_prefix ID  PU  SM  PL  LB  Flowcell    Lane
1   DPCh_plate1_A01_S1_L1_R HTYYCBBXX.1.CH_plate1_A01   HTYYCBBXX.1.TCAATCCG+TTCGCAGT   CH_plate1_A01   ILLUMINA    Lib-1   HTYYCBBXX   1
2   DPCh_plate1_A01_S1_L2_R HTYYCBBXX.2.CH_plate1_A01   HTYYCBBXX.2.TCAATCCG+TTCGCAGT   CH_plate1_A01   ILLUMINA    Lib-1   HTYYCBBXX   2
3   DPCh_plate1_A01_S1_L3_R HTYYCBBXX.3.CH_plate1_A01   HTYYCBBXX.3.TCAATCCG+TTCGCAGT   CH_plate1_A01   ILLUMINA    Lib-1   HTYYCBBXX   3
4   DPCh_plate1_A01_S1_L4_R HTYYCBBXX.4.CH_plate1_A01   HTYYCBBXX.4.TCAATCCG+TTCGCAGT   CH_plate1_A01   ILLUMINA    Lib-1   HTYYCBBXX   4
5   DPCh_plate1_A01_S1_L5_R HTYYCBBXX.5.CH_plate1_A01   HTYYCBBXX.5.TCAATCCG+TTCGCAGT   CH_plate1_A01   ILLUMINA    Lib-1   HTYYCBBXX   5
6   DPCh_plate1_A01_S1_L6_R HTYYCBBXX.6.CH_plate1_A01   HTYYCBBXX.6.TCAATCCG+TTCGCAGT   CH_plate1_A01   ILLUMINA    Lib-1   HTYYCBBXX   6
7   DPCh_plate1_A01_S1_L7_R HTYYCBBXX.7.CH_plate1_A01   HTYYCBBXX.7.TCAATCCG+TTCGCAGT   CH_plate1_A01   ILLUMINA    Lib-1   HTYYCBBXX   7
8   DPCh_plate1_A01_S1_L8_R HTYYCBBXX.8.CH_plate1_A01   HTYYCBBXX.8.TCAATCCG+TTCGCAGT   CH_plate1_A01   ILLUMINA    Lib-1   HTYYCBBXX   8
9   DPCh_plate1_A02_S2_L1_R HTYYCBBXX.1.CH_plate1_A02   HTYYCBBXX.1.CGCTACAT+CGAGACTA   CH_plate1_A02   ILLUMINA    Lib-1   HTYYCBBXX   1

In the exercise we will see how to use such a file to assign read groups when mapping.

Aligning reads with bwa

There are several programs for aligning reads to a reference genome. We focus on bwa which is an industry standard aligner written by Heng Li and Richard Durbin (Li and Durbin 2009). The name stands for “Burrows-Wheeler Aligner.” We won’t go deeply into the guts of this alignment algorithm, but we will briefly state that nearly all alignment methods rely on pre-processing the reference genome into a data-structure (like a suffix tree) that provides an index which makes it fast to find the place in a genome where a query sequence matches. Such indexes can take up a large amount of computer memory. The Burrows-Wheeler Transform (BWT) provides a way of decreasing the size of such indexes. The BWT is an operation that takes a sequence of characters (in this case DNA bases) and re-orders them so that similar characters tend to appear together in long runs of the same character. Sequences with long runs of the same character can be easily compressed to take up less space using run length encoding, for example. (We have already seen an example of run length encoding in the CIGAR string , in which runs of the same kind of alignment characteristic (Matches, Deletions, etc) where encoded by their type and the length of the run.) The remarkable thing about the BWT is that it is invertible: if someone gives you a sequence and says that it is the result of applying the BWT to an original sequence, you can actually recover the original sequence without any other information! The program bwa uses the BWT to compress the index used for alignment so that it can easily fit into memory—even for very large genomes—on most any computer system, even a low-powered laptop produced in the last half decade.

bwa index path-to-genome

bwa index genome/GCA_002872995.1_Otsh_v1.0_genomic.fna.gz

% ls -hl
total 3.8G
-rw-rw-r-- 1 eanderson eanderson 704M Mar  5 00:08 GCA_002872995.1_Otsh_v1.0_genomic.fna.gz
-rw-rw-r-- 1 eanderson eanderson 4.4M Mar  5 00:47 GCA_002872995.1_Otsh_v1.0_genomic.fna.gz.amb
-rw-rw-r-- 1 eanderson eanderson 2.1M Mar  5 00:47 GCA_002872995.1_Otsh_v1.0_genomic.fna.gz.ann
-rw-rw-r-- 1 eanderson eanderson 2.3G Mar  5 00:47 GCA_002872995.1_Otsh_v1.0_genomic.fna.gz.bwt
-rw-rw-r-- 1 eanderson eanderson 579M Mar  5 00:47 GCA_002872995.1_Otsh_v1.0_genomic.fna.gz.pac
-rw-rw-r-- 1 eanderson eanderson 1.2G Mar  5 00:58 GCA_002872995.1_Otsh_v1.0_genomic.fna.gz.sa

bwa mem   path-to-reference-genome   path-to-read1-fastq-file   path-to-read2-fastq-file

bwa mem genome/GCA_002872995.1_Otsh_v1.0_genomic.fna.gz fastqs/DPCh_plate1_A01_S1_L1_R1.fq.gz fastqs/DPCh_plate1_A01_S1_L1_R2.fq.gz

ID=HTYYCBBXX.1.CH_plate1_A01
SM=CH_plate1_A01
LB=Lib-1
PU=HTYYCBBXX.1.TCAATCCG+TTCGCAGT
PL=ILLUMINA

bwa mem -R "@RG\tID:$ID\tSM:$SM\tLB:$LB\tPU:$PU\tPL:$PL"  genome.fasta  file_R1.fq file_R2.fq

19.4 Processing alignment output with samtools

Doing the alignment with bwa mem is only the first step of getting data ready to do variant calling. Once the SAM format data comes out of bwa mem, a number of things must happen to it before it can be used for variant calling:

1. It must be converted to BAM format (the compressed binary companion format to SAM)
2. Additional information about mate pairs might need to be attached to sequences
3. You must sort the BAM files from the original ordering produced by bwa mem into a coordinate-ordered format, in which the reads are sorted by where they appear within the reference genome.
4. If you started with multiple pairs of FASTQ files for a single individual (i.e. from different lanes or from different sequencing runs or libraries) you might need to merge those all into a single BAM file before variant calling.
5. You might want to mark PCR duplicates as such in the BAM file.
6. Finally when you have your BAM file all ready for variant calling, you typically will need to index it for rapid access by the variant caller program.

Phew! That is a lot of steps! Fortunately, there is a single “go-to” program (one-stop shopping!) that can manage all of that for you. It is the program samtools, brought to you by the same people who developed the SAM/BAM formats (Li et al. 2009) and the bwa aligner (and bcftools for that matter). If you want to get anywhere in bioinformatics, it is important to become good friends with samtools.

The program samtools includes a large number of different subcommands. We will cover just a few of them, here, that are used in our typical paired-end workflow for whole genome sequencing, and in a somewhat cursory fashion, and we will disscuss only the most commonly used options. All readers are encouraged to completely read the online manual page for samtools for a complete account of its uses.

Or, you can read the usage guidelines for any subcommand from the command line:

conda activate bwasam

# This lists all the different subcommands for you, categorized by
# their uses
samtools  

# if you follow samtools by any subcommand with nothing else, it gives you information
# about all the options of that subcommand.  For example:
samtools view

samtools sort

The usage patterns of the subcommands can broadly be categorized into two groups:

1. Subcommands that explicitly require an argument that gives the name of an output file. These subcommands are usually described like this:

   samtools subcommand [options] out.bam in.bam

   where an output file is explicitly called for. These commands do not write output to stdout! They must be given an output file name. It is not possible to pipe output from one of these subcommands to another program or utility.

2. Subcommands that don’t explicitly require an argument giving an output file:

   samtools subcommand [options] in.bam

   These subcommands write their output to stdout, and that output can be piped into other programs or can be redirected into a file (some also provide a -o file option to write output to file rather than to stdout.)

In short, if the syntax of the subcommand explicitly calls for an output file on the command line (not as part of a -o file option), then that subcommand does not write its output to stdout and you can’t pipe it into the next part of your pipeline.

On the other hand, if you can write output from samtools to stdout, you might be interested in knowing how you can pipe that into the input for another samtools subcommand. Rejoice! There is a syntactically economical and lovely way to do that. Any subcommand which expects a single input file will be described like:

samtools subcommand [options]  in.bam

You would typically put the path to the input file in the place of in.bam there. However, if you put - in place of in.bam then samtools will take input for in.bam from stdin.

So, for example, you can do something like this:

bwa mem genome.fna R1.fq R2.fq |
  samtools view -u -  |   # convert the SAM output from bwa mem into BAM format
  samtools sort -l 9  - -o output_file.bam  # take stdin as the input, sort it, and write (with the best
                                        # compression possible: -l 9) the output to output_file.bam

This is HUGE! It means that you don’t even have to save the initial SAM file that bwa mem makes, you can proceed directly to the sorted BAM file you want without eating up disk space on the intermediate files.

mkdir -p results/sam
wget https://eriqande.github.io/eca-bioinf-handbook/downloads/s001---1.sam
mv s001---1.sam results/sam/

mamba create -n bwasam -c bioconda bwa=0.7.17 samtools=1.11

ssh scompile
srun --partition shas-interactive --pty /bin/bash

conda activate bwasam

samtools view

samtools view -b results/sam/s001---1.sam > results/sam/s001---1.bam

du -h results/sam/*

head results/sam/s001---1.bam

# pipe output to head to look at the first 10 alignments
samtools view results/sam/s001---1.bam | head

samtools view -h results/sam/s001---1.bam | head

# this is a handy way to get the last N (12 in this case) lines of the
# SAM header:
samtools view -H results/sam/s001---1.bam | tail -n 12

samtools sort

samtools sort -l 9 -o results/sam/s001---1.srt.bam results/sam/s001---1.bam

du -h results/sam/*

samtools view -H results/sam/s001---1.srt.bam | head

samtools view -H results/sam/s001---1.srt.bam | tail -n 10

# remove the .bam .srt.bam files:
rm results/sam/s001---1.bam results/sam/s001---1.srt.bam

# now, directly make a sorted BAM file at results/sam/s001---1.srt.bam
# by piping samtools view output into samtools sort.  Note the use
# of -u for uncompressed BAM output, and the - at the end of the
# line, instead of a file name, to mean take
# input from stdin instead of a file
samtools view -u results/sam/s001---1.sam | \
  samtools sort -l 9 -o results/sam/s001---1.srt.bam -

samtools index results/sam/s001---1.srt.bam

du -h results/sam/*

samtools view results/sam/s001---1.srt.bam CM031202.1:1001-1100

samtools merge [options] output-bam-name.bam  sorted-input-bam-1.bam sorted-input-bam-2.bam ...

samtools flagstats results/sam/s001---1.srt.bam

samtools idxstats results/sam/s001---1.srt.bam | less

samtools stats results/sam/s001---1.srt.bam | less

 BONEYARD BELOW HERE

Preprocess ?

Will have a bit about sequence pre-processing (with WGS data it already comes demultiplexed, so maybe we can hold off on this until we get to RAD data). No, we need to talk about trimming and maybe slicing. Perhaps put that in a separate chapter of “preliminaries”

Quick notes to self on chaining things:

When using samtools markdup you can do like this:

# first step
bwa mem chinook-play/chinook-genome-idx/GCA_002872995.1_Otsh_v1.0_genomic.fna.gz chr32-160-chinook/laned/DPCh_plate1_A01_S1_L1_R1.fq chr32-160-chinook/laned/DPCh_plate1_A01_S1_L1_R2.fq | \
    samtools view -b -1 - | \
    samtools fixmate -m -O BAM - fixed.bam
    
# unfortunately, fixmate can't write to stdout.
# then, after this, we have to sort it into coordinate order and
# markdup them.
samtools sort fixed.bam | samtools markdup - marked.bam

# after that, fixed.bam could get tossed.

Merging BAM files

There is a lot of discussion on biostars about how samtools does not reconstruct the @RG dictionary. But I think that this must be a problem with an older version. The newer version works just fine. That said, Picard’s MergeSamFiles seems to be just about as fast (in fact, faster. For a comparably sized file it took 25 minutes, and gives informative output telling you where it is at). And samtools merge is at well over 30 and only about 3/4 of the way through. Ultimately it took 37 minutes. Might have been on a slower machine…

However, if you have sliced your fastqs and mapped each separately, then samtools let’s you not alter duplicate read group IDs, and so you can merge those all together faithfully, as I did in the impute project. Cool.