使用bcftools和awk对多个文件类型运行循环以细分文件。

Question

亲爱的堆栈溢出社区，

我有100个.VCF文件(一种txt文件)。在"ID“列中有不同的结构变体调用：

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

(旁边还有几个，如护身符:00:13:467，帆布:增益: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生成了一个文件路径的txt文件。这个文件看起来如下(仅前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个输出文件，按ID类型除以，例如对于Manta_INS，我想要一个.txt文件，它的下面一行将从上面的示例中提取：

2 14477084护卫舰: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个细分文件。

(例如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

它工作得非常好(除了画布调用中有一个冒号)。但是，我想输入一个包含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

这给了我一条错误消息：“没有这样的文件或目录”，在并行内部为我的任何文件类型。

我正在通过远程终端开发一台高性能计算机。

你的帮助将不胜感激。

非常感谢

Answer 1

您可以这样写，这对于一个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运行什么命令，所以请让它成为3个示例之一)。