3

给定 8 个样本 (A1-A8) 的以下数据：

A1 A2 A3 A4 A5 A6 A7 A8
T T T T T T T C 
T C T T T T T C
A A A G G A A A

每列是一个样本，每一行是一个标记，可能编码为 A、C、G、T，我希望计算每一行的 4 个等位基因中任何一个的起源的概率。例如，上面第 1 行数据的输出应该是

   A C G T
A1 0 0 0 1/7
A2 0 0 0 1/7 
A3 0 0 0 1/7
A4 0 0 0 1/7
A5 0 0 0 1/7
A6 0 0 0 1/7
A7 0 0 0 1/7
A8 0 1 0 0

由于第 1 行中有 7 个样本具有 T，因此每个样本的概率为 1/7。由于只有 A8 拥有 C，因此将 C 分配给 A8 的概率为 100%。对于第 3 行，输出应为

   A C G T
A1 1/6 0 0 0
A2 1/6 0 0 0 
A3 1/6 0 0 0
A4 1/2 0 0 0
A5 1/2 0 0 0
A6 1/6 0 0 0
A7 1/6 0 0 0
A8 1/6 0 0 0

总输出应该是 i 8x4 矩阵的列表，其中 i 等于行数。

一个可重做的代码示例是：

states <- c("A1","A2","A3","A4","A5","A6","A7","A8") # Define the names of the states
A1 <- c("T","T","A") # Set the alleles for state A1 across 3 SNPs
A2 <- c("T","C","A") # Set the alleles for state A2 across 3 SNPs
A3 <- c("T","T","A") # Set the alleles for state A3 across 3 SNPs
A4 <- c("T","T","G") # Set the alleles for state A4 across 3 SNPs
A5 <- c("T","T","G") # Set the alleles for state A5 across 3 SNPs
A6 <- c("T","T","A") # Set the alleles for state A6 across 3 SNPs
A7 <- c("T","T","A") # Set the alleles for state A7 across 3 SNPs
A8 <- c("C","C","A") # Set the alleles for state A8 across 3 SNPs
theemissionmatrix <- matrix(t(c(A1,A2,A3,A4,A5,A6,A7,A8)), 8, 3, byrow = TRUE) # Create an 8 x 3 matrix
rownames(theemissionmatrix) <- states
theemissionmatrix # Print out the data matrix
   [,1] [,2] [,3]
A1 "T"  "T"  "A" 
A2 "T"  "C"  "A" 
A3 "T"  "T"  "A" 
A4 "T"  "T"  "G" 
A5 "T"  "T"  "G" 
A6 "T"  "T"  "A" 
A7 "T"  "T"  "A" 
A8 "C"  "C"  "A" 

test <- cbind(theemissionmatrix[,1]=="A",theemissionmatrix[,1]=="C",theemissionmatrix[,1]=="G",theemissionmatrix[,1]=="T")
colnames(test) <- c("A","C","G","T")

test
    [,1]  [,2]  [,3]  [,4]
A1 FALSE FALSE FALSE  TRUE
A2 FALSE FALSE FALSE  TRUE
A3 FALSE FALSE FALSE  TRUE
A4 FALSE FALSE FALSE  TRUE
A5 FALSE FALSE FALSE  TRUE
A6 FALSE FALSE FALSE  TRUE
A7 FALSE FALSE FALSE  TRUE
A8 FALSE  TRUE FALSE FALSE

经过这一步，我不确定如何将每列的总计数相加并除以得到总概率。

4

3 回答 3

2

这里有一些替代方案（从@zx8754 的答案中的“df”开始）：

选项 1：melt+ table+prop.table

有NaN值时不能很好地显示。

library(reshape2)
dfL <- melt(as.matrix(df))
Levs <- c("A", "C", "G", "T") 
dfL$value <- factor(dfL$value, Levs) ## Just to be sure

prop.table(table(dfL[c(2, 3, 1)]), c(2, 3))
# , , Var1 = 1
# 
#     value
# Var2 A         C G         T
#   A1   0.0000000   0.1428571
#   A2   0.0000000   0.1428571
#   A3   0.0000000   0.1428571
#   A4   0.0000000   0.1428571
#   A5   0.0000000   0.1428571
#   A6   0.0000000   0.1428571
#   A7   0.0000000   0.1428571
#   A8   1.0000000   0.0000000
# 
# , , Var1 = 2
# 
#     value
# Var2 A         C G         T
#   A1   0.0000000   0.1666667
#   A2   0.5000000   0.0000000
# ..... OUTPUT TRUNCATED

选项 2：melt+ by+table

可以很容易地显示0在NaN其他地方。

dfL <- melt(as.matrix(df))
Levs <- c("A", "C", "G", "T") 
dfL$value <- factor(dfL$value, Levs) ## Just to be sure

by(dfL[-1], dfL[1], FUN = function(x) {
  A <- prop.table(table(x), 2)
  A[is.nan(A)] <- 0
  A
})
# Var1: 1
#     value
# Var2         A         C         G         T
#   A1 0.0000000 0.0000000 0.0000000 0.1428571
#   A2 0.0000000 0.0000000 0.0000000 0.1428571
#   A3 0.0000000 0.0000000 0.0000000 0.1428571
#   A4 0.0000000 0.0000000 0.0000000 0.1428571
#   A5 0.0000000 0.0000000 0.0000000 0.1428571
#   A6 0.0000000 0.0000000 0.0000000 0.1428571
#   A7 0.0000000 0.0000000 0.0000000 0.1428571
#   A8 0.0000000 1.0000000 0.0000000 0.0000000
# ------------------------------------------------------------------------ 
# Var1: 2
#     value
# Var2         A         C         G         T
#   A1 0.0000000 0.0000000 0.0000000 0.1666667
#   A2 0.0000000 0.5000000 0.0000000 0.0000000
# ..... OUTPUT TRUNCATED

选项 3：lapply+table在对数据进行一些重组后

完全坚持基础 R，这是另一种选择....

Levs <- c("A", "C", "G", "T")
out <- data.frame(N = names(df), t(df), row.names=NULL)
Rows <- setdiff(names(out), "N")
out[Rows] <- lapply(out[Rows], function(x) factor(x, Levs))
Tables <- lapply(seq_along(Rows), function(x) {
  A <- prop.table(table(out[, 1], out[, Rows[x]]), 2)
  A[is.nan(A)] <- 0
  A
})
于 2013-11-18T16:49:44.433 回答
2

这是一个基本的 R 方式，依赖于split、table和sweep：

res <- lapply(split(as.matrix(df), 1:nrow(df)), factor, levels=unique(unlist(df)))
lapply(res, function(row) sweep(sapply(levels(row), '==', row), 1, table(row)[row], FUN='/'))
于 2013-11-18T17:07:49.393 回答
1

尝试这个：

#data
df <- read.table(text="
A1 A2 A3 A4 A5 A6 A7 A8
T T T T T T T C 
T C T T T T T C
A A A G G A A A", header=TRUE, as.is=T)

#ACGT
allele <- c("A","C","G","T")

#get counts: loop samples loop alleles
lapply(1:nrow(df),function(sample){
  sapply(c("A","C","G","T"),
         function(x){
           p <- as.numeric(df[sample,]==x) / sum(df[sample,]==x)
           #check if it is `not a number`
           ifelse(is.nan(p),0,p)
           })
  })
于 2013-11-18T16:30:27.460 回答