Here are the R packages you'll need. Some will install with install.packages(). The packages listed under Bioconductor need to be installed with the dedicated installer. That's described here. If you need to do anything further, installation will be described in the recipes in which the packages are used:

• Bioconductor: Following are the packages:
  • Biostrings
  • GenomicRanges
  • gmapR
  • karyoploteR
  • rtracklayer
  • systemPipeR
  • SummarizedExperiment
  • VariantAnnotation
  • VariantTools
• rrBLUP

Bioconductor is huge and has its own installation manager. You can install these packages with the following code (further information is available at https://www.bioconductor.org/install/):

if (!requireNamespace("BiocManager"))
    install.packages("BiocManager")
BiocManager::install()

Normally, in R, a user will load a library and use the functions directly...

A key bioinformatics task is to take an alignment of high-throughput sequence reads, typically stored in a BAM file, and compute a list of variant positions. Of course, this is ably handled by many external command-line programs and tools and usually results in a VCF file of variants, but there are some really powerful packages in Bioconductor that can do the whole thing, and in a fast and efficient manner, by taking advantage of BiocParallel's facilities for parallel evaluation—a set of tools designed to speed up work with large datasets in Bioconductor objects. Using Bioconductor tools allows us to keep all of our processing steps within R, and in this section, we'll go through a whole pipeline—from reads to lists of genes carrying variants—using purely R code and a number of Bioconductor...

A draft genome assembly of a previously unsequenced genome can be a rich source of biological knowledge, but when genomics resources such as gene annotations aren't available, it can be tricky to proceed. Here, we'll look at a first stage pipeline for finding potential genes and genomic loci of interest absolutely de novo and without information beyond the sequence. We'll use a very simple set of rules to find open reading frames—sequences that begin with a start codon and end with a stop codon. The tools for doing this are encapsulated within a single function in the Bioconductor package, systemPipeR. We'll end up with yet another GRanges object that we can integrate into processes downstream that allow us to cross-reference other data, such as RNAseq, as we saw in the Finding unannotated transcribed...

One of the most rewarding and insightful things we can do is visualize data. Very often, we want to see on a chromosome or genetic map where some features of interest lie in relation to others. These are sometimes called chromosome plots, and sometimes ideograms, and in this section, we'll see how to create one of these using the karyoploteR package. The package takes as input the familiar GRanges objects and creates detailed plots from configuration. We'll take a quick look at some different plot styles and some configuration options for ironing out the bumps in your plots when labels spill off the page or overlap each other.

For this recipe, you&apos...

In pipelines where we've called variants, we'll often want to do subsequent analysis steps that need further filtering or classification based on features of the individual variants, such as the depth of coverage in the alternative allele. This is best done from a VCF file, and a common protocol is to save a VCF of all variants from the actual calling step and then experiment with filtering that. In this section, we'll look at taking an input VCF and filtering it to retain variants in which the alternative allele is the major allele in the sample.

We'll need a tabix index VCF file; I provide one in the datasets/ch2/sample.vcf.gz file...

Very often, you'll want to look in more detail at data that falls in a particular genomic region of interest, whether that be the SNPs and variants in a gene or the genes in a particular locus. This extremely common task is handled very well by the extremely powerful GRanges and SummarizedExperiment objects, which are a little fiddly to set up but have very flexible subsetting operations that make the effort well worth it. We'll look at a few ways to set up these objects and a few ways we can manipulate them to get interesting information.

In this recipe, we need the GenomicRanges, SummarizedExperiment, and rtracklayer Bioconductor packages. We&apos...

A powerful application of being able to find many thousands of genetic variants in many samples using high-throughput sequencing is genome-wide association studies (GWAS) of genotype and phenotypes. GWAS is a genomic analysis set of genetic variants in different individuals or genetic lines to see whether any particular variant is associated with a trait. There are numerous techniques for doing this, but all rely on gathering data on variants in particular samples and working out each sample's genotype before cross-referencing with the phenotype in some way or other. In this recipe, we'll look at the sophisticated mixed linear model described by Yu et al in 2006 (Nature Genetics, 38:203-208). Describing the workings of the unified mixed linear model is beyond the scope of the recipe, but it is a suitable model for...

It is often of interest to know how often a sequence occurs in a sample of interest—that is, to estimate whether, in your particular sample, a locus has been duplicated or its copy number has increased. The locus could be anything from a gene at Kbp scale or a large section of DNA at Mbp scale. Our approach in this recipe will be to use HTS read coverage after alignment to estimate a background level of coverage and then inspect the coverage of our region of interest. The ratio of the coverage in our region of interest to the background level will give us an estimate of the copy number in the region. The recipe here is the first step. The background model we use is very simple—we calculate only a global mean, but we'll discuss some alternatives later. Also, this recipe does not cover ploidy—the number...