组学内样本关系验证

方法检测数据链接错误omicsPrint使用人工生成的数据和实验数据进行打包。另外还有几个小插图，显示了在不同环境下对进一步实验数据的使用，即450k DNA甲基化和估算基因型。

创建玩具数据

这里我们从集合中生成一个包含100个随机抽取的整数的单个向量;1,2,3，代表来自单个个体的100个SNP呼叫。通过随机交换一定比例的snp产生了另外三个个体。交换5个snp会引入一些不匹配，模拟同一个体被测量两次(重复)但测量误差的情况。交换50%的SNPs将类似于父母和后代之间的基因型差异。交换所有snp将导致类似比较两个不相关的个体的情况。

Swap <- function(x, frac=0.05) {n <- length(x) k <- floor(n*frac) x1 <- sample(1:n,k) x2 <- sample(1:n,k) ##可能是重叠x[x2] <- x[x1] x} x1 <- 1 + rbinom(100, size=2, prob=1/3) x2 <- Swap (x1, 0.05) ##相关例如，复制x3 <- Swap (x1, 0.5) ##相关例如，parent off spring x4 <- Swap (x1, 1) ##不相关x <- cbind(x1, x2, x3, x4)

现在x包含了四个个体使用的100个snp头我们可以检查前六个snp。

头(x)

## x1 x2 x3 x4 ## [1，] 1 1 2 2 ## [2，] 2 2 2 1 ## [3，] 1 1 3 1 ## [4，] 1 1 2 1 ## [5，] 1 1 2 2 ## [6，] 2 2 2 2 2

运行`allelesharing`算法

我们使用状态标识(IBS)的snp集来推断样本关系。看到（?？?），描述了这种方法应用于遗传数据。简单地说，在每个样本对之间，计算ibs向量，这是个体之间遗传距离的测量。接下来，通过其均值和方差对向量进行汇总。均值为2，方差为0表示样本相同。

data <- allelessharing (x, verbose=TRUE)

##哈希关系

##修剪100个snp…

## 0个snp被删除，因为呼叫速率低!

## 0个样本被删除，因为调用的snp太少!

##使用100个多态snp确定等位基因共享。

##运行“平方”IBS算法!

snp集可能包含无信息的snp，质量差的snp，甚至可能缺少snp。实现以下修剪步骤，以产生信息量最大的snp集(阈值可以调整，见alleleSharing ?)：

如果> %的样本中缺失一个SNP，则去除(callRate = 0.95）
一个样本中至少有2/3的SNPs被称为(coverageRate = 2/3）
如果一个SNP违反了Hardy-Weinberg平衡(= 0）.
如果次要等位基因频率低于给定的阈值(Maf = 0）

Hardy-Weinberg检验统计量使用\ (\ x ^ 2 \)-test和Bonferonni多重测试校正。

数据

## var colnames。x colnames。yrelation ## 1 2.00 0.0000000 x1 x1 identical ## 2 1.96 0.0589899 x2 x1 unrelated ## 3 1.62 0.3187879 x3 x1 unrelated ## 4 1.40 0.3636364 x4 x1 unrelated ## 5 2.00 0.0000000 x2 x2 identical ## 6 1.60 0.3434343 x3 x2 unrelated ## 7 1.38 0.3793939 x4 x2 unrelated ## 8 2.00 0.0000000 x3 x3 identical ## 9 1.40 0.4242424 x4 x3 unrelated ## 10 2.00 0.0000000 x4 x4 identical

默认情况下，除了self-self关系外，不假定任何关系。

输出为data.frame包含所有与IBS在snp集上的均值和方差的成对比较，以及报告的样本关系，包括标识符。

报告不匹配并提供图形汇总

由于我们提供了一个已知关系的列表，并假设大多数是正确的，我们可以构建一个分类器来发现错误分类的关系。使用线性判别分析生成混淆矩阵，然后使用该矩阵图形化地表示分类边界，并生成带有错误分类样本对的输出文件。

mismatches <- inferRelations(data)

IBS均值和方差与分类边界的散点图，用于样本之间的成对比较，而不使用人工生成的数据指定样本关系。

不匹配

## var colnames。x colnames。yrelation predicted ## 2 1.96 0.0589899 x2 x1 unrelated identical

有一个错误分类的样本，即我们引入的复制，但没有先验地指定为存在的关系。样本之间的真实关系x1和样本x2是同一关系。此外，我们看到两个样本对的平均IBS为1.96，方差为0.06，这表明这些对也具有相当数量的等位基因。如果已知，可以在分析之前指定这种关系。

关系<- expand。Grid (idx = colnames(x)， idy= colnames(x)) relations$relation_type <- "不相关"关系$relation_type[relations$idx == relations$idy] <- "相同"关系$relation_type[c(2,5)] <- "相同" ##复制关系$relation_type[c(3,7,9,10)] <- "父-子"关系

# # idx补贴relation_type x1相同的# # # # 1 x1 2 x2 x1相同的# # 3 x3 x1父的后代# # 4 x4 x1无关的# # 5 (x1, x2)相同的# # 6 x2相同# # 7 x3 x2父的后代# # 8 x4 x2无关的# # 9 x1 x3父的后代# # 10 x2 x3父的后代# # 11 x3 x3相同的# # 12 x4 x3无关的x1 x4无关的# # # 13 # 14 # #无关的x2 x4 15 x3 x4无关的# # 16 x4 x4相同

现在重新运行已知关系提供的等位基因共享算法。

data <- allelesshare (x, relations=relations)

##哈希关系

##修剪100个snp…

## 0个snp被删除，因为呼叫速率低!

## 0个样本被删除，因为调用的snp太少!

##使用100个多态snp确定等位基因共享。

##运行“平方”IBS算法!

数据

## var colnames。x colnames。yrelation ## 1 2.00 0.0000000 x1 x1 identical ## 2 1.96 0.0589899 x2 x1 identical ## 3 1.62 0.3187879 x3 x1 parent offspring ## 4 1.40 0.3636364 x4 x1 unrelated ## 5 2.00 0.0000000 x2 x2 identical ## 6 1.60 0.3434343 x3 x2 parent offspring ## 7 1.38 0.3793939 x4 x2 unrelated ## 8 2.00 0.0000000 x3 x3 identical ## 9 1.40 0.4242424 x4 x3 unrelated ## 10 2.00 0.0000000 x4 x4 identical

mismatches <- inferRelations(data)

IBS均值和方差与分类边界的散点图，用于样本之间的成对比较，并使用人工生成的数据指定样本关系。

不匹配

## [1] mean var colnames。x colnames。yrelation predicted ## <0 rows> (or 0-length row.names)

所有错误分类的关系都得到了解决。

跨组学数据类型样本关系验证

前面的例子展示了如何在单个组学数据类型中执行样本关系验证。如果从不同的组学数据类型中获得了第二组snp(并且这些snp部分重叠)，omicsPrint可用于验证组学类型之间的关系，例如，确定两个组学数据类型是否确实为同一个体生成，以排除或检测样本混淆。

在这个人为的例子中，从不同的组型中选择80个snp的随机子集作为snp集。第一次运行没有提供已知的关系。

rownames(x) <- paste0("rs"， 1:100) y <- x[sample(1:100, 80)，] data <- allelessharing (x, y)

##哈希关系

##修剪100个snp…

## 0个snp被删除，因为呼叫速率低!

## 0个样本被删除，因为调用的snp太少!

##修剪80个snp…

## 0个snp被删除，因为呼叫速率低!

## 0个样本被删除，因为调用的snp太少!

##使用80个多态性snp确定等位基因共享。

##运行“矩形”IBS算法!

请注意，在这两种数据类型上执行剪接，并自动使用一组重叠的snp(80)，前提是行名x而且y是相同的(这也适用于示例关系，其中关系标识符idx而且补贴的栏名应该匹配x而且y）.

数据

## var colnames。x colnames。yrelation ## 1 2.000 0.00000000 x1 x1 identical ## 2 1.950 0.07341772 x2 x1 unrelated ## 3 1.625 0.31329114 x3 x1 unrelated ## 4 1.425 0.34873418 x4 x1 unrelated ## 5 1.950 0.07341772 x1 x2 unrelated ## 6 2.000 0.00000000 x2 x2 identical ## 7 1.600 0.34430380 x3 x2 unrelated ## 8 1.400 0.36962025 x4 x2 unrelated ## 9 1.625 0.31329114 x1 x3 unrelated ## 10 1.600 0.34430380 x2 x3 unrelated ## 11 2.000 0.00000000 x3 x3 identical ## 12 1.375 0.41455696 x4 x3 unrelated ## 13 1.425 0.34873418 x1 x4 unrelated ## 14 1.400 0.36962025 x2 x4 unrelated ## 15 1.375 0.41455696 x3 x4 unrelated ## 16 2.000 0.00000000 x4 x4 identical

mismatches <- inferRelations(data)

IBS均值和方差与分类边界的散点图，用于样本之间的成对比较，而不使用人工数据指定样本关系。

不匹配

## var colnames。x colnames。yrelation predicted ## 2 1.95 0.07341772 x2 x1 unrelated identical ## 5 1.95 0.07341772 x1 x2 unrelated identical

现在开始提供已知的关系。

data <- allelessharing (x, y, relations)

##哈希关系

##修剪100个snp…

## 0个snp被删除，因为呼叫速率低!

## 0个样本被删除，因为调用的snp太少!

##修剪80个snp…

## 0个snp被删除，因为呼叫速率低!

## 0个样本被删除，因为调用的snp太少!

##使用80个多态性snp确定等位基因共享。

##运行“矩形”IBS算法!

数据

## var colnames。x colnames。yrelation ## 1 2.000 0.00000000 x1 x1 identical ## 2 1.950 0.07341772 x2 x1 identical ## 3 1.625 0.31329114 x3 x1 parent offspring ## 4 1.425 0.34873418 x4 x1 unrelated ## 5 1.950 0.07341772 x1 x2 identical ## 6 2.000 0.00000000 x2 x2 identical ## 7 1.600 0.34430380 x3 x2 parent offspring ## 8 1.400 0.36962025 x4 x2 unrelated ## 9 1.625 0.31329114 x1 x3 parent offspring ## 10 1.600 0.34430380 x2 x3 parent offspring ## 11 2.000 0.00000000 x3 x3 identical ## 12 1.375 0.41455696 x4 x3 unrelated ## 13 1.425 0.34873418 x1 x4 unrelated ## 14 1.400 0.36962025 x2 x4 unrelated ## 15 1.375 0.41455696 x3 x4 unrelated ## 16 2.000 0.00000000 x4 x4 identical

mismatches <- inferRelations(data)

IBS均值和方差与分类边界的散点图，用于样本之间的成对比较，并使用人工数据指定样本关系。

不匹配

## [1] mean var colnames。x colnames。yrelation predicted ## <0 rows> (or 0-length row.names)

因此，提供已知的、真实的关系不会产生错误分类的样本关系。

一个使用真实世界甲基化数据的例子`SummarizedExperiment`

在这里，我们将展示如何在公开的DNA甲基化数据上验证样本关系。这里使用的数据集包含一对同卵双胞胎。我们将从GEO中提取beta值矩阵GSE100940，印刷中的纸张．

方法从GEO中提取数据GEOquery包。

library(GEOquery) library(summarizeexperiment) file <- tempfile(fileext = ".txt.gz") download.file("ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE100nnn/GSE100940/matrix/GSE100940_series_matrix.txt.gz"， file) gset <- getGEO(filename=file, gettgpl =FALSE)

接下来，我们将返回的对象转换为SummarizedExperiment：

se <- maksummarizedexperimentfromexpressionset (gset) se

##类:rangedsummarizeexperimental ## dim: 485577 24 ##元数据(3):experimentData注释protocolData ## assays(1): exprs ## rownames(485577): cg00000029 cg00000108…rs966367 rs9839873 ## rowData names(0): ## colnames(24): GSM2696882 GSM2696883…GSM2696904 GSM2696905 ## colData names(33): title geo_accession…data_row_count gender.ch1

样本数据可以从SummarizedExperiment-object使用colData函数，我们可以看到每个样本属于哪一对双胞胎source_name_ch1字段。利用这些知识，我们可以构建一个预期关系表:

r <- expand.grid(idx=colnames(se)， idy=colnames(se)) r$Xpair <- sapply(strsplit(as.character(colData(se)[r$idx， "source_name_ch1"])， split = "_")， head, 1) r$Ypair <- sapply(strsplit(as.character(colData(se)[r$idy， "source_name_ch1"])， split = "_")， head, 1) r$relation_type <- "不相关" r$relation_type[r$Xpair == r$idy] <- "孪生" r$relation_type[r$idx == r$idy] <- "相同" head(r)

## idx idy Xpair配对关系## 1 GSM2696882 GSM2696882 M01 M01相同## 2 GSM2696883 GSM2696882 M01 M01不相关## 3 GSM2696884 GSM2696882 M02 M01不相关## 5 GSM2696886 GSM2696882 M03 M01不相关## 6 GSM2696887 GSM2696882 M03 M01不相关

阵列上的几个探针包含在不同群体中频繁出现的snp（?？?；?？?）．我们可以使用这些来验证预期的关系。我们已经在这个包中提供了这些数据。

现在我们选择了东亚人群中可能受多态性SNPS影响的cpg，因为这些样本来自韩国:

data(hm450.manifest.pop. gonl) cpgs <- names(hm450.manifest.pop. gonl)GoNL[mcols(hm450.manifest.pop.GoNL)$MASK.snp5.EAS]) se <- se[cpgs，]

接下来，使用我们的增强型k均值算法将beta值转换为基因型:

dnamCalls <- beta2genotype(se, assayName = "exprs") dim(dnamCalls)

## [1] 821 24

dnamCalls [1:5, 1:5)

## GSM2696882 GSM2696883 GSM2696884 GSM2696885 GSM2696886 ## cg09762182 3 3 22 22 ## cg24345856 3 3 22 3 3 ## cg12213037 22 3 3 3 3 3 3

基于DNA甲基化的基因型调用可以直接提供给等位基因共享算法进行组内样本匹配:

data <- allelessharing (dnamCalls, relations = r, verbose = TRUE)

##哈希关系

##修剪821个snp…

## 0个snp被删除，因为呼叫速率低!

## 0个样本被删除，因为调用的snp太少!

##利用821个多态snp确定等位基因共享。

##运行“平方”IBS算法!

## 25 / 300(8.33%)…

mismatches <- inferRelations(data)

具有分类边界的IBS均值和方差散点图，用于由同卵双胞胎组成的样本之间的成对比较。

不匹配

## var colnames。x colnames。yrelation predicted ## 2 1.995128 0.004854282 GSM2696883 GSM2696882 twin identical ## 49 1.997564 0.002433083 GSM2696885 GSM2696884 twin identical ## 92 1.991474 0.008463801 GSM2696887 GSM2696886 twin identical ## 131 1.995128 0.004854282 GSM2696889 GSM2696888 twin identical ## 166 1.995128 0.004854282 GSM2696891 GSM2696890 twin identical ## 197 1.995128 0.004854282 GSM2696893 GSM2696892 twin identical ## 224 1.996346 0.003645168 GSM2696895 GSM2696894 twin identical ## 247 2.000000 0.000000000 GSM2696897 GSM2696896 twin identical ## 266 1.995128 0.004854282 GSM2696899 GSM2696898 twin identical ## 281 1.992692 0.007263599 GSM2696901 GSM2696900 twin identical ## 292 1.993910 0.006060426 GSM2696903 GSM2696902 twin identical ## 299 1.993910 0.006060426 GSM2696905 GSM2696904 twin identical

这对双胞胎被认为是一模一样的。这并不意外，因为它们是同卵的。

SessionInfo

sessionInfo ()

## R版本4.2.1(2022-06-23)##平台:x86_64-pc-linux-gnu(64位)##运行在Ubuntu 20.04.5 LTS ## ##矩阵产品:默认## BLAS: /home/biocbuild/bbs-3.16-bioc/R/lib/libRblas。/home/biocbuild/bbs-3.16-bioc/R/lib/libRlapack。所以## ## locale: ## [1] LC_CTYPE=en_US。UTF-8 LC_NUMERIC= c# # [3] LC_TIME=en_GB LC_COLLATE= c# # [5] LC_MONETARY=en_US。utf - 8 LC_MESSAGES = en_US。UTF-8 ## [7] LC_PAPER=en_US。UTF-8 LC_NAME= c# # [9] LC_ADDRESS=C lc_phone = c# # [11] LC_MEASUREMENT=en_US。UTF-8 LC_IDENTIFICATION=C ## ##附加的基本包:## [1]stats4 stats graphics grDevices utils datasets methods ##[8]基础## ##其他附加包:[1] SummarizedExperiment_1.28.0 GenomicRanges_1.50.0 ## [3] GenomeInfoDb_1.34.0 IRanges_2.32.0 ## [5] S4Vectors_0.36.0 MatrixGenerics_1.10.0 ## [9] Biobase_2.58.0 BiocGenerics_0.44.0 ## [11] BiocStyle_2.26.0 omicsPrint_1.18.0 ## [13] MASS_7.3-58.1 ## ##通过命名空间加载(并且没有附加):# # # # [1] lattice_0.20-45 tidyr_1.2.1 [3] assertthat_0.2.1 digest_0.6.30 # # [5] utf8_1.2.2 R6_2.5.1 # # [7] evaluate_0.17 highr_0.9 # # [9] pillar_1.8.1 zlibbioc_1.44.0 # # [11] rlang_1.0.6 data.table_1.14.4 # # [13] jquerylib_0.1.4 R.utils_2.12.1 # # [15] R.oo_1.25.0 Matrix_1.5-1 # # [17] rmarkdown_2.17 readr_2.1.3 # # [19] stringr_1.4.1 rcurl_1.98 - 1.9 # # [21] DelayedArray_0.24.0 compiler_4.2.1 # # [23] xfun_0.34 pkgconfig_2.0.3 # # [25] htmltools_0.5.3 tidyselect_1.2.0 # # [27] tibble_3.1.8[41] stringi_1.7.8 cachem_1.0.6 ## [43] XVector_0.38.0 MultiAssayExperiment_1.24.0 ## [45] limma_3.54.0 xml2_1.3.3 ## [47] bslib_0.4.0 ellipsis_0.3.2 ## [51] tools_4.2.1 RaggedExperiment_1.22.0 ## [53] glue_1.6.2 purrr_0.3.5 ## [55]hms_1.1.2 fastmap_1.1.0 ## [57] yaml_2.3.6 BiocManager_1.30.19 ## [59] knitr_1.40 sass_0.4.2