r - 使用变量对象调用函数
问题描述
我有一个只有一列的多行数据框，该列具有可变长度的字符串，范围从 30000 到 200000 个字符（DNA 序列）。[以下是 150 个字符的示例]
TTCCCCAAACAGCAACTTTAAGGAGCAGCTTCCTTTATGATCCCTGATTGCCTCCCCTTTGTTCCCATAACAAGTAGTTTAAATTTTCTGTTAAAGTCCAAACCACATATTTACAATACCTCGCACC
这是完整的数据集：https ://drive.google.com/open?id=1f9prtKW5NnS-BLI5lqsl4FEi4PvRfxGR
我在 R 中有一个代码，它根据每行的长度将每行分成 20 个 bin，并计算每个 bin 的 G 和 C 的出现次数，并返回一个 20 列的矩阵。这是代码：
library(data.table)
data <- fread("string.fa", header = F)
loopchar <- function(data){ bins <- sapply(seq(1, nchar(data), nchar(data)/20), function(x) substr(data, x, x + nchar(data)/20 - 1))output <- (str_count(bins, c("G"))/nchar(bins) + str_count(bins, c("C"))/nchar(bins))*100}
result <- data.frame(t(apply(data,1,loopchar)))
但是，现在我想做一些不同的事情。而不是nchar(data)/20，我希望子字符串段 (20) 与我拥有的列表不同。所以现在对于我的数据框，第一行应该分为 22 个 bin/segments，代码是nchar(data)/22.
第二行应分为 21 个 bin，代码为nchar(data)/21，以此类推。我希望该函数不断更改数据的 bin 数量。我的带有字符串的数据数据框和带有 bin 的数字向量列表的长度相同。
做这个的最好方式是什么？
解决方案
将 Bioconductor 的一些库用于此类任务更为自然。就我而言，我使用Biostrings，但也许您可以找到另一种方法。
数据
您的文件太大，所以我创建了一个文本文件（在内存中），其中包含每一行的随机 DNA：
# set seed to create reproducible example
set.seed(53101614)
# create an example text file in memory
temp <- tempfile()
writeLines(
sapply(1:100, function(i){
paste(sample(c("A", "T", "C", "G"), sample(100:6000),
replace = T), collapse = "")
}),
con = temp
)
# read lines from tmp file
dna <- readLines(temp)
# unlink file
unlink(temp)
数据预处理
创建Biostrings::DNAStringSet对象
使用Biostrings::DNAStringSet()函数我们可以读取character向量来创建DNAStringSet对象。请注意，我假设所有记录都在标准 DNA 字母表中，即每个字符串仅包含A, T, C, G符号。如果它不适用于您的情况，请参阅Biostrings文档。
dna <- DNAStringSet(dna, use.names = F)
# inspect the output
dna
A DNAStringSet instance of length 100
width seq
[1] 2235 GGGCTTCCGGTGGTTGTAGGCCCATAAGGTGGGAAATATACA...GAAACGTCGACAAGATACAAACGAGTGGTCAACAGGCCAGCC
[2] 1507 ATGCGGTCTATCTACTTGTTCGGCCGAACCTTGAGGGCAGCC...AACGCTTTGTACCTGTCCCAGAGTCAGAAGTAACAGTTTAGC
[3] 1462 CATTGGAGTACATAGGGTATTCCCTCTCGTTGTATAACTCCA...TCCTACTTGCGAAGGCAGTCGCACACAAGGGTCTATTTCGTC
[4] 1440 ATGCTACGTTGGTAGGGTAACGCAGACTAGAACCACACGGGA...ATAAAGCCGTCACAAGGAATGTTAGCACTCAATGGCTCGCTA
[5] 3976 AAGCGGAAGTACACGTACCCGCGTAGATTACGTATAGTCGCC...TTACGCGTTGCTCAAATCGTTCGGTGCAGTTTTATAGTGATG
... ... ...
[96] 4924 AGTAAGCAGATCCAGAGTACTGTGAAAGACGTCAGATCCCGA...TATAATGGGTTGCGTGTTTGATTCTGCCATGAATCCTATGTT
[97] 5702 CCTGAAGAGGACGTTTCCCCCTACATCCAGTAGTATTGGTGT...TCTGCTTTGCGCGGCGGGGCCGGACTGTCCATGGCTCACTTG
[98] 5603 GCGGCTGATTATTGCCCGTCTGCCTGCATGATCGAGCAGAAC...CTCTTTACATGCTCATAGGAATCGGCAACGAAGGAGAGAGTC
[99] 3775 GGCAAGACGGTCAGATGTTTTGATGTCCGGGCGGATATCCTT...CGCTGCCCGTGACAATAGTTATCATAAGGAGACCTGGATGGT
[100] 407 TGTCGCAACCTCTCTTGCACGTCCAATTCCCCGACGGTTCTA...GCGACATTCCGGAGTCTGCGCAGCCTATGTATACCCTACAGA
创建随机 N 个 bin 的向量
set.seed(53101614)
k <- sample(100, 100, replace = T)
# inspect the output
head(k)
[1] 37 32 63 76 19 41
创建Views对象是由N = k[i]块表示的每个 DNA 序列
IRanges::Views使用容器解决您的问题要容易得多。这东西非常快速和美丽。
首先，我们将每个测序的 DNA 划分为多个k[i]范围：
seqviews <- lapply(seq_along(dna), function(i){
seq = dna[[i]]
seq_length = length(seq)
starts = seq(1, seq_length - seq_length %% k[i], seq_length %/% k[i])
Views(seq, start = starts, end = c(starts[-1] - 1, seq_length))
}
)
# inspect the output for k[2] and seqviews[2]
k[2]
seqviews[2]
32
Views on a 1507-letter DNAString subject subject: ATGCGGTCTATCTACTTG...GTCAGAAGTAACAGTTTAG
views:
start end width
[1] 1 47 47 [ATGCGGTCTATCTACTTGTTCGGCCGAACCTTGAGGGCAGCCAGCTA]
[2] 48 94 47 [ACCGCCGGAGACCTGAGTCCACCACACCCATTCGATCTCCATGGTTG]
[3] 95 141 47 [GCGCTCTCCGAGGTGCCACGTCAAGTTGTACTACTCTCTCAGACCTC]
[4] 142 188 47 [TTGTTAGAAGTCCCGAGGTATATGCGCAATACCTCAACCGAAGCGCC]
[5] 189 235 47 [TGATGAGCAAACGTTTCTTATAGTCGCGACCTTGTCCCGAGGACTTG]
... ... ... ... ...
[28] 1270 1316 47 [AGGCGAGGGCAGGGCACATGTTTCTACAGTGAGGCGTGATCCGCTCC]
[29] 1317 1363 47 [GAGGCAAGCTCGTGAACTGTCGTGGCAAGTTACTTATGAGGATGTCA]
[30] 1364 1410 47 [TGGGCAGATGCAACAGACTGCTATTGGCGGGAGAGAGGCATCGACAT]
[31] 1411 1457 47 [ACCGTCTCAAGTACCACAGCTGAGAGGCTCTCGTGGAGATGCGCACA]
[32] 1458 1507 50 [TGAGTCGTAACGCTTTGTACCTGTCCCAGAGTCAGAAGTAACAGTTTAGC]
之后，我们检查所有序列是否已被划分为所需数量的块：
all(sapply(seq_along(k), function(i) k[i] == length(seqviews[[i]])))
[1] TRUE
重要观察
在我们继续之前，有一个关于你的函数的重要观察。
您的函数生成 N 个长度可变的块（因为它生成的索引是float而不是integers，所以substr()当您调用它时，会将提供的索引四舍五入到最接近的整数。
例如，从集合中提取第一条记录，并使用您的dna代码将此序列拆分为 37 个 bin将产生以下结果：
dna_1 <- as.character(dna[[1]])
sprintf("DNA#1: %d bp long, 37 chunks", nchar(dna_1))
[1] "DNA#1: 2235 bp long, 37 chunks"
bins <- sapply(seq(1, nchar(dna_1), nchar(dna_1)/37),
function(x){
substr(dna_1, x, x + nchar(dna_1)/37 - 1)
}
)
bins_length <- sapply(bins, nchar)
barplot(table(bins_length),
xlab = "Bin's length",
ylab = "Count",
main = "Bin's length variability"
)
我在我的代码中使用的方法，虽然length(dna[[i]]) %% k[i] != 0（提醒），产生相等长度k[i] - 1的箱，只有最后一个箱的长度等于：length(dna[i]) %/% k[i] + length(dna[[i]] %% k[i]
bins_length <- sapply(seqviews, length)
barplot(table(bins_length),
xlab = "Bin's length",
ylab = "Count",
main = "Bin's length variability"
)
GC含量计算
如上所述，Biostrings::letterFrequency()应用于可以IRanges::Views让您轻松计算 GC 内容：
找出每个 DNA 序列中每个 bin 的 GC 频率
GC <- lapply(seqviews, letterFrequency, letters = "GC", as.prob = TRUE)
GC <- lapply(GC, "*", 100)
head(GC[[1]])
G|C
[1,] 53.33333
[2,] 46.66667
[3,] 50.00000
[4,] 55.00000
[5,] 60.00000
[6,] 45.00000
绘制 DNA 的 GC 含量1:9
par(mfrow = c(3, 3))
invisible(
lapply(1:9, function(i){
plot(GC[[i]],
type = "l",
main = sprintf("DNA #%d, %d bp, %d bins", i, length(dna[[i]]), k[i]),
xlab = "N bins",
ylab = "GC content, %",
ylim = c(0, 100)
)
abline(h = 50, lty = 2, col = "red")
}
)
)
推荐阅读
- spring-boot - 我已经在 Spring Boot 代码中实现了 JWT 令牌安全性。如何在我的代码中的任何地方获取 jwt 令牌？需要保存审核
- c - 如何替换字符串中的字符
- excel - 如何解码 Excel VBA 中选择的二维码？
- java - 如何在将对象添加到 HashSet 之前检查对象是否属于某个类？
- c# - 非对称加密和产品密钥
- r - 如何在不删除格式的情况下将格式化的“.sas7bdat”文件导入“R”？
- postman - How to create signup test that works every time without changing user signup data?
- c - I have to read a file and use Linked Lists to store the data in C
- mongodb - $lookup to an array in MongoDB
- c# - 使用对象和方法将文本存储到由分隔符分隔的列表中