主页 / unix / 问题 / 746286
Accepted
Anna1364
Asked: 2023-05-19 01:21:24 +0800 CST

根据某些条件过滤行

  • 772

我有一些 vcf 文件，我想过滤掉一些变体。这只是我的 vcf 的一小部分：文件开头有一些标题行（以 ## 开头），然后是变体（每个变体一行）。

##fileformat=VCFv4.2
##source=combiSV-v2.2
##fileDate=Mon May  8 11:32:53 2023
##contig=<ID=chrM,length=16571>
##contig=<ID=chr1,length=249250621>    
##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record">
##INFO=<ID=SVCALLERS,Number=.,Type=String,Description="SV callers that support this SV">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=DR,Number=1,Type=Integer,Description="# High-quality reference reads">
##FORMAT=<ID=DV,Number=1,Type=Integer,Description="# High-quality variant reads">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample
1   10862   id.1    N   <INS>   .   PASS    SVTYPE=INS;SVLEN=101;END=10862;SVCALLERS=cutesv,SVIM    GT:DR:DV    1/1:0:26
1   90258   id.2    N   <INS>   .   PASS    SVTYPE=INS;SVLEN=118;END=90258;SVCALLERS=SVIM,NanoSV    GT:DR:DV    1/1:0:9
1   90259   id.3    N   <INS>   .   PASS    SVTYPE=INS;SVLEN=36;END=90259;SVCALLERS=Sniffles    GT:DR:DV    0/1:44:7
1   185824  id.4    N   <DEL>   .   PASS    SVTYPE=DEL;SVLEN=80;END=186660;SVCALLERS=Sniffles,cutesv    GT:DR:DV    1/1:0:15
1   186241  id.5    N   <DEL>   .   PASS    SVTYPE=DEL;SVLEN=418;END=186662;SVCALLERS=SVIM,NanoSV   GT:DR:DV    1/1:2:12
1   526111  id.6    N   <DEL>   .   PASS    SVTYPE=DEL;SVLEN=624;END=526735;SVCALLERS=Sniffles,cutesv   GT:DR:DV    0/1:8
2   91926078    id.3958 N   <BND>   .   PASS    SVTYPE=BND;SVLEN=.;END=;SVCALLERS=Sniffles,NanoSV   GT:DR:DV    0/1:60:15

在保留标题行的同时，我想删除 SVLEN < 100 的行和仅包含一个 SVCALLERS 的行（这是两个都必须满足的两个条件，换句话说，我只想保留 SVLEN > 100 和至少两个 SVCALLERS 的行). 另外有一些行 ALT 是 BND 并且 vcf 没有为这种类型的变体提供任何 SVLEN，如果该行包含 BND，我只想保留它，如果它被两个调用者支持。示例：我想删除此变体，因为 SVLEN 小于 100，并且只有一个 SVCALLERS 检测到它

SVTYPE=INS;SVLEN=36;END=90259;SVCALLERS=Sniffles    GT:DR:DV    0/1:44:7
    1   185824  id.4    N   <DEL>   .   PASS

或者这一行也是，虽然有两个调用者但是SVLEN小于100：

SVTYPE=DEL;SVLEN=80;END=186660;SVCALLERS=Sniffles,cutesv    GT:DR:DV    1/1:0:15
    1   186241  id.5    N   <DEL>   .   PASS

有简单的方法吗？谢谢

我的最终文件应如下所示：

##fileformat=VCFv4.2
##source=combiSV-v2.2
##fileDate=Mon May  8 11:32:53 2023
##contig=<ID=chrM,length=16571>
##contig=<ID=chr1,length=249250621>    
##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record">
##INFO=<ID=SVCALLERS,Number=.,Type=String,Description="SV callers that support this SV">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=DR,Number=1,Type=Integer,Description="# High-quality reference reads">
##FORMAT=<ID=DV,Number=1,Type=Integer,Description="# High-quality variant reads">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample
1   10862   id.1    N   <INS>   .   PASS    SVTYPE=INS;SVLEN=101;END=10862;SVCALLERS=cutesv,SVIM    GT:DR:DV    1/1:0:26
1   90258   id.2    N   <INS>   .   PASS    SVTYPE=INS;SVLEN=118;END=90258;SVCALLERS=SVIM,NanoSV    GT:DR:DV    1/1:0:9
1   186241  id.5    N   <DEL>   .   PASS    SVTYPE=DEL;SVLEN=418;END=186662;SVCALLERS=SVIM,NanoSV   GT:DR:DV    1/1:2:12
1   526111  id.6    N   <DEL>   .   PASS    SVTYPE=DEL;SVLEN=624;END=526735;SVCALLERS=Sniffles,cutesv   GT:DR:DV    0/1:8
2   91926078    id.3958 N   <BND>   .   PASS    SVTYPE=BND;SVLEN=.;END=;SVCALLERS=Sniffles,NanoSV   GT:DR:DV    0/1:60:15
text-processing

2 个回答

  1. Best Answer

    这是一个 perl 方法：

    $ perl -F'\t' -lane '
  if(/^#/){ print; next }; 
  $F[7] =~ /\bSVLEN=(\d+)/; 
  $svlen=$1; 
  $F[7] =~ /\bSVCALLERS=([^;]+)/; 
  @callers=split(/,/,$1); 
  print if $svlen > 100 && scalar(@callers) > 1' file.vcf 
##fileformat=VCFv4.2
##source=combiSV-v2.2
##fileDate=Mon May  8 11:32:53 2023
##contig=<ID=chrM,length=16571>
##contig=<ID=chr1,length=249250621>    
##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record">
##INFO=<ID=SVCALLERS,Number=.,Type=String,Description="SV callers that support this SV">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=DR,Number=1,Type=Integer,Description="# High-quality reference reads">
##FORMAT=<ID=DV,Number=1,Type=Integer,Description="# High-quality variant reads">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample
1   10862   id.1    N   <INS>   .   PASS    SVTYPE=INS;SVLEN=101;END=10862;SVCALLERS=cutesv,SVIM    GT:DR:DV    1/1:0:26
1   90258   id.2    N   <INS>   .   PASS    SVTYPE=INS;SVLEN=118;END=90258;SVCALLERS=SVIM,NanoSV    GT:DR:DV    1/1:0:9
1   186241  id.5    N   <DEL>   .   PASS    SVTYPE=DEL;SVLEN=418;END=186662;SVCALLERS=SVIM,NanoSV   GT:DR:DV    1/1:2:12
1   526111  id.6    N   <DEL>   .   PASS    SVTYPE=DEL;SVLEN=624;END=526735;SVCALLERS=Sniffles,cutesv   GT:DR:DV    0/1:8

    解释

    • perl -F'\t' -lane:-a使 perl 工作起来有点像awk它自动将每个输入行拆分为由-F（空白，默认情况下，如 awk）给出的字符并将其保存在数组中@F。因此，由于数组从 开始计数0，因此第一个字段是$F[0[，第二个字段$F[1]依此类推。接下来，从每个输入行中删除尾随换行符，并向每个调用-l添加一个，意思是“逐行读取输入文件并将给定的脚本应用于每一行。\nprint-n-e
    • if(/^#/){ print; next }; : 如果这是标题行，只需打印它并转到下一行。
    • $F[7] =~ /\bSVLEN=(\d+)/; $svlen=$1;: 匹配第8个字段之后最长的一串数字SVLEN=，保存为$svlen. 这将是长度。这确保我们只匹配单词边界，所以如果你的文件中\b有类似的东西，这不会失败。NOTSVLEN=
    • $F[7] =~ /\bSVCALLERS=([^;]+)/; @callers=split(/,/,$1);：现在在第 8 个字段中搜索字符串SVCALLERS=，取其后最长的一段非;字符，然后将其拆分,到数组中@callers。现在有用于此 CNV 的 CNV 呼叫者列表。
    • print if $svlen > 100 && scalar(@callers) > 1：如果长度超过 100 并且调用者的数量（scalar(@array)给出数组中的元素数量）超过 ，我们现在打印该行1。

    如果您更喜欢轻度打高尔夫球，那么这里有一个更简洁、更不清晰的基本内容：

    perl -F'\t' -lane '$F[7]=~/\bSVLEN=(\d+)/;$s=$1;$F[7]=~/\bSVCALLERS=([^;]+)/; /^#/ || ($s>100&&scalar(split(/,/,$1)) > 1) || next; print' file.vcf

    如果你还想保持没有 SVLEN 的线路，只要他们至少有两个呼叫者，使用这个：

    $ perl -F'\t' -lane 'if(/^#/){ print; next }; $F[7] =~ /\bSVLEN=([.\d]+)/; $svlen=$1; $F[7] =~ /\bSVCALLERS=([^;]+)/; next unless ($svlen > 100 || $svlen == ".") && scalar(split(/,/,$1)) > 1; print' file.vcf 
##fileformat=VCFv4.2
##source=combiSV-v2.2
##fileDate=Mon  May 8   11:32:53    2023
##contig=<ID=chrM,length=16571>
##contig=<ID=chr1,length=249250621> 
##INFO=<ID=END,Number=1,Type=Integer,Description="End   position    of  the variant described   in  this    record">
##INFO=<ID=SVCALLERS,Number=.,Type=String,Description="SV   callers that    support this    SV">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=DR,Number=1,Type=Integer,Description="#    High-quality    reference   reads">
##FORMAT=<ID=DV,Number=1,Type=Integer,Description="#    High-quality    variant reads">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample
1   10862   id.1    N   <INS>   .   PASS    SVTYPE=INS;SVLEN=101;END=10862;SVCALLERS=cutesv,SVIM    GT:DR:DV    1/1:0:26
1   90258   id.2    N   <INS>   .   PASS    SVTYPE=INS;SVLEN=118;END=90258;SVCALLERS=SVIM,NanoSV    GT:DR:DV    1/1:0:9
1   186241  id.5    N   <DEL>   .   PASS    SVTYPE=DEL;SVLEN=418;END=186662;SVCALLERS=SVIM,NanoSV   GT:DR:DV    1/1:2:12
1   526111  id.6    N   <DEL>   .   PASS    SVTYPE=DEL;SVLEN=624;END=526735;SVCALLERS=Sniffles,cutesv   GT:DR:DV0/1:8
2   91926078    id.3958 N   <BND>   .   PASS    SVTYPE=BND;SVLEN=.;END=;SVCALLERS=Sniffles,NanoSV   GT:DR:DV0/1:60:15
    • 5

  2. 这会忽略标题行，只处理数据行。

    start cmd:> awk -F '=|;| +' '$11<100 || $15 !~ "," { next; }; { print $0; }' input
1   10862   id.1    N   <INS>   .   PASS    SVTYPE=INS;SVLEN=101;END=10862;SVCALLERS=cutesv,SVIM    GT:DR:DV    1/1:0:26
1   90258   id.2    N   <INS>   .   PASS    SVTYPE=INS;SVLEN=118;END=90258;SVCALLERS=SVIM,NanoSV    GT:DR:DV    1/1:0:9
1   186241  id.5    N   <DEL>   .   PASS    SVTYPE=DEL;SVLEN=418;END=186662;SVCALLERS=SVIM,NanoSV   GT:DR:DV    1/1:2:12
1   526111  id.6    N   <DEL>   .   PASS    SVTYPE=DEL;SVLEN=624;END=526735;SVCALLERS=Sniffles,cutesv   GT:DR:DV    0/1:8

    更新 1

    对于调整后的问题

    start cmd:> awk -F '=|;| +' '/^#/ { print; next; }; $5 != "<BND>" && $11<100 || $15 !~ "," { next; }; $5 == "<BND>" && $15 !~ "," { next; }; { print $0; }' input
##fileformat=VCFv4.2
##source=combiSV-v2.2
##fileDate=Mon May  8 11:32:53 2023
##contig=<ID=chrM,length=16571>
##contig=<ID=chr1,length=249250621>
##INFO=<ID=END,Number=1,Type=Integer,Description="End position of the variant described in this record">
##INFO=<ID=SVCALLERS,Number=.,Type=String,Description="SV callers that support this SV">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=DR,Number=1,Type=Integer,Description="# High-quality reference reads">
##FORMAT=<ID=DV,Number=1,Type=Integer,Description="# High-quality variant reads">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  Sample
1   10862   id.1    N   <INS>   .   PASS    SVTYPE=INS;SVLEN=101;END=10862;SVCALLERS=cutesv,SVIM    GT:DR:DV    1/1:0:26
1   90258   id.2    N   <INS>   .   PASS    SVTYPE=INS;SVLEN=118;END=90258;SVCALLERS=SVIM,NanoSV    GT:DR:DV    1/1:0:9
1   186241  id.5    N   <DEL>   .   PASS    SVTYPE=DEL;SVLEN=418;END=186662;SVCALLERS=SVIM,NanoSV   GT:DR:DV    1/1:2:12
1   526111  id.6    N   <DEL>   .   PASS    SVTYPE=DEL;SVLEN=624;END=526735;SVCALLERS=Sniffles,cutesv   GT:DR:DV    0/1:8
2   91926078    id.3958 N   <BND>   .   PASS    SVTYPE=BND;SVLEN=.;END=;SVCALLERS=Sniffles,NanoSV   GT:DR:DV    0/1:60:15
    • 1

相关问题