☰

更新1

tag-icon • tag-icon text-processing command-line bioinformatics
更新1

我有一些.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包含 1 行的行。这是必须满足的两个标准，换句话说，我只想保留SVLEN> 100 且至少有 2 SVCALLERS) 的行。

此外，还有一些行 是 ，ALT并且BND该文件没有提供任何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

答案1

这是一个 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它自动在由（空格，默认情况下，如 awk）给出的字符上分割每个输入行-F并将其保存在数组中@F。因此，由于数组从 开始计数0，因此第一个字段是$F[0[，第二个字段$F[1]以此类推。接下来，从每个输入行中删除尾随换行符，并在每个调用中-l添加 a ，这意味着“逐行读取输入文件并将 给定的脚本应用于每一行。\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

答案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

相关内容