awk или другие инструменты биоинформатики для фильтрации vcf - PullRequest
1 голос
/ 29 января 2020

Я пытаюсь отфильтровать некоторые строки в файле vcf, вот пример строк:

1   10505   rs548419688 A   T   100 PASS    AC=1;AF=0.000199681;AN=5008;NS=2504;DP=9632;EAS_AF=0;AMR_AF=0;AFR_AF=0.0008;E
UR_AF=0;SAS_AF=0;AA=.|||;VT=SNP
1   10506   rs568405545 C   G   100 PASS    AC=1;AF=0.000199681;AN=5008;NS=2504;DP=9676;EAS_AF=0;AMR_AF=0;AFR_AF=0.0008;E
UR_AF=0;SAS_AF=0;AA=.|||;VT=SNP
1   10511   rs534229142 G   A   100 PASS    AC=1;AF=0.000199681;AN=5008;NS=2504;DP=9869;EAS_AF=0;AMR_AF=0.0014;AFR_AF=0;E
UR_AF=0;SAS_AF=0;AA=.|||;VT=SNP
1   10539   rs537182016 C   A   100 PASS    AC=3;AF=0.000599042;AN=5008;NS=2504;DP=9203;EAS_AF=0;AMR_AF=0.0014;AFR_AF=0;E
UR_AF=0.001;SAS_AF=0.001;AA=.|||;VT=SNP
1   10542   rs572818783 C   T   100 PASS    AC=1;AF=0.000199681;AN=5008;NS=2504;DP=9007;EAS_AF=0.001;AMR_AF=0;AFR_AF=0;EU
R_AF=0;SAS_AF=0;AA=.|||;VT=SNP

Скажем, я хочу извлечь строки с AMR_AF больше 0,5, но не могу вычислить узнать, как использовать регулярные выражения Awk для выполнения такой работы. Пробовал vcftools, но это не сработало.

Ответы [ 3 ]

1 голос
/ 29 января 2020

Вы можете разбить строку в выбранном вами поле и проверить, больше ли значение элемента c элемента сразу после разбиения, чем ваш порог.

Более подробно, разбив вход yes,foo=2,bar=0.23,baz=1 на ,bar= даст массив, содержащий yes,foo=2 и 0.23,baz=1. В Awk, если вы сравните второй элемент с 0.2, он будет просто конвертировать столько, сколько может из начала значения, в число, а затем выполнить сравнение чисел c.

Таким образом

awk '{ split($0, x, /[\t;]AMR_AF=/) } x[2]>0.5' file.vcf

должен делать то, что вы хотите. Мы разбиваем строку на x и исследуем значение чисел c, равное x[2].

. [\t;] в регулярном выражении допускает либо табуляцию, либо точку с запятой перед именем поля; чтобы быть совершенно общим, возможно, вам даже следует использовать (^|[\t;]), чтобы также разрешить совпадение в начале строки.

Если вы хотите параметризовать это, возможно, попробуйте

awk -v field="AMR_AF" -v thres=0.5 '{ split($0, x, "(^|[\t;])" field "=")) } x[2]>thres' file.vcf

Recall что Awk обрабатывает скрипт для каждой строки ввода сверху вниз, где каждый оператор скрипта имеет форму

[ условие ] [{ action } ]

Как показывают квадратные скобки, обе части являются необязательными - если условие отсутствует, действие выполняется безоговорочно; если действие отсутствует, по умолчанию используется { print $0 }. Таким образом, наш скрипт сначала безоговорочно разделит строку, а затем условно напечатает ее, если x[2] больше порога.

GNU Awk может разбивать многосимвольный разделитель полей, поэтому вы также можете использовать -F '[\t;]AMR_AF=' .

awk -F '[\t;]AMR_AF=' '$2>0.5' file.vcf
1 голос
/ 29 января 2020

Не могли бы вы попробовать следующее.

awk 'match($0,/AMR_AF=[0-9]+\.[0-9]+|AMR_AF=[0-9]+/) && substr($0,RSTART+7,RLENGTH-7)>0.5'  Input_file

Объяснение: Использование match функции awk для соответствия регулярному выражению AMR_AF= digits.digits ИЛИ AMR_AF=digits и всякий раз, когда это регулярное выражение получает совпадения в строке, оно устанавливает RSTART и RLENGTH переменные. && (условие AND), чтобы проверить, что значение подстроки от RSTART+7 до RLENGTH-7 больше 0,5, затем вывести эту строку.

0 голосов
/ 31 января 2020

Использование bcftools :

bcftools view -i 'INFO/AMR_AF > 0.5' myFile.vcf

Дополнительные параметры см. В руководствах bcftools .

...