问题标签 [vcf-variant-call-format]

For questions regarding programming in ECMAScript (JavaScript/JS) and its various dialects/implementations (excluding ActionScript). Note JavaScript is NOT the same as Java! Please include all relevant tags on your question; e.g., [node.js], [jquery], [json], [reactjs], [angular], [ember.js], [vue.js], [typescript], [svelte], etc.

0 投票
0 回答
64 浏览

python - AttributeError:“列表”对象没有属性“拆分”vcf

我正在编写一个 python 代码来修改一个 vcf 文件,但是我有几个小时无法修复的错误。

以下是涉及的部分代码:

终端输出:

返回的错误:

谢谢你。谢谢你。

0 投票
1 回答
74 浏览

google-cloud-platform - GRCh37 WGS 上的 Google DeepVariant 管道,外显子组模型未完成

我有一个 hg19 对齐的 BAM,我希望为其生成一个 DeepVariant VCF。我使用 samtools 提取标题并确保 hg19 参考 FASTA 索引包含相同的重叠群和位置。我最初的目标是在这个 WGS BAM 上只运行一个外显子组模型,使用以下模型和区域:

MODEL=gs://deepvariant/models/DeepVariant/0.7.2/DeepVariant-inception_v3-0.7.2+data-wes_standard

--regions gs://deepvariant/exome-case-study-testdata/refseq.coding_exons.b37.extended50.bed

不幸的是,脚本抗议说 BED 和 BAM / FASTA 参考之间的匹配为 0。我决定运行相同的外显子组模型,但没有指定区域。这是我的脚本:

BAM 有对应的 BAI,FA 有 FAI 文件。DeepVariant QuickStart 表明这些设置将在 1-2 小时内生成 VCF,但我的管道现在已经运行了 7 多个小时。暂存文件夹现在有一个 call_variants,其中包含 32 个 GZ 文件中的 31 个。Genomics 管道视图显示了 11 个运行 call_variant 的管道,因此我怀疑它正在处理最后一个文件,准备将所有文件合并到一个 VCF 中。

我只是不明白为什么要花这么长时间。我故意排除了抢占式实例,并且文档说外显子组管道应该只需要 20 分钟(WGS 需要 1-2 小时)。为什么会这么慢?

0 投票
1 回答
216 浏览

python - 访问配置列表元素的前缀以从字典中获取值

所以我正在尝试使用snakemake构建管道,但我在访问配置文件中字典“small_reference”中的值时遇到了一些问题。根据示例,我想使用不同的参考进行对齐。

配置文件:

根据样本的前 3 个字符,我想选择不同的参考。当 config["samples"] 只是一个字符串时,我编写了一个可以解决问题的函数。但现在我想考虑运行文件夹,所以我有一个样本列表。

工作功能:

当我刚刚更改配置文件时,我首先遇到了这个错误Duplicate output file pattern in rule(所以在运行完整管道时)

测试规则:

现在运行测试规则时出现此错误:

示例:C014038p --> C01 --> /mnt/storage/projects/hiv_data/refs/BRCA/BRCA12_PALB2.fasta

如何根据管道正在分析的样本的前缀获得正确的“small_reference”?

0 投票
0 回答
35 浏览

linux - 循环文件中的行并替换linux命令行中的值

当我有多个具有相似扩展名的文件时,例如:

我还有另一个文件:

我想在 linux 中运行一个命令行,它对 test.vcf 文件执行操作,并且一直使用扩展名为 .txt 的 4 个文件之一,我可以遍历文件并只写一个命令而不是写4个不同的命令:

此命令行将在 test.vcf 上使用 files1.txt 并输出 test_output1.vcf 等,因此我可以使用 $i 替换文件的编号。

现在假设我只有一个文件“file.txt”,并且这个文件有多个列:

我想循环遍历列 V2 中的每个值,类似地,而不是编写 4 个命令行,我可以只编写一个命令。

对此有何建议?谢谢

0 投票
3 回答
68 浏览

awk - 列之间是否有用于字符串减法的 Linux 命令?

我正在将一些 SNP 列处理成 vcf 格式。

输入列如下:

预期输出:

0 投票
1 回答
319 浏览

r - 未来的系统命令/Rshiny 中的承诺

我在闪亮的应用程序中有以下 server.R 代码,其中系统命令在未来运行,它提供了一个output.vcf文件。创建此文件后,进度条将被删除并运行第二个系统命令以转换out.vcfout.txt

使用系统命令是因为 R 无法在 32Gb 机器上读取巨大的向量。因此,一些系统命令用于处理数据。

第一个系统命令中产生的输出,即out.vcf必须渲染到downloadHandler ,第二个命令的输出out.txt必须返回到renderDataTable

有人可以提出一种有效的方法吗?可能在 内部 运行两个系统命令future()并将输出返回到downloadHandlerrenderDataTable

0 投票
1 回答
86 浏览

bash - 自动 awk 命令删除 VCF 的 INFO 列中的重复字段,保持第一次出现

我的 VCF 文件中存在重复字段会导致其他程序出现问题。VCF 文件是制表符分隔的文件。下面给出了 INFO 列的一个单元格。细胞的结构是:

info1=x;info2=y;.....;info1=.;info2=.;.....

我需要一个脚本来删除第一次出现后重复的重复字段并将单元格替换为:

info1=x;info2=y;.....

实际细胞...

0 投票
2 回答
616 浏览

genetics - 将文件输入 plink 时遇到问题

我只是将文件输入 plink 时遇到问题(我对编码很陌生)。目前,我只是想在plink中制作一个床文件

我最近运行了一个 HD SNP Array 来寻找拷贝数变异。在我尝试将 ped 和 map 文件输入到 plink 之前,一切似乎都很好。最初,我收到“.ped 文件中缺少一半调用”的错误消息,但是当我检查文件行时,我看不到任何问题,即在错误消息中注明的位置,有一个 0/​​0。

经过多次尝试解决这个问题,我被建议尝试使用 vcf 格式的文件并将这些文件输入 plink。我现在收到一条不同的错误消息:

.vcf 文件第 533890 行的变体 bp 坐标无效

我在整个互联网上进行了搜索,但无法找到任何解决方案

这是我试图运行的代码行:

我知道我想制作一个床文件,而且我很确定我没有犯任何语法错误。我当前的错误消息如下:

检测到 257665 MB RAM;为主工作区保留 128832 MB --vcf:533k 变体已完成。错误:.vcf 文件第 533890 行的变体 bp 坐标无效。

任何帮助将不胜感激!

0 投票
1 回答
85 浏览

string - bash 文件子字符串仅在匹配项中添加/替换

我有一个包含不同条目的VCF 文件,我需要根据多个匹配替换(或添加如果不存在)子字符串。例如

我需要查看一些条件,以便为变异患者的每个组合仅隔离一条特定的行(两者都可以重复,但它们的组合是唯一的),例如,以第 4 行为例:

  • $2 == "47746675" && $3 == "C" && $4 == "G"
  • pz_name=GHARTxI16uuT15921

在这个特定的行中,我需要:

  • status=something;缺席时添加
  • 如果存在则替换status=<something-else>status=something

如何在 bash 中使用某种就地替换来完成这一切?可能吗?替代品性能有效的方法建议将不胜感激!

非常感谢您的帮助!

0 投票
3 回答
633 浏览

gzip - 如何在 Nim 中打开和读取 .gz 文件(最好是逐行)

我只是坐下来编写我的第一个Nim脚本来解析.vcf(变体调用格式)文件。这种文件格式存储来自测序数据的基因突变。

对于脚本语言,我PerlPythonNim. 我意识到我Nim还很年轻,但我什至找不到一个关于如何打开和读取.gz( gzip) 文件的清晰示例(最好是逐行)。

谁能提供一个简单的例子来逐行打开和读取gzip文件?Nim

Python,我习惯了以下(超级简单)代码:

我看过相关的问题,但不清楚。这些帖子也相对较旧,我希望/怀疑有更好的事情发生。这是我发现的:

  1. 逐行读取 gzip 压缩文件
  2. 文件、文件流和 GZFileStream
  3. 从 Nim 中的 tar.gz 存档中读取文件

真的很感激。

PS我也认为如果有人Nim在 StackOverflow 中创建标签会很有用。我没有创建标签的声誉。