Выполнение цикла для файлов разных типов с использованием bcftools и awk для разделения файлов - PullRequest
1 голос
/ 29 октября 2019

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

У меня есть 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

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

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

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

(например, person 1.vcf -> person1_MantaINS.txt, person1_MantaDEL.txt, person1_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 через удаленный терминал.

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

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

1 Ответ

0 голосов
/ 05 ноября 2019
doit() {
    vcf="$1"
    out="$2"
    T="$3"
    bcftools view "$vcf" |
       awk -v T=${T} '{split($3,a,/\:/);if(a[1]==T) print $0}'  > "$out"
}
export -f doit
parallel doit {1} {1.}_{2}.txt {2} \
:::: list_of_vcf_files \
::: MantaINS MantaINV MantaDEL MantaBND MantaDUP Canvas

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

...