查找与其他几列相比唯一的列中的字母，分段计数并打印它们的位置

Question

我正在尝试为 Skov 方法 https://github.com/LauritsSkov/Introgression-detection 准备一个输入文件，以查看海豚种群中的幽灵祖先。

我有一个如下所示的数据集：

chr  pos ind0 ind1 ind2 ind3 ind4 ind5 ind6 ind7 ind8 ind9 ind10

MRVK01001299.1 972    C    C    T    N    C    C    T    N    N    C     C

MRVK01001299.1 973    G    G    G    N    G    G    G    N    N    G     G

MRVK01001299.1 997    C    T    T    T    T    T    T    T    T    T     T

MRVK01001299.1 999    A    T    T    N    T    T    T    T    T    T     T

MRVK01001299.1 1018   A    C    T    N    T    C    C    T    T    T     T

MRVK01001299.1 1086   A    T    T    T    T    T    T    T    T    T     T

MRVK01001299.1 2125   C    C    T    N    C    C    T    N    N    C     C

MRVK01001299.1 2456   G    G    G    N    G    G    G    N    N    G     G

数据是在每个个体中发现的支架、位置和等位基因（A、T、C、G、N）。我正在将一个特定的个体 (ind0) 与其他个体进行比较，以查看他是否具有私人等位基因簇。如果是这样，这可能表明他的祖先来自我们没有抽样的人群（幽灵人群）。 chrm 和 pos 列表示宽吻海豚参考基因组中的位置。

这是dput() 在我的数据集的前 8 行（共 43,500 行）中针对前 3 个人（有 22 个人）的结果：

structure(list(chr = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = "MRVK01001299.1", class = "factor"), pos = c(1972L, 1973L, 2097L, 8281L, 8357L, 8653L, 8746L, 9255L), ind0 = structure(c(2L, 3L, 2L, 5L, 2L, 1L, 2L, 2L), .Label = c("A", "C", "G", "N", "T" ), class = "factor"), ind1 = structure(c(2L, 3L, 2L, 5L, 2L, 1L, 2L, 2L), .Label = c("A", "C", "G", "N", "T"), class = "factor"), ind2 = structure(c(5L, 3L, 5L, 5L, 5L, 1L, 5L, 2L), .Label = c("A", "C", "G", "N", "T"), class = "factor")), class = "data.frame", row.names = c(NA, 8L)) 

请注意，上表是包含私有等位基因的示例，但第一行不包含任何等位基因。

我想确定在 ind0 中唯一找到字母的位置 (pos)。 “N”不会被算作不同的字母。例如，我们将为位置 997、999 和 1086 设置一个唯一值。

然后，我想数一下ind0 有多少次在位置 (pos) 列中有 1000 个系列的私人信件。并打印每个箱的私有等位基因的位置。所以这将是（需要输出文件）：

0 2 997, 999

1000 1 1086

2000 0

etc

因为我们有两个位置（位置 997、999），ind0 的唯一值在 0 和 1000 之间，1 在 1000 和 2000 之间（位置 1086），0 在 2000 和 3000 之间。最远的值将高于 20,000,000 .

@zx87754 提供了以下代码，我想在私有等位基因的位置添加“print”。我在论坛上没有找到帮助，除了 print 或 awk 功能在这里可能会有所帮助。

确定哪些 SNP ind0 具有私有等位基因 - “N”不作为缺失数据算作不同的等位基因：

res1 <- df1[ rowSums(df1$ind0 == df1 [, -c(1:3)]) == 0 &
               apply(df1 [, -c(1:3)], 1, function(i) length(unique(i[ i != "N" ]))) == 1, ]

Count the number of positions with a private allele for ind0 in bins of 1000 bp:

res2<-table(cut(res1$pos, c(seq(0, 23092000, by=1000))))

非常感谢，

Answer 1

#DATA
set.seed(42)
d = data.frame(chr = "A",
               pos = 1:20,
               ind0 = sample(x = c("A", "C", "T", "G", "N"),
                             size = 20,
                             replace = TRUE,
                             prob = c(.24, .24, .24, .24, .04)))
d$ind1 = d$ind0
d$ind1[sample(1:20, 5)] = sample(x = c("A", "C", "T", "G", "N"),
                                 size = 5,
                                 replace = TRUE,
                                 prob = c(.24, .24, .24, .24, .04))
d$ind2 = d$ind1


# Subset the columns ind1 thorugh the end and get the unique value for each row
all_rows = apply(X = d[!colnames(d) %in% c("chr", "pos", "ind0")],
                 MARGIN = 1,
                 FUN = function(x){
                     if (length(unique(x)) == 1 & x[1] != "N") {
                         x[1]
                     } else {
                         NA
                     }
                 })

# Identify the row index where ind0 differs from all_rows
inds = which(d[, "ind0"] != all_rows & !is.na(d[, "ind0"]))

# Use the inds to extract the data that you need
d$pos[inds]
#> [1]  5 19 20
d[inds,]
#>    chr pos ind0 ind1 ind2
#> 5    A   5    G    A    A
#> 19   A  19    T    C    C
#> 20   A  20    G    A    A

^{由reprex package (v0.3.0) 于 2019 年 8 月 30 日创建}

对于分箱和计数，您可以这样做

positions = d$pos[inds]
bins = seq(0, 20, 5)
ans = aggregate(list(pos = positions), list(b = cut(positions, bins)), c)
ans$n = lengths(ans$pos)
ans
#        b    pos n
#1   (0,5]      5 1
#2 (15,20] 19, 20 2