如何计算R中每列中残基（核苷酸）覆盖率的频率？

Question

我有一个 fasta 文件，其数据集如下：

>sequence_1
ACCTGC--A
>sequence_2
ACC-GCTTA
>sequence_3
ACCTGCTTA

最终目标是产生指示核苷酸覆盖百分比的结果。因此，我会得到 100%、100%、100%、66.66%、100%、100%、66.66%、66.66$、100%。

我是 R 的新手，我现在的想法是将序列作为字符串读取，然后转换为列表，然后遍历所有列表的第一个元素等等。

有更简单的方法吗？

Answer 1

这是一个简单的解决方案 strsplit。

Biostrings 使阅读 .fasta 文件变得微不足道。请注意，此方法要求序列已经对齐。

#BiocManager::install("Biostrings")
library(Biostrings)
data <- readDNAStringSet("data.fasta")
split <- strsplit(as.character(data),"")
logical.mat <- sapply(split, function(x) x != "-")
rowSums(logical.mat) / ncol(logical.mat)
#[1] 1.0000000 1.0000000 1.0000000 0.6666667 1.0000000 1.0000000 0.6666667 0.6666667 1.0000000

Answer 2

没有使用 bio-packages 的经验..所以这里有一个 data.table 方法...

library( data.table )
DT <- data.table( text = c( "ACCTGC--A", "ACC-GCTTA", "ACCTGCTTA" ))

#         text
# 1: ACCTGC--A
# 2: ACC-GCTTA
# 3: ACCTGCTTA

#split text to nucleoid
#get maximum number of nucleoid
n_max = length( tstrsplit(gsub("(.)", "\1 ", DT$text), " ") )
#create a new column for each nucleoid
DT[, paste0( "nucl_", 1:n_max ) := tstrsplit(gsub("(.)", "\1 ", text), " ") ]
#         text nucl_1 nucl_2 nucl_3 nucl_4 nucl_5 nucl_6 nucl_7 nucl_8 nucl_9
# 1: ACCTGC--A      A      C      C      T      G      C      -      -      A
# 2: ACC-GCTTA      A      C      C      -      G      C      T      T      A
# 3: ACCTGCTTA      A      C      C      T      G      C      T      T      A
                           
#how to procees? possibly melt and summarise?
DT.melt <- melt(DT, id.vars = "text", measure.vars = patterns("^nucl") )
ans <- DT.melt[, .N, by = .(variable, value) ][, perc := N / sum(N), by = .(variable)][]

#     variable value N      perc
# 1:    nucl_1     A 3 1.0000000
# 2:    nucl_2     C 3 1.0000000
# 3:    nucl_3     C 3 1.0000000
# 4:    nucl_4     T 2 0.6666667
# 5:    nucl_4     - 1 0.3333333
# 6:    nucl_5     G 3 1.0000000
# 7:    nucl_6     C 3 1.0000000
# 8:    nucl_7     - 1 0.3333333
# 9:    nucl_7     T 2 0.6666667
# 10:   nucl_8     - 1 0.3333333
# 11:   nucl_8     T 2 0.6666667
# 12:   nucl_9     A 3 1.0000000

如何计算R中每列中残基（核苷酸）覆盖率的频率？

How to calculate frequency of residue(nucleotide) coverage in each column in R?

r

sequence

bioinformatics

fasta