根据特定列的 len() 提取文本行答案

【问题标题】：Extracting lines of text depending on the len() of a particular column根据特定列的 len() 提取文本行
【发布时间】：2019-02-21 15:49:24
【问题描述】：

我正在尝试编写一个简单的脚本来从 VCF 文件中提取特定数据，该文件显示基因组序列中的变体。

脚本需要从文件中提取标头以及 SNV，同时省略任何插入删除。变体显示在 2 列中，即 ALT 和 REF。每列由空格分隔。 Indels 在 ALT 或 REF 中将有 2 个字符，SNV 将始终有 1 个。

到目前为止，我已经提取了标题（始终以 ## 开头），但没有提取任何变体数据。

original_file = open('/home/user/Documents/NA12878.vcf', 'r')
extracted_file = open('NA12878_SNV.txt', 'w+')

for line in original_file:
   if '##' in line:
       extracted_file.write(line)

# Extract SNVs while omitting indels 
# Indels will have multiple entries in the REF or ALT column
# The ALT and REF columns appear at position 4 & 5 respectively

for line in original_file:
    ref = str.split()[3]
    alt = str.split()[4]
    if len(ref) == 1 and len(alt) == 1:
        extracted_file.write(line)

original_file.close()
extracted_file.close()

【问题讨论】：

您确定正确拆分数据吗？如果您的变量名为line，为什么还要使用str.split()？应该是line.split()
你目前得到什么输出与你想要什么输出？
我认为您是对的，数据可能没有正确拆分，我已将代码修改为line.split()。代码运行，但我没有得到我想要的输出。理想情况下，提取的文件将包含完整的标题，以及不包含 indel 的数据行。我当前的输出只是标题行，以 # 开头。 VCF 文件的其余部分包含不带# 字符的纯文本行。这些行不见了。
您可以在此处附上您的文件吗？
我已经把文件上传到这个谷歌驱动器，VCF 文件可以以纯文本形式打开。我也上传了我的输出。 https://drive.google.com/open?id=1kzFZOxliWmbCcezsmMfNt0EBUdcfdPud

标签： python bioinformatics genome vcf-variant-call-format

【解决方案1】：

original_file = open('NA12878.vcf', 'r')
extracted_file = open('NA12878_SNV.txt', 'w+')
i=0

for line in original_file:
    if '##' in line:
        extracted_file.write(line)
    else:
        ref = line.split('  ')[3]
        alt = line.split('  ')[4]
        if len(ref) == 1 and len(alt) == 1:
            extracted_file.write(line)

# Extract SNVs while omitting indels 
# Indels will have multiple entries in the REF or ALT column
# The ALT and REF columns appear at position 4 & 5 respectively

original_file.close()
extracted_file.close()

有两个问题：

第二个循环永远不会执行，因为在第一个循环中您已经到达 VCF 文件的末尾。您可以查看here 如何在同一个文本文件上重新开始新循环。
您没有正确分隔由制表符分隔的行。

因此我将代码设置为仅使用一个循环执行并将选项卡作为拆分参数。

【讨论】：

太棒了，太棒了！非常感谢！
我会使用re.split(r'\t+', line) 和import re。

【解决方案2】：

Adirmola 给出的答案很好，但是你可以通过应用一些修改来提高代码质量：

# Use "with" context managers to open files.
# The closing will be automatic, even in case of problems.
with open("NA12878.vcf", "r") as vcf_file, \
        open("NA12878_SNV.txt", "w") as snv_file:
    for line in vcf_file:
        # Since you have specific knowledge of the format
        # you can be more specific when identifying header lines
        if line[:2] == "##":
            snv_file.write(line)
        else:
            # You can avoid doing the splitting twice
            # with list unpacking (using _ for ignored fields)
            # https://www.python.org/dev/peps/pep-3132/
            [_, _, _, ref, alt, *_] = line.split("\t")  # "\t" represents a tab
            if len(ref) == 1 and len(alt) == 1:
                snv_file.write(line)

我在您的文件上使用 python 3.6 对此进行了测试，最终得到 554 个 SNV。此处使用的某些语法（尤其是列表解包）可能不适用于较旧的 python 版本。

【讨论】：