library(comapr)
library(GenomicRanges)
#> Loading required package: stats4
#> Loading required package: BiocGenerics
#>
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:stats':
#>
#> IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#>
#> Filter, Find, Map, Position, Reduce, anyDuplicated, aperm, append,
#> as.data.frame, basename, cbind, colnames, dirname, do.call,
#> duplicated, eval, evalq, get, grep, grepl, intersect, is.unsorted,
#> lapply, mapply, match, mget, order, paste, pmax, pmax.int, pmin,
#> pmin.int, rank, rbind, rownames, sapply, setdiff, sort, table,
#> tapply, union, unique, unsplit, which.max, which.min
#> Loading required package: S4Vectors
#>
#> Attaching package: 'S4Vectors'
#> The following object is masked from 'package:utils':
#>
#> findMatches
#> The following objects are masked from 'package:base':
#>
#> I, expand.grid, unname
#> Loading required package: IRanges
#> Loading required package: GenomeInfoDb
library(BiocParallel)
Install via BiocManager as follow:
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("comapr")
comapr
lets you interrogate the genotyping results for a sequence of markers
across the chromosome and detect meiotic crossovers from genotype shifts. In
this document, we demonstrate how genetic distances are calculated from genotyping
results for a group of samples using functions available from comapr
.
The package includes a small data object that contains 70 markers and their genotyping results for 22 samples. The 22 samples are the BC1F1 progenies generated by
Therefore, the genotype shifts detected in BC1F1 samples represent crossovers that have happened during meiosis for the F1 parents.
BiocParallel::register(SerialParam())
BiocParallel::bpparam()
#> class: SerialParam
#> bpisup: FALSE; bpnworkers: 1; bptasks: 0; bpjobname: BPJOB
#> bplog: FALSE; bpthreshold: INFO; bpstopOnError: TRUE
#> bpRNGseed: ; bptimeout: NA; bpprogressbar: FALSE
#> bpexportglobals: FALSE; bpexportvariables: FALSE; bpforceGC: FALSE
#> bpfallback: FALSE
#> bplogdir: NA
#> bpresultdir: NA
data(snp_geno_gr)
data(parents_geno)
In order to detect crossovers from the 22 samples’ genotype results (BC1F1
samples), we need to format the result into Homo_ref
, Homo_alt
and Het
encodings. This can be simply done through comparing with the parents’ genotypes.
Here, the genotype noted as “Fail” will be converted to NA. Please also note that we supplied “Homo_ref” as one of the fail options because “Homo_ref” is not in the possible genotypes of markers in BC1F1 samples.
corrected_geno <- correctGT(gt_matrix = mcols(snp_geno_gr),
ref = parents_geno$ref,
alt = parents_geno$alt,
fail = "Fail",
wrong_label = "Homo_ref")
mcols(snp_geno_gr) <- corrected_geno
The corrected_geno
matrix
head(mcols(snp_geno_gr)[,1:5])
#> DataFrame with 6 rows and 5 columns
#> X92 X93 X94 X95 X96
#> <character> <character> <character> <character> <character>
#> 1 Homo_alt Homo_alt Homo_alt Homo_alt Homo_alt
#> 2 Homo_alt Homo_alt Homo_alt Homo_alt Homo_alt
#> 3 Homo_alt Homo_alt Homo_alt Homo_alt Homo_alt
#> 4 Homo_alt Homo_alt Homo_alt Homo_alt Homo_alt
#> 5 NA Homo_alt Homo_alt Homo_alt Homo_alt
#> 6 Het Homo_alt Homo_alt Homo_alt Homo_alt
Note that there are missing values in this resulting matrix that can be resulted from:
NA
In this step, we try to identify markers that have NA
genotype across many
samples or samples that have a lot markers failed for removal. We use the
countGT
function for find bad markers/samples.
genotype_counts <- countGT(mcols(snp_geno_gr))
genotype_counts$plot
The number of markes and samples are saved in a list returned by countGT
genotype_counts$n_markers
#> [1] 30 33 33 33 33 33 31 31 32 33 33 33 32 33 33 33 32 32 33 32 33 33
genotype_counts$n_samples
#> [1] 22 22 22 22 18 22 22 22 22 22 22 21 22 22 22 22 22 22 22 22 22 22 22 22 22
#> [26] 22 22 22 16 22 22 21 22
We now filter out markers/samples by using function filterGT
. min_markers
specifies at least how many markers a sample needs to be kept. Likewise for
min_samples
.
A printed message contains information about how much markers or samples have been filtered.
corrected_geno <- filterGT(snp_geno_gr,
min_markers = 30,
min_samples = 2)
#> filter out 0 marker(s)
#> filter out 0 sample(s)
Sample duplicates are identified by finding samples that share exactly same genotypes
across all available markers. findDupSamples
can be applied and a threshold
value is provided and used as a cut-off on the percentage of same genotype markers
the duplicated samples should share.
dups <- findDupSamples(mcols(corrected_geno),
threshold = 0.99)
dups
#> X98
#> [1,] "X98"
#> [2,] "X99"
Now we remove the duplicated samples.
mcols(corrected_geno) <- mcols(corrected_geno)[,colnames(mcols(corrected_geno))!="X98"]
#corrected_geno
Crossovers are detected and counted through examining the patterns of genotypes
along the chromosome. When there is a shift from one genotype block to another,
a crossover is observed. This is done through calling countCOs
function which
returns a GRange
object with crossover counts for the list of marker intervals.
The crossover count values in the columns can be non-integer when one observed crossover can not be determined to be completely distributed to the marker interval in the corresponding row. The observed crossover is then distributed to the adjacent intervals proportionally to their interval base pair sizes.
marker_gr_cos <- countCOs(corrected_geno)
marker_gr_cos[1:5,1:5]
#> GRanges object with 5 ranges and 5 metadata columns:
#> seqnames ranges strand | X92 X93 X94
#> <Rle> <IRanges> <Rle> | <numeric> <numeric> <numeric>
#> [1] 1 6655965-21638463 * | 0.9358742 0 0
#> [2] 1 21638465-22665059 * | 0.0641258 0 0
#> [3] 1 22665061-34590735 * | 0.0000000 0 0
#> [4] 1 35033642-38996025 * | 0.0000000 0 0
#> [5] 2 4248665-5348752 * | 0.0000000 0 0
#> X95 X96
#> <numeric> <numeric>
#> [1] 0 0
#> [2] 0 0
#> [3] 0 0
#> [4] 0 0
#> [5] 0 0
#> -------
#> seqinfo: 4 sequences from an unspecified genome; no seqlengths
The genetic distances of marker intervals are calcuated based on the crossover
rates via applying mapping function, Kosambi
or Haldane
by calling the
calGeneticDist
function. The returned genetic distances are in unit of
centiMorgan.
dist_gr <- calGeneticDist(marker_gr_cos,
mapping_fun = "k")
dist_gr[1:5,]
#> GRanges object with 5 ranges and 1 metadata column:
#> seqnames ranges strand | kosambi_cm
#> <Rle> <IRanges> <Rle> | <numeric>
#> [1] 1 6655965-21638463 * | 14.36279
#> [2] 1 21638465-22665059 * | 5.08472
#> [3] 1 22665061-34590735 * | 9.64156
#> [4] 1 35033642-38996025 * | 9.64156
#> [5] 2 4248665-5348752 * | 4.77638
#> -------
#> seqinfo: 4 sequences from an unspecified genome; no seqlengths
Alternatively, instead of returning the genetic distances in supplied marker
intervals, we can specify a bin_size
which tells calGeneticDist
to return
calculated genetic distances for equally sized chromosome bins.
dist_bin_gr <- calGeneticDist(marker_gr_cos,bin_size = 1e6)
dist_bin_gr[1:5,]
#> GRanges object with 5 ranges and 1 metadata column:
#> seqnames ranges strand | kosambi_cm
#> <Rle> <IRanges> <Rle> | <numeric>
#> [1] 1 1-998753 * | 0
#> [2] 1 998754-1997506 * | 0
#> [3] 1 1997507-2996258 * | 0
#> [4] 1 2996259-3995011 * | 0
#> [5] 1 3995012-4993763 * | 0
#> -------
#> seqinfo: 4 sequences from mm10 genome
With genetic distances calculated, we can do a sum of all genetic distances across all marker intervals. We can see that we got the same total genetic distances for marker based intervals and the equally binned intervals.
sum(dist_bin_gr$kosambi_cm)
#> [1] 168.0926
sum(dist_gr$kosambi_cm)
#> [1] 168.0926
comapr
also includes functions for visulising genetic distances of marker
intervals or binned intervals.
plotGeneticDist(dist_bin_gr,chr = "1")
We can also plot the cumulative genetic distances of certain chromosomes
plotGeneticDist(dist_bin_gr,chr=c("1"),cumulative = TRUE)
Multiple chromosomes
plotGeneticDist(dist_bin_gr,cumulative = TRUE)
comapr
implements a whole genome plot function too, that takes all chromosomes
available in the result and plot a cumulative genetic distances by cumulatively
summing all intervals across all chromosomes.
plotWholeGenome(dist_bin_gr)