Я написал сценарий, который принимает три пользовательских ввода:
- Имя текстового файла со строками, соответствующими позиции хромосомы и базовой пары, например
extract.txt. - Имя файла VCF для поиска значений в текстовом файле и извлечения, т.е.
test.vcf. - Желаемое выходное имя нового файла VCF, который содержит извлеченные строки из
test.vcf
Для тех, кто не знаком с файлами VCF, есть много сложной информации заголовка, но единственная строка заголовка, которая имеет значение в этом контексте проблемы, - это последняя строка заголовка, которая начинается с #CHROM. Игрушечный пример одного из этих файлов выглядит следующим образом (Примечание: в файле VCF эти столбцы разделены табуляцией, но я не уверен, как отформатировать это здесь):
#CHROM POS ID REF ALT QUAL FILTER INFO FORMAT
Chr06 5567134 . T C 999 PASS DP=13782;VDB=0.0302;AF1=0.7154;AC1=1312;DP4=1987,2142,4337,4552;MQ=39;FQ=999;PV4=0.49,1.5e-10,0,1 GT:PL:DP:SP:GQ
Chr06 5567140 . G A 999 PASS DP=13537;VDB=0.0304;AF1=0.7489;AC1=1374;DP4=1744,1858,4383,4621;MQ=39;FQ=999;PV4=0.8,0.068,0,1 GT:PL:DP:SP:GQ
Chr06 5567195 . G T 999 PASS DP=12016;VDB=0.0284;AF1=0.384;AC1=704;DP4=3311,4518,1537,1850;MQ=37;FQ=999;PV4=0.0026,0,0,1 GT:PL:DP:SP:GQ
Chr06 5567224 . TAAAAA TAGAAACAAAAA 999 PASS INDEL;DP=13190;VDB=0.0352;AF1=0.1229;G3=0.7854,0.1806,0.03405;HWE=1.71e-05;AC1=225;DP4=4930,5006,542,685;MQ=40;FQ=999;PV4=0.00035,1,0,1 GT:PL:DP:SP:GQ
Chr06 5567247 . T A 999 PASS DP=14383;VDB=0.03;AF1=0.1233;G3=0.7772,0.1986,0.02415;HWE=0.022;AC1=226;DP4=6484,6134,587,654;MQ=41;FQ=999;PV4=0.0062,9.3e-07,0,0.16 GT:PL:DP:SP:GQ
Chr06 5567444 . TAAAAAA TAAAAAAA 999 PASS INDEL;DP=10303;VDB=0.0235;AF1=0.6037;G3=0.1784,0.4375,0.3841;HWE=0.0162;AC1=1107;DP4=1996,2224,2158,2446;MQ=46;FQ=999;PV4=0.7,1,6.3e-30,1 GT:PL:DP:SP:GQ
Chr06 5567497 . AG A 999 PASS INDEL;DP=5243;VDB=0.028;AF1=0.07142;G3=0.873,0.08398,0.04299;HWE=0.000541;AC1=131;DP4=2297,2010,0,146;MQ=46;FQ=999;PV4=0,0,0.00097,0 GT:PL:DP:SP:GQ
Chr06 5567499 . TAAA TGGCCAAATAAGCCACTCAAAAGAAATACAGCCAAAAACATCTACAAA 999 PASS INDEL;DP=5243;VDB=0.0273;AF1=0.3508;G3=0.4662,0.2686,0.2652;HWE=1.95e-12;AC1=643;DP4=1993,1638,195,545;MQ=46;FQ=999;PV4=0,1,5.1e-18,0 GT:PL:DP:SP:GQ
Chr06 5567583 . G C 999 PASS DP=12372;VDB=0.0279;AF1=0.09276;AC1=170;DP4=4794,5928,512,569;MQ=42;FQ=999;PV4=0.095,1,5.5e-31,1 GT:PL:DP:SP:GQ
Chr06 5567628 . G T 999 PASS DP=12657;VDB=0.0244;AF1=0.1049;AC1=192;DP4=5230,6194,197,578;MQ=40;FQ=999;PV4=1.2e-29,9.8e-06,0,1 GT:PL:DP:SP:GQ
Игрушечный пример txt файл элементов для извлечения выглядит так:
Chr06\t5567140
Chr06\t5567583
Chr06\t5567224
Я создал этот текстовый файл, чтобы явно содержать разделитель VCF (табуляция) между значениями хромосомы и позиции.
У меня есть bash скрипт, который работает, когда я тестирую его на своем локальном компьютере OSX.
Вот мой bash скрипт:
#!/bin/bash
# Reset getopts
OPTIND=1
usage="$(basename "$0") [-h] [-i] [-v] [-o] -- program to extract loci from VCF files
where:
-h show this help text
-i input text file containing loci to extract
CHR and BP columns must be separated by appropriate delimiter
-v VCF file to extract loci from
-o output file name
"
# Read in variables
while getopts ":hi:v:o:" opt; do
case "$opt" in
h) echo "$usage"
exit
;;
i)
input=$OPTARG
;;
v)
vcf=$OPTARG
;;
o)
output=$OPTARG
;;
esac
done
echo
echo "File containing SNP positions: '$input'"
echo "VCF file to extract from: '$vcf'"
echo "Output filename: '$output'"
LC_ALL=C grep '#CHROM' $vcf > $output
while IFS= read -r line
do
LC_ALL=C grep "$line" $vcf >> $output
done < $input
Когда я загружаю этот файл на свой Linux сервер на бегать в масштабе он ломается. Команда grep в while l oop, т.е. LC_ALL=C grep "$line" $vcf >> $output не может найти что-либо в более крупном файле VCF.
Устранение неполадок. Я выяснил, что сценарий работает на сервере, если я изменю текстовый файл так, чтобы записи больше не были Chr06\t5567140, а теперь просто 5567140. Это просто доказывает, что записи содержатся в файле, который я ищу, но что-то не удается расширить шаблон CHR-tab-BP для соответствия.
Очевидно, я мог бы просто изменить все свои текстовые файлы и покончить с проблема, но это сделает мой инструмент менее универсальным, и теперь мне просто любопытно, почему этот сценарий работал локально, но не работает на сервере.
Любая помощь или комментарии о том, как улучшить мой существующий код, приветствуются. Дайте мне знать, если что-то неясно по этому вопросу.
Спасибо!
РЕДАКТИРОВАТЬ:
Я получил скрипт, работающий. Мне просто пришлось изменить форматирование текстового файла, чтобы использовать фактический символ табуляции вместо буквального \t. Тот факт, что использование \t работало на OSX, но не на Linux, предполагает, что OSX grep может интерпретировать этот символ, а GNU grep - нет (как @Sundeep указал в комментариях ниже). Поместите это изменение здесь на случай, если у кого-то возникнет аналогичная проблема в будущем.