Genotype diagnostics for diversity outbred mice

Array intensities

Another useful indicator is the distribution of intensities on the MegaMUGA arrays. For speed of access, we store the array intensities in an .fst file, with the R package fst, as a matrix of SNPs x samples, with the two alleles at a SNP being stored in adjacent rows. See the R script geneseek2fst.R for the details of how we converted the FinalReport.txt files from GeneSeek into this format.

We’ll read in the full set of intensities and take the sum of the two allele intensities at each SNP.

int <- fst::read.fst("diag_cache/intensities.fst")
int <- int[seq(1, nrow(int), by=2),-(1:2)] + int[-seq(1, nrow(int), by=2),-(1:2)]
int <- int[,ind_ids(svenson)]

In the following plot, we display the distributions of array intensities (after a \(\log_{10}(x+1)\) transformation).

In the top panel, the arrays are sorted by the proportion of missing genotype data for the sample, and the curves connect various quantiles of the intensities. Hover over the top panel, and the corresponding histogram is shown below.

n <- names(sort(percent_missing, decreasing=TRUE))
iboxplot(log10(t(int[,n])+1), orderByMedian=FALSE, chartOpts=list(ylab="log10(SNP intensity + 1)"))

The first nine arrays are the ones that are missing > 20% of the genotypes; they all (well, except for F326) show a reduced median and long right tail. The next three arrays are not so extreme, but still show greater variability in array intensity than the bulk of the arrays, and F326 looks a bit like these. There are then a set of samples showing a spike of low intensity values.

For this particular set of arrays, a plot of the 1 %ile vs the 99 %ile is quite revealing. In the following, the orange points are those with > 20% missing genotypes, the pink points are the samples with 5-20% missing genotypes, and the blue points are the samples with 2-5% missing genotypes.

qu <- apply(int, 2, quantile, c(0.01, 0.99), na.rm=TRUE)
group <- (percent_missing >= 19.97) + (percent_missing > 5) + (percent_missing > 2) + 1
labels <- paste0(colnames(qu), " (", round(percent_missing), "%)")
iplot(qu[1,], qu[2,], indID=labels, group=group,
      chartOpts=list(xlab="1 %ile of array intensities",
                     ylab="99 %ile of array intensities",
                     pointcolor=c("#ccc", "slateblue", "Orchid", "#ff851b")))

The samples with >20% missing data (in orange) have somewhat reduced 1st percentile and elevated 90th percentile, except for one orange point (F326) that is close to the bulk of the samples. The samples with 5-9% missing data (in pink) are similar but less extreme. Many of the samples with 2-5% missing data show decreased 1st percentile, but many fit in with the bulk of the samples.

We’re again tempted to omit the 12 samples with >5% missing data, but we’ll continue to leave them in, for now.

Genotype frequencies

It can be useful to look at the genotype frequencies for the raw SNP data. We’ll focus on the autosomal markers.

We first pull out the SNP genotypes. They’ve been coded according to the allele frequencies in the founder strains, with 1 and 3 corresponding to the homozygotes for the major and minor alleles, respectively; let’s call them AA and BB (that is, A is the major allele). We also pull out the founder genotypes and count the number of founders with the minor allele. In studying the genotype frequencies, we’ll omit markers where any of the founders’ genotypes are missing (coded as 0).

g <- do.call("cbind", svenson$geno[1:19])
fg <- do.call("cbind", svenson$founder_geno[1:19])
g <- g[,colSums(fg==0)==0]
fg <- fg[,colSums(fg==0)==0]
fgn <- colSums(fg==3)

We now calculate the genotype frequencies, by individual and by marker. For the genotype frequencies by individual, we’ll calculate separate frequencies for the four groups of markers, split according to their minor allele frequency in the founder strains.

gf_ind <- vector("list", 4)
for(i in 1:4) {
    gf_ind[[i]] <- t(apply(g[,fgn==i], 1, function(a) table(factor(a, 1:3))/sum(a != 0)))
}

The following triangle plots show the genotype frequency distributions for the mice, among the four groups of markers with common minor allele frequency (MAF) in the founder strains. These plots make use of the fact that for a point within an equilateral triangle, the sum of the distances to the three sides is a constant.

par(mfrow=c(2,2), mar=c(0.6, 0.6, 2.6, 0.6))
for(i in 1:4) {
    triplot(c("AA", "AB", "BB"), main=paste0("MAF = ", i, "/8"))
    tripoints(gf_ind[[i]], pch=21, bg="lightblue")
    tripoints(c((1-i/8)^2, 2*i/8*(1-i/8), (i/8)^2), pch=21, bg="violetred")

    if(i>=3) { # label mouse with lowest het
      wh <- which(gf_ind[[i]][,2] == min(gf_ind[[i]][,2]))
      tritext(gf_ind[[i]][wh,,drop=FALSE] + c(0.02, -0.02, 0),
              names(wh), adj=c(0, 1))
    }

    # label other mice
    if(i==1) {
        lab <- rownames(gf_ind[[i]])[gf_ind[[i]][,2]>0.3]
    }
    else if(i==2) {
        lab <- rownames(gf_ind[[i]])[gf_ind[[i]][,2]>0.48]
    }
    else if(i==3) {
        lab <- rownames(gf_ind[[i]])[gf_ind[[i]][,2]>0.51]
    }
    else if(i==4) {
        lab <- rownames(gf_ind[[i]])[gf_ind[[i]][,2]>0.6]
    }

    for(ind in lab) {
        if(grepl("^F", ind) && i != 3) {
            tritext(gf_ind[[i]][ind,,drop=FALSE] + c(-0.01, 0, +0.01), ind, adj=c(1,0.5))
        } else {
            tritext(gf_ind[[i]][ind,,drop=FALSE] + c(0.01, 0, -0.01), ind, adj=c(0,0.5))
        }
    }
}

The majority of individuals are tightly clustered around the expected distribution (in pink). There are a few outliers with elevated heterozygosity, but these are all among the samples with highest rates of missing genotypes. Mouse F313 has slightly reduced heterozygosity in markers with MAF=3/8 or 4/8, but not by a worrisome amount.

Crossover counts

Counts of inferred crossovers can be a useful diagnostic for identifying problem samples, which may show an excessive number of apparent crossovers.

We first reconstruct the haplotypes of the DO mice, at each position calculating the probability of each of the 36 possible genotypes, given the observed marker data.

pr <- calc_genoprob(svenson, error_prob=0.002, map_function="c-f", cores=0)

There are a couple of different ways to get estimates of the numbers of crossovers. The simplest is to identify the most probable genotype at each position, and then count exchanges in such data. Here we use a cut-off of 0.5, leaving as missing any genotypes where the largest probability is <0.5.

m <- maxmarg(pr, minprob=0.5, cores=0)
nxo <- count_xo(m, cores=0)

The function count_xo() returns counts of crossovers on each chromosome (as columns) in each mouse. To get estimates of the total numbers of crossovers, genome-wide, we take the sums of each row:

totxo <- rowSums(nxo)

Here’s a plot of the number of crossovers vs the mouse index, colored by generation. (We have omitted the 9 samples with > 20% missing genotypes, as these samples show 501–2292 crossovers, well above the others.)

There is a lot of variation in the number of crossovers, and some overlap between generation 8 (in red, with an average of 304 crossovers) and generation 11 (in blue, with an average of 357 crossovers).

iplot(seq_along(totxo)[percent_missing < 19.97],
      totxo[percent_missing < 19.97],
      group=(svenson$covar$ngen=="11")[percent_missing < 19.97]+1,
      chartOpts=list(xlab="Mouse", ylab="Number of crossovers",
                     margin=list(left=80,top=40,right=40,bottom=40,inner=5),
                     axispos=list(xtitle=25,ytitle=50,xlabel=5,ylabel=5)))

We’re looking for mice with an excessive number of crossovers. Mice with very few crossovers would also deserve follow up. And note that we need to take the DO generation of an individual into account. But none of the mice seem unusual, in terms of numbers of crossovers, except for the 9 mice with > 20% missing genotypes which we omitted from this figure.

Here are the crossover counts for those 9 mice:

tmp <- cbind(percent_missing=round(percent_missing), total_xo=totxo)[percent_missing >= 19.97,]
tmp[order(tmp[,1]),]

##      percent_missing total_xo
## M387              20      501
## M411              21      501
## M412              26      562
## F326              31     1041
## M394              38      948
## M419              42      901
## M392              68     2292
## M405              73     2049
## M388              74     2187

Genotyping error LOD scores

Lincoln and Lander (1992) proposed the calculation of genotyping error LOD scores, to identify potential genotyping errors. For each SNP in each mouse, we calculate a LOD score statistic that measures the evidence that the particular genotype is in error, by comparing the observed genotype to the multipoint genotype probabilities, making use of the founder strains’ genotypes at the SNP. These are calculated using the function calc_errorlod().

e <- calc_errorlod(svenson, pr, cores=0)

The output is a list of matrices (the components being chromosomes). We can combine these into a single large matrix (mice × SNPs), and then get estimated genotyping error rates for each mouse, taking error LOD > 2 as a cutoff.

e <- do.call("cbind", e)
errors_ind <- rowSums(e>2)/n_typed(svenson)*100
lab <- paste0(names(errors_ind), " (", myround(percent_missing,1), "%)")
iplot(seq_along(errors_ind), errors_ind, indID=lab,
      chartOpts=list(xlab="Mouse", ylab="Percent genotyping errors", ylim=c(0, 4.1),
                     axispos=list(xtitle=25, ytitle=50, xlabel=5, ylabel=5)))

There are 9 mice with estimated genotyping error rates > 1%, and these are exactly the mice with > 20% missing genotypes. There are 3 mice with error rates in the 0.5-1% range, and these are the mice with around 9% missing genotypes.

The mean error rate of mice with <5% missing genotype data is 9 in 10,000.

Apparent genotyping errors

A more direct approach to look for genotyping errors is to get inferred SNP genotypes and compare them to the observed ones. We use the result of maxmarg(), calculated above, and then predict_snpgeno(). The output is again a list of matrices, one per chromosome, but we’ll combine them into a single matrix.

snpg <- predict_snpgeno(svenson, m, cores=0)
snpg <- do.call("cbind", snpg)

We’ll similarly grab the raw SNP genotypes.

gobs <- do.call("cbind", svenson$geno)
gobs[gobs==0] <- NA

The apparent error rates derived by comparing the observed and predicted SNP genotypes are basically the same as using a cutoff of 0 on the error LOD score. The rates differ most for the samples missing a lot of genotypes.

par(pty="s")
err_direct <- rowMeans(snpg != gobs, na.rm=TRUE)*100
errors_ind_0 <- rowSums(e > 0)/n_typed(svenson)*100
par(mar=c(4.1,4.1,0.6, 0.6))
grayplot(errors_ind_0, err_direct,
         xlab="Percent errors (error LOD > 0)",
         ylab="Percent errors (obs vs predicted)",
         xlim=c(0, 21.2), ylim=c(0, 21.2))
abline(0,1,lty=2, col="gray60")

Here are just the samples with estimated error rates < 0.4% (that is, omitting 13 samples with error rates > 0.5%).

par(pty="s")
par(mar=c(4.1,4.1,0.6, 0.6))
grayplot(errors_ind_0, err_direct,
         xlab="Percent errors (error LOD > 0)",
         ylab="Percent errors (obs vs predicted)",
         xlim=c(0, 0.4), ylim=c(0, 0.4))
abline(0,1,lty=2, col="gray60")

Using a cutoff of 0 on the genotyping error LOD score, the average error rate in mice with <5% missing genotypes is 12 per 10,000.

Missing data

It can also be useful to look at the proportion of missing genotypes by marker. Markers with a lot of missing data were likely difficult to call, and so the genotypes that were called may contain a lot of errors.

pmis_mar <- n_missing(svenson, "marker", "proportion")*100

Here’s a histogram of the proportion of missing genotypes by marker. The tick marks below the histogram indicate individual values.

par(mar=c(5.1,0.6,0.6, 0.6))
hist(pmis_mar, breaks=seq(0, 100, length=201),
     main="", yaxt="n", ylab="", xlab="Percent missing genotypes")
rug(pmis_mar)

The vast majority of markers have very little missing data: there are 41,931 markers (60%) that are missing no data, and 64,003 markers (92%) that are missing <2%.

But some markers are missing a lot of data, including 45 that are missing more than 50% of genotypes. 240 markers are missing more than 30% of genotypes, while 742 markers are missing more than 15% of genotypes. Recall that there are a total of 69,339, markers, so 742 markers is just 1.1% of the total.

For the markers with lots of missing genotypes, it’s not necessarily the case that the remaining genotypes are full of errors, but in studying the allele intensities at these SNPs, it does appear that for the bulk of such markers, the genotypes are not being called appropriately.

Genotype frequencies

The genotype frequencies by marker can be informative, just as genotype frequencies by individual were.

gf_mar <- t(apply(g, 2, function(a) table(factor(a, 1:3))/sum(a != 0)))
gn_mar <- t(apply(g, 2, function(a) table(factor(a, 1:3))))

Here are scatter plots of the genotype frequencies by marker, split according to the minor allele frequency in the 8 founders. A pink dot is placed at the expected frequencies.

par(mfrow=c(2,2), mar=c(0.6, 0.6, 2.6, 0.6))
for(i in 1:4) {
    triplot(c("AA", "AB", "BB"), main=paste0("MAF = ", i, "/8"))
    z <- gf_mar[fgn==i,]
    z <- z[rowSums(is.na(z)) < 3,]
    tripoints(z, pch=21, bg="gray80", cex=0.6)
    tripoints(c((1-i/8)^2, 2*i/8*(1-i/8), (i/8)^2), pch=21, bg="violetred")
}

The bulk of the markers seem well behaved, but there are a number of markers with unusual genotype frequencies. There are 6 markers that show no homozygotes for the major allele. (These sit on the left edge.) There are 9 markers that are monomorphic (100% AA genotypes; lower-right vertex). And there are 35 markers with some of each homozygote but no heterozygotes (bottom edge).

Further, there are a bunch of markers where the minor allele appears to be private to one founder strain (upper-left panel) but show a high frequency of minor alleles in the DO offspring. For example, there are 11 markers where the minor allele is private to one founder line but has >50% allele frequency in the DO offspring, and 56 such private SNPs have >30% minor allele frequency in the DO offspring. But it turns out that all but two of these are private to WSB and reside on chromosome 2. This is due to meiotic drive of the WSB allele at the R2D2 locus; see Didion et al. (2016).

The other two are also interesting: JAX00463240 on chr 18 (private to CAST) and UNC9169851 on chr 5 (private to NZO) both have >30% minor allele frequency in the DO offspring. Both seem to be in regions where one of the founder alleles is unusually frequent.

Genotyping errors

Let’s look at the apparent genotyping errors by marker (again starting with a cutoff of error LOD score > 2), plotted against the amount of missing data.

errors_mar <- colSums(e>2)/n_typed(svenson, "marker")*100
grayplot(pmis_mar, errors_mar,
         xlab="Proportion missing", ylab="Proportion genotyping errors")

Most of the markers show no evidence for genotyping errors. Only 2,716 (3.9% of markers) show at least one apparent error, but 325 show >5% errors. Markers with higher rates of missing genotypes tend to show higher errors rates, though this is difficult to detect in the above figure.

Array intensities

Scatterplots of SNP allele intensities are quite revealing about the cause, particularly when we look at the scatterplots colored by both the observed SNP genotypes and by the predicted SNP genotypes, given the multipoint SNP information.

Let us first use the inferred 36-state genotypes to get predicted SNP genotypes, which we then combine into a single wide matrix.

snpg <- predict_snpgeno(svenson, m, cores=0)
snpg <- do.call("cbind", snpg)

Here is a function to plot the SNP intensities, plus several helper functions.

# get data for one marker
read_intensities <-
    function(file="../RawData/intensities.fst", marker=NULL, snps=NULL)
{
    if(is.null(snps)) {
        snps <- fst::read.fst(file, column="snp")[,1]
    }

    if(is.null(marker)) { # if marker is NULL, return the snp names
        return(snps)
    } else {
        marker <- as.character(marker)
    }

    n <- match(marker, snps)
    if(is.na(n)) stop("marker ", marker, " not found")

    reorg_intensities( fst::read.fst(file, from=n, to=n+1) )
}


# reorganize SNP intensity data (for one marker)
reorg_intensities <-
    function(intensities, marker=NULL)
{
    if(!is.null(marker)) {
        marker <- as.character(marker)
        wh <- which(intensities[,1] == marker)
        if(sum(wh==0)) stop("Marker ", marker, " not found")
        intensities <- intensities[wh,]
    }

    result <- data.frame(X=as.numeric(intensities[1,-(1:2)]),
                         Y=as.numeric(intensities[2,-(1:2)]))
    rownames(result) <- colnames(intensities)[-(1:2)]

    result
}

# combine intensities and genotypes
grab_gni <-
    function(marker=NULL, cross=svenson, intensities_file="../RawData/intensities.fst",
             intensities_data=NULL, drop_bad_samples=TRUE)
{
    if(is.null(intensities_data)) {
        # get intensities
        int <- read_intensities(intensities_file, marker=marker)
    } else {
        int <- reorg_intensities(intensities_data, marker)
    }
    marker <- as.character(marker)

    # get genotypes
    chr <- qtl2::find_markerpos(cross, marker)$chr
    g <- cross$geno[[chr]][,marker]
    gg <- setNames(rep(0, nrow(int)), rownames(int))
    gg[names(g)] <- g

    if(drop_bad_samples) {
        int <- int[names(g),,drop=FALSE]
        gg <- gg[names(g)]
    }

    cbind(int, g=gg)
}

# load intensities and plot them, colored by genotype calls
#
# drop_bad_samples: if TRUE, don't plot points for samples that are not in the cross object
#                   (e.g., were omitted previously as being bad DNAs)
#
plot_intensities <-
    function(marker, cross=svenson, intensities_file="diag_cache/intensities.fst",
             intensities_data=NULL, drop_bad_samples=TRUE, geno=NULL, ...)
{
    if(is.character(marker) || is.factor(marker)) {
        marker <- as.character(marker)
        gni <- grab_gni(marker, cross=cross, intensities_file=intensities_file,
                        intensities_data=intensities_data, drop_bad_samples=drop_bad_samples)
    } else {
        gni <- marker
    }

    if(!is.null(geno)) {
        if(is.logical(geno)) geno <- geno + 1 # FALSE/TRUE -> 1/2
        gni[names(geno),"g"] <- geno
    }

    internal_plot <-
        function(pch=21, bg=broman::brocolors("f2"),
                 xlab="allele 1", ylab="allele 2",
                 xlim=c(0, max(gni$X, na.rm=TRUE)),
                 ylim=c(0, max(gni$Y, na.rm=TRUE)),
                 ...)
        {
            grayplot(gni$X, gni$Y, pch=pch, bg=bg[match(gni$g, c(1:3,0))],
                     xlab=xlab, ylab=ylab,
                     xlim=xlim, ylim=ylim, ...)
        }

    internal_plot(...)

    invisible(gni)
}

The scatterplots below show examples of the allele intensities for four SNP markers. Each point is a mouse sample, and the x and y axes are for the array intensities of the two alleles at the SNP. In the left panels, the points are colored by the observed SNP genotypes, with yellow and blue corresponding to the two homozygotes and green to the heterozygotes; gray points indicate missing genotypes. In the right panels, the points are colored by the predicted SNP genotypes; gray points indicate samples for which the SNP genotype could not be predicted.

snpg[is.na(snpg)] <- 0

markers <- c("JAX00279019",
             "UNC12329705",
             "backupUNC140137407",
             "UNC20478577")

par(mfrow=c(4,2), mar=c(3.1,3.1,2.1,1.1))
for(i in seq_along(markers)) {
    mar <- markers[i]

    plot_intensities(mar, mgp.x=c(1.4,0.3,0), mgp.y=c(1.8,0.3,0))
    title(main=mar, line=1.05, cex.main=0.9)
    title(main="(observed genotypes)", cex.main=0.8, line=0.25)

    u <- par("usr")
    mtext(side=3, adj=0, LETTERS[i*2-1])

    plot_intensities(mar, geno=snpg[,mar],
                     mgp.x=c(1.4,0.3,0), mgp.y=c(1.8,0.3,0))
    title(main=mar, line=1.05, cex.main=0.9)
    title(main="(predicted genotypes)", cex.main=0.8, line=0.25)

    mtext(side=3, adj=0, LETTERS[i*2])
}

The top panels (A and B) are for a well-behaved SNP, JAX00279019. There are three clear genotype clusters, and the observed and predicted genotypes match.

In panels C and D (UNC12329705), the genotype calling algorithm went awry. There are three distinct clusters of samples, but the homozygotes for the major allele appear to have been called heterozygotes.

In panels E and F (backupUNC140137407), there appears to be an extra cluster of samples with distinct genotype, but the genotype calling algorithm classified them as homozygotes and they instead (based on surrounding markers) appear to be heterozygotes.

Panels G and H (UNC20478577) are for a case where the genotype clusters are overlapping and the genotyping calling algorithm had difficulty distinguishing them.