- html - 出于某种原因,IE8 对我的 Sass 文件中继承的 html5 CSS 不友好?
- JMeter 在响应断言中使用 span 标签的问题
- html - 在 :hover and :active? 上具有不同效果的 CSS 动画
- html - 相对于居中的 html 内容固定的 CSS 重复背景?
给定一个核苷酸序列,我正在编写一些 Julia 代码来生成(屏蔽的)kmer 计数的稀疏向量,我希望它尽可能快地运行。
这是我当前的实现,
using Distributions
using SparseArrays
function kmer_profile(seq, k, mask)
basis = [4^i for i in (k - 1):-1:0]
d = Dict('A'=>0, 'C'=>1, 'G'=>2, 'T'=>3)
kmer_dict = Dict{Int, Int32}(4^k=>0)
for n in 1:(length(seq) - length(mask) + 1)
kmer_hash = 1
j = 1
for i in 1:length(mask)
if mask[i]
kmer_hash += d[seq[n+i-1]] * basis[j]
j += 1
end
end
haskey(kmer_dict, kmer_hash) ? kmer_dict[kmer_hash] += 1 : kmer_dict[kmer_hash] = 1
end
return sparsevec(kmer_dict)
end
seq = join(sample(['A','C','G','T'], 1000000))
mask_str = "111111011111001111111111111110"
mask = BitArray([parse(Bool, string(m)) for m in split(mask_str, "")])
k = sum(mask)
@time kmer_profile(seq, k, mask)
这段代码在我的 M1 MacBook Pro 上运行大约需要 0.3 秒,有什么方法可以让它运行得更快吗?
函数kmer_profile
使用大小为length(mask)
的滑动窗口来计算每个掩码kmer出现在核苷酸序列中的次数。掩码是二进制序列,掩码kmer是在掩码为零的位置处掉落核苷酸的kmer。例如。 kmer ACGT
和掩码 1001
将生成掩码 kmer AT
。
为了生成 kmer 哈希,该函数将每个 kmer 视为基数 4 的数字,然后将其转换为(基数 10)64 位整数,用于索引到 kmer 向量。
k 的大小等于掩码字符串中 1 的数量,并且隐式限制为 31,以便 kmer 哈希可以适合 64 位整数类型。
最佳答案
有几种可能的优化可以使此代码更快。
首先,可以将 Dict 转换为数组,因为基于数组的索引比基于字典的索引更快,而且这里可以实现这一点,因为键是 ASCII 字符。 p>
此外,通过预先计算代码并将结果放入临时数组,可以提取一次序列代码,而不是length(mask)
次。
此外,基于mask
的条件和循环携带的依赖使事情变得很慢。事实上,处理器无法(轻松)预测该情况,导致其停顿几个周期。循环携带的依赖性使事情变得更糟,因为处理器在此停顿期间几乎无法执行其他指令。这个问题可以通过基于mask
和basis
预先计算因子来解决。结果是更快的无分支循环。
完成上述优化后,最大的瓶颈是sparsevec
。事实上,它也花费了最初实现的近一半时间!优化这一步很困难,但并非不可能。由于 Julia 实现中的随机访问,速度很慢。首先可以通过对键值对进行排序来加快速度。由于更适合缓存的执行,它的速度更快,并且还可以帮助处理器的预测单元。这是一个复杂的话题。有关其工作原理的更多详细信息,请阅读 Why is processing a sorted array faster than processing an unsorted array? .
这是最终优化的代码:
function kmer_profile_opt(seq, k, mask)
basis = [4^i for i in (k - 1):-1:0]
d = zeros(Int8, 128)
d[Int64('A')] = 0
d[Int64('C')] = 1
d[Int64('G')] = 2
d[Int64('T')] = 3
seq_codes = [d[Int8(e)] for e in seq]
j = 1
premult = zeros(Int64, length(mask))
for i in 1:length(mask)
if mask[i]
premult[i] = basis[j]
j += 1
end
end
kmer_dict = Dict{Int, Int32}(4^k=>0)
for n in 1:(length(seq) - length(mask) + 1)
kmer_hash = 1
j = 1
for i in 1:length(mask)
kmer_hash += seq_codes[n+i-1] * premult[i]
end
haskey(kmer_dict, kmer_hash) ? kmer_dict[kmer_hash] += 1 : kmer_dict[kmer_hash] = 1
end
sorted_kmer_pairs = sort(collect(kmer_dict))
sorted_kmer_keys = [e[1] for e in sorted_kmer_pairs]
sorted_kmer_values = [e[2] for e in sorted_kmer_pairs]
return sparsevec(sorted_kmer_keys, sorted_kmer_values)
end
此代码比我机器上的初始实现快两倍。很大一部分时间仍然花在排序算法上。
代码还可以进一步优化。一种方法是使用并行排序算法。另一种方法是用移位替换 premult[i]
乘法,假设 premult[i]
已修改为包含指数,则移位速度更快。我预计该代码比原始代码快大约 4 倍。主要瓶颈应该是大字典的创建。进一步提高其性能非常困难(尽管仍然有可能)。
关于performance - 从核苷酸序列生成 kmer 计数向量的最快方法 (Julia),我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/75263837/
我想用一个向量执行以下操作。 a = np.array(np.arange(0, 4, 1)) 我想得到一个乘法,结果是一个矩阵 | 0 1 2 3 4 -| - - - - - - - 0
正如标题所述,我正在尝试使用 gsub,其中我使用向量作为“模式”和“替换”。目前,我的代码如下所示: names(x1) names(x1) [1] "2110023264A.Ms.Amp
所以当我需要做一些线性代数时,我更容易将向量视为列向量。因此,我更喜欢 (n,1) 这样的形状。 形状 (n,) 和 (n,1) 之间是否存在显着的内存使用差异? 什么是首选方式? 以及如何将 (n,
我不明白为什么 seq() 可以根据元素中是否存在小数点输出不同的类,而 c() 总是创建一个 num向量,无论是否存在小数。 例如: seqDec <- seq(1, 2, 0.5) # num v
机器学习与传统编程的一个重要区别在于机器学习比传统编程涉及了更多的数学知识。不过,随着机器学习的飞速发展,各种框架应运而生,在数据分析等应用中使用机器学习时,使用现成的库和框架成为常态,似乎越来越不需
寻找有关如何将 RegEnable 用作向量的示例/建议。此外,我想控制输入和使能信号成为 Vector 中寄存器索引的函数。 首先,我如何声明 RegEnable() 的 Vector,其次如何迭代
假设我有一个包含变量名称的向量 v1,我想为每个变量分配一个值(存储在单独的向量中)。我如何在没有迭代的情况下做到这一点? v1 <- c("a","b","c") v2 <- c(1,2,3) 我想
R 提供了三种类型来存储同质对象列表:向量、矩阵 和数组。 据我所知: 向量是一维数组的特殊情况 矩阵是二维数组的特例 数组还可以具有任意维度级别(包括 1 和 2)。 在向量上使用一维数组和在矩阵上
我正在绕着numpy/scipy中的所有选项转圈。点积、乘法、matmul、tensordot、einsum 等 我想将一维向量与二维矩阵(这将是稀疏csr)相乘并对结果求和,这样我就有了一个一维向量
我是一个 IDL 用户,正在慢慢切换到 numpy/scipy,并且有一个操作我在 IDL 中非常经常做,但无法用 numpy 重现: IDL> a = [2., 4] IDL> b = [3., 5
在python计算机图形工具包中,有一个vec3类型用于表示三分量向量,但是我如何进行以下乘法: 三分量向量乘以其转置结果得到 3*3 矩阵,如下例所示: a = vec3(1,1,1) matrix
我正在构建一款小型太空射击游戏。当涉及到空间物理学时,我曾经遇到过数学问题。 用文字描述如下:有一个最大速度。因此,如果您全速行驶,您的飞船将在屏幕上一遍又一遍地移动,就像在旧的小行星游戏中一样。如果
我正在尝试在 python 中实现 Vector3 类。如果我用 c++ 或 c# 编写 Vector3 类,我会将 X、Y 和 Z 成员存储为 float ,但在 python 中,我读到鸭式是要走
我是 Spark 和 Scala 的新手,我正在尝试阅读有关 MLlib 的文档。 http://spark.apache.org/docs/1.4.0/mllib-data-types.html上的
我有一个包含四个逻辑向量的数据框, v1 , v2 , v3 , v4 是对还是错。我需要根据 boolean 向量的组合对数据帧的每一行进行分类(例如, "None" , "v1 only" , "
我正在创建一个可视化来说明主成分分析的工作原理,方法是绘制一些实际数据的特征值(为了说明的目的,我将子集化为二维)。 我想要来自 this fantastic PCA tutorial 的这两个图的组
我有以下排序向量: > v [1] -1 0 1 2 4 5 2 3 4 5 7 8 5 6 7 8 10 11 如何在不遍历整个向量的情况下删除 -1、0 和 11
有什么方法可以让 R 对向量和其他序列数据结构使用基于零的索引,例如在 C 和 python 中。 我们有一些代码在 C 中进行一些数值处理,我们正在考虑将其移植到 R 中以利用其先进的统计功能,但是
我有一个函数可以查询我的数据库中最近的 X 个条目,它返回一个 map 向量,如下所示: [{:itemID "item1" :category "stuff" :price 5} {:itemI
我有 ([[AA ww me bl qw 100] [AA ee rr aa aa 100] [AA qq rr aa aa 90]] [[CC ww me bl qw 100] [CC ee rr
我是一名优秀的程序员,十分优秀!