问题标签 [htslib]

htslib我的模块顶级目录中有一个已编译的 C 库，这是./configure; make在运行之前首先使用 python 外部构建的python setup.py install。我正在尝试使用htslib. 我有 cython 脚本，例如script.pyx它包装了支持的头文件cpp_header.h，但这些是通过安装脚本使用 C++ 编译的。python 模块的设置运行没有问题，但在运行时htslib符号丢失，我得到ImportError: foo.cpython-37m-x86_64-linux-gnu.so: undefined symbol: hts_close. 我的目录结构是：

在我的设置脚本中，我列出了一些需要的库htslib，一些库library_dirs，include_dirs以及extra_compile_args构建Extension：

此外，我还添加htslib了以下package_data选项setup：

在cpp_header.h文件中，htslib 标头包含为：

我在unix上尝试这个。有人知道我怎样才能让这个工作吗？

我正在使用 c++ 处理 VCF 文件，为此我使用来自 htslib 的 vcf 库（https://github.com/samtools/htslib/blob/develop/htslib/vcf.h）。我知道可能有一些更好的库，但我也在使用 htslib 也有库的其他文件格式，所以我想坚持使用 htslib。

我找到了一些代码示例，可以在文件中打开读取并创建正确的结构。标头并在此处使用 VCF 文件中的一些信息：https ://gist.github.com/gatoravi/cad922bdf2b625a91126和http://wresch.github.io/2014/11/18/process-vcf-file-with -htslib.html

但是，如果我们坚持第一个示例，我已经将代码“解码”为以下代码，并附上我对代码的注释：

在上面的这段代码中，我在 while 循环中对多个 std::cout 进行了注释，以通过我的注释清楚地说明某些功能是什么——即“摆脱”是染色体。据我所知，vcf 库的名称“rid”或“nfmt”都是预定义的。运行此代码，我可以打印多个内容，例如染色体名称、位置等。但我有几个问题：

我的 VCF 文件具有 #CHROM POS ID REF ALT QUAL FILTER INFO FORMAT 的一般结构，其中有几行仅显示前 6 列的小示例：

我的问题是在使用库时我希望打印出参考（第 4 列）和替代（第 5 列），所以对于第 1 行：REF = A & ALT = G，对于第 5 行：REF = T & ALT = TATGTTATG。

谁能帮助我准确理解提取这两个字段需要做什么？我在库描述中看不到如何使用“test_record->”来提取这些？

我希望我的问题有点道理。感谢您的时间和帮助。

I am using htslib for extracting all the information contained in a VCF file in C++.

Currently, thanks to the VCF specification and the documentation in the file vcf.h, I have successfully extracted all the metadata information in the header (Meta-Information Lines), and most of the information contained in each row of the body of the file (Data Lines).

However, I don't know how to extract the genotype information (sample columns).

I am using example files from the 1000G project. This is an example of two rows of the file, it shows the Format field and two samples (The file has more than 1000 samples per each row, I would like to extract the data for all of them):

I know that this is a heavy task that would take some computation time. I have extracted the names of each column (HG00096, HG00077...), but I don't know how to extract the information of each sample either as a full string (e.g., "0|0:0.050:-0.48,-0.48,-0.48"), as a set (array, map, vector...) of key-value pairs (e.g., [("GT", "0|0"), ("DS", "0.050"), ("GL", "-0.48,-0.48,-0.48")), or simply as an array of values (e.g., ["0|0", "0.050", "-0.48,-0.48,-0.48"]. I would like to do this for each sample.

I have been reading the documentation in the vcf.h file and I think that the function bcf_get_genotypes(hdr,line,dst,ndst) may be suitable for this, but I don't know for sure how to use it for extracting the values as strings. Also, I think that this information may be stored inside the 'p' pointer of 'bcf_fmt_t', but I don't know for sure, it just contains an array set of uint8_t values and I don't know if a string (or char array) can be extracted in the way I want.

Is there a way of doing that I am trying to do?

我无法使用以下任一命令安装带有 conda 的 htslib v1.12：

使用conda install -c bioconda/label/cf201901 htslib给了我 htslib v1.9。

有谁知道如何使用 conda 安装 v1.12？谢谢！

我正在使用 Rcpp 将 C++ 程序包装在 R 包中。我的 C++ 程序需要以下标头：

在编译之前，我通常在 Ubuntu 中加载以下模块：

我通常使用 GCC/9.3.0 在 Ubuntu 中使用以下标志编译 C++ 脚本：

由于我是通过 Rcpp 从 R 访问程序，所以我不知道如何加载 HTSlib 模块。当我尝试“清理并重建”包时，我收到以下错误：

我有两个问题：

1. 从 C++ 源代码构建 R 包时如何加载 C++ 模块？

2. 从 C++ 源代码构建 R 包时如何包含编译标志？

我创建了一个最小的头文件、.R 文件和 C++ 源脚本。该脚本打开一个 bam 文件并输出读取的染色体名称和位置。这些文件并不代表我想要运行的实际程序（这里包含的程序太长太复杂），但是当我尝试使用 Rcpp 构建包时会产生相同的错误。

C++源文件：

头文件：

R文件：

R 文件位于包的“R”目录中，而 C++ 脚本和头文件位于“src”目录中。

我有一个 Rcpp 函数，它读取大型 BAM 文件（1-20GB，使用htslib）并创建几个非常长std::vector的 s（最多 80M 个元素）。在阅读之前不知道元素的数量，所以我不能使用Rcpp::IntegerVectorand Rcpp::CharacterVector。据我了解，当我Rcpp::wrap进一步使用它们时，会创建副本。在这种情况下，有没有办法加快数据从 C++ 到 R 的传输？是否有可以在 Rcpp 函数中创建的数据结构，对push_back元素尽可能快std::vector，并通过引用传递给 R？

以防万一，这是我目前创建它们的方式：

以下是我包装和退回它们的方式：

编辑 1（2021.10.19）：

感谢大家的意见，我需要更多的时间来检查是否stringfish可以使用，但我从 cpp11 包 vignettes 中运行了一个稍微修改的测试，以将其与std::vector. 这是代码和结果（表明它std::vector<int>仍然更快，尽管它必须Rcpp::wrap在返回时被 ped）：

编辑2：

std::vector<std::string>比在这些测试条件下稍慢cpp11::writable::strings，但更节省内存：

解决方案（2022.01.12）：

...对于那些有类似问题的人。在这种特殊情况下，我不需要std::vector在 R 中使用数据。因此XPtr很容易解决了我的问题，将 BAM 加载时间缩短了近两倍。创建指针：

然后存储为data.frame属性：

并在其他地方重复使用，如下所示：