Извлечение строк текста в зависимости от len() определенного столбца

Я пытаюсь написать простой скрипт для извлечения определенных данных из файла VCF, который отображает варианты в последовательностях генома.

Сценарий должен извлечь заголовок из файла, а также SNV, исключая при этом любые значения. Варианты отображаются в 2 столбцах, ALT и REF. Каждый столбец отделен пробелом. Инделс будет иметь 2 символа в ALT или REF, 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()
Источник

2 ответа

Решение
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()

Было два вопроса:

  1. Второй цикл никогда не выполняется, потому что вы уже дошли до конца файла VCF в первом цикле. Здесь вы можете увидеть, как начать новый цикл в том же текстовом файле.
  2. Вы не разделяли строку правильно, это разделитель табуляции.

Поэтому я настроил выполнение кода только с одним циклом и поместил вкладку в качестве параметра split.

Источник

Ответ, данный 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.

Источник
Другие вопросы по тегам