Запуск цикла по нескольким типам файлов с использованием bcftools и awk для разделения файлов

Уважаемое сообщество переполнения стека,

У меня есть 100 файлов .VCF (тип файла txt). В столбце «ID» есть различные вызовы структурных вариантов:

MantaINS
MantaINV
MantaDEL
MantaBND
MantaDUP
Canvas:REF
Canvas:GAIN
Canvas:LOSS

(рядом с числом, например, MantaINS:00:13:467, Canvas:Gain:594:31:23 и т. д.)

Файлы выглядят так (но они намного больше, тысячи записей в файле)

CHROM POS ID REF ALT QUAL FILTER INFO FORMAT
1 2827693 MantaDEL:0:2:5000 CCGTGGATGCGGGGACCCGCATCCCCTCTCCCTTCACAGCTGAGTGACCCACATCCCCTCTCCCCTCGCA C . PASS SVTYPE=DEL;END=2827680;BKPTID=Pindel_LCS_D1099159;HOMLEN=1;HOMSEQ=C;SVLEN=-66 GT:GQ 1/1:13.9
2 321682  MantaBND:5:7:1:0 6 PASS IMPRECISE;SVTYPE=DEL;END=321887;SVLEN=-105;CIPOS=-56,20;CIEND=-10,62 GT:GQ 0/1:12
2 14477084 MantaINS:88:22:00:3 12 PASS IMPRECISE;SVTYPE=DEL;END=14477381;SVLEN=-297;MEINFO=AluYa5,5,307,+;CIPOS=-22,18;CIEND=-12,32 GT:GQ 0/1:12
3 9425916 MantaDEL:5:333:000 23 PASS IMPRECISE;SVTYPE=INS;END=9425916;SVLEN=6027;CIPOS=-16,22;MIINFO=L1HS,1,6025,- GT:GQ 1/1:15
3 2658945 MantaDUP:5:22:000 23 PASS IMPRECISE;SVTYPE=INS;END=9425916;SVLEN=6027;CIPOS=-16,22;MIINFO=L1HS,1,6025,- GT:GQ 1/1:15
6 1325462 MantaINV:3:000:000 23 PASS IMPRECISE;SVTYPE=INS;END=9425916;SVLEN=6027;CIPOS=-16,22;MIINFO=L1HS,1,6025,- GT:GQ 1/1:15
6 5783961 CavnasREF:7:943:1453 23 PASS IMPRECISE;SVTYPE=INS;END=9425916;SVLEN=6027;CIPOS=-16,22;MIINFO=L1HS,1,6025,- GT:GQ 1/1:15
7 9425916 CanvasGAIN:9:323:123 23 PASS IMPRECISE;SVTYPE=INS;END=9425916;SVLEN=6027;CIPOS=-16,22;MIINFO=L1HS,1,6025,- GT:GQ 1/1:15
8 9425916 CanvasLOSS:2:932:123 23 PASS IMPRECISE;SVTYPE=INS;END=9425916;SVLEN=6027;CIPOS=-16,22;MIINFO=L1HS,1,6025,- GT:GQ 1/1:15

Каждый файл находится в отдельной папке, и я создал текстовый файл с путями к файлам для всех 100 vcfs. Этот файл выглядит так (только первые 4):

genomes/by_date/2015-09-03/batch1/patient30/patient30.SV.vcf.gz   
genomes/by_date/2016-03-05/batch1/patient4/patient4.SV.vcf.gz    
genomes/by_date/2018-10-14/batch1/patient16/patient16.SV.vcf.gz   
genomes/by_date/2018-012-28/batch1/patient100/patient100.SV.vcf.gz
genomes/by_date/2018-03-14/batch1/patient1/patient1.SV.vcf.gz

Я хочу разделить файлы по типу структурного варианта, который находится в столбце идентификатора, поэтому для каждого входного файла vcf я получаю 8 выходных файлов, разделенных по типу идентификатора, например. для Manta_INS мне нужен файл .txt, в котором будет только строка ниже, взятая из приведенного выше примера:

2 14477084 MantaINS:88:22:00:3 12 PASS IMPRECISE;SVTYPE=DEL END=14477381 SVLEN=-297;MEINFO=AluYa5,5,307,+;CIPOS=-22,18;CIEND=-12,32 GT:GQ 0/1:12

т.е. для каждого входного vcf я хотел бы, чтобы на выходе было 8 разделенных файлов.

(например, человек 1.vcf -> человек1_MantaINS.txt, человек1_MantaDEL.txt, человек1_MantaINV.txt и т. д.)

В одном файле VCF я запустил:

for T in   MantaINS MantaINV MantaDEL MantaBND MantaDUP Canvas
do
   bcftools view person1.vcf  | awk -v T=${T} '{split($3,a,/\:/);if(a[1]==T) print $0}'  > ${T}.txt
done

Что отлично работает (кроме вызовов Canvas, в которых есть двоеточие). Однако я хочу ввести список из 100 файлов для запуска одного и того же цикла.

Я устал:

for T in   MantaINS MantaINV MantaDEL MantaBND MantaDUP Canvas:REF Canvas: GAIN Canvas:LOSS
    do
       parallel -j6 "bcftools view {}  | awk -v T=${T} '{split($3,a,/\:/);if(a[1]==T) print $0}'  > $basename{}.txt" :::: paths_to_files.txt
    done

Что дает мне сообщение об ошибке: «нет такого файла или каталога» из параллельной для любого из моих типов файлов.

Я работаю над HPC через удаленный терминал.

Ваша помощь будет принята с благодарностью.

Большое спасибо


loops awk vcf-vcard gnu-parallel bcftools
person tacrolimus    schedule 29.10.2019    source источник
comment
Ваше требование неясно, и, возможно, поэтому вы проголосовали за закрытие вашего вопроса. Пожалуйста, отредактируйте свой вопрос, чтобы показать один входной файл, а затем имена файлов, необходимые для вывода. Нам нужен точный результат, а не словесное описание ;-). Удачи.   -  person shellter    schedule 29.10.2019
comment
@shellter Большое спасибо за отзыв. Я обновил вопрос, чтобы попытаться быть более ясным.   -  person tacrolimus    schedule 29.10.2019
comment
Хм, все еще есть некоторые проблемы с вашим вопросом. Не говорите, что выдает мне сообщение об ошибке, но включите фактический текст сообщения об ошибке, которое вы получаете. Как и мы, мы не можем сказать, является ли ваша ошибка ошибкой оболочки или чем-то связана с bcftools или parallel. Возвращаясь к предыдущему вопросу в своем вопросе, вы упоминаете person1. Это предполагает, что есть person2 ... personN? Почему бы не обернуть for T внутри другого цикла for P in person1 ... personN? Я бы удалил все дополнительные попытки решить вашу проблему и оставил бы только ту версию, которая, по вашему мнению, действительно должна работать или является вашей лучшей попыткой....   -  person shellter    schedule 29.10.2019
comment
И действительно ли между вашими данными есть пустые строки? Если нет, не могли бы вы сократить ваши данные. Если в нем есть пустые строки, укажите это под примерными данными в виде комментария в скобках). Извините, я не пытаюсь вас обыграть, просто сделать проблему ясной для других — это само по себе небольшое искусство, так что это только обратная связь. Удачи.   -  person shellter    schedule 29.10.2019
comment
@shellter Я попытался сделать это яснее, как и просили. Нет, вы правы, эта обратная связь полезна, поскольку чем яснее я могу быть, тем выше шансы на успех. Я ценю, что более опытные члены этого сообщества, вероятно, устали от плохо заданных вопросов и т. Д. И т. Д. Большое спасибо.   -  person tacrolimus    schedule 29.10.2019
comment
Отлично! Если вы не хотите использовать цикл с переносом for P и должны использовать параллель, я рекомендую добавить тег для gnu-parallel. Мне нужно провести некоторое исследование, чтобы понять вашу строку parallel cmd, а сегодня я занят. ИЛИ я бы написал awk-скрипт, который фильтрует как awk -v T="$T" '/MantaINS/{print $0 > "MantaINS "T".out}' InFile . Это должно быть достаточно быстро. Удачи.   -  person shellter    schedule 29.10.2019
comment
Вы можете иметь все ваши MantaINS похожие цели в одном скрипте, и когда каждая строка инициирует совпадение, она будет записываться в соответствующий файл outFile. Поэкспериментируйте с одним, прежде чем возиться со всем проектом ;-). Нужно идти. Удачи.   -  person shellter    schedule 29.10.2019
comment
Ваш вызов parallel имеет только одну двойную кавычку (), поэтому что-то было введено неправильно. Обратите внимание, что переменные оболочки будут расширены внутри двойных кавычек до того, как параллель увидит их ($3, $0, $basename и т. д.).   -  person jhnc    schedule 30.10.2019
comment
@shellter спасибо, я посмотрю ваш скрипт   -  person tacrolimus    schedule 30.10.2019
comment
@jhnc извините, я добавил вторые кавычки - я работаю в шлюзовой камере, поэтому мне приходится вручную вводить любые проблемы, которые у меня возникают!   -  person tacrolimus    schedule 30.10.2019

Ответы (1)


arrow_upward
0
arrow_downward

Вы пишете, что это отлично работает для одного VCF-файла:

for T in   MantaINS MantaINV MantaDEL MantaBND MantaDUP Canvas
do
   bcftools view person1.vcf  | awk -v T=${T} '{split($3,a,/\:/);if(a[1]==T) print $0}'  > ${T}.txt
done

тогда это также должно работать:

doit() {
    vcf="$1"
    out="$2"
    T="$3"
    bcftools view "$vcf" |
       awk -v T=${T} '{split($3,a,/\:/);if(a[1]==T) print $0}'  > "$out"
}

for T in   MantaINS MantaINV MantaDEL MantaBND MantaDUP Canvas
do
   doit person1.vcf person1_${T}.txt ${T}
done

и если это работает, то это также должно работать:

export -f doit
parallel doit {1} {1.}_{2}.txt {2} \
:::: list_of_vcf_files \
::: MantaINS MantaINV MantaDEL MantaBND MantaDUP Canvas

Если это не то, что вы хотите, покажите 3 полных примера команд, которые вы хотите выполнить.

(Мне также неясно, какую именно команду вы хотите запустить для Canvas:GAIN, поэтому, пожалуйста, пусть это будет один из трех примеров).

person Ole Tange    schedule 05.11.2019
comment
Уважаемый Оле, Большое спасибо, что нашли время ответить на мой вопрос. Я попытался запустить скрипт, и я продолжаю получать сообщение об ошибке, связанное с отсутствием такого файла или каталога - ошибка, которую он дает, - это пути к файлам из list_of_vcf_files, добавленные к тегам MantaINS Manta INV и т. Д. Например, путь к файлу может быть genomes/18-11-10/HB39013/LI9034.sv.vcf, и теперь ошибка говорит: $mypath/genomes/18-11-10/HB39013/LI9034.sv.vcf_MantaINS.txt — такого файла нет или каталог. - person tacrolimus; 12.11.2019
comment
Файлы VCF, которые я хочу проанализировать, хранятся в отдельной папке только для чтения, куда я хотел бы, чтобы вывод был, если это вообще полезная информация? - person tacrolimus; 12.11.2019