根据某些标准过滤行

Question

根据某些标准过滤行

Ann*_*364 1 command-line text-processing bioinformatics

我有一些.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

Run Code Online (Sandbox Code Playgroud)

在保留标题行的同时，我想删除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

Run Code Online (Sandbox Code Playgroud)

或者这一行，虽然有两个调用者，但 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

Run Code Online (Sandbox Code Playgroud)

有简单的方法吗？谢谢

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

##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

Run Code Online (Sandbox Code Playgroud)

Answer 1

ter*_*don 7

这是一个 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

Run Code Online (Sandbox Code Playgroud)

解释

perl -F'\t' -lane：这-a使得 perl 的工作方式有点像awk它自动在由（空格，默认情况下，如 awk）给出的字符上分割每个输入行-F并将其保存在数组中@F。因此，由于数组从开始计数0，因此第一个字段是$F[0[，第二个字段$F[1]以此类推。接下来，-l从每个输入行中删除尾随换行符，并\n在每个print调用中添加 a ，这-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

Run Code Online (Sandbox Code Playgroud)

如果您还想保留没有 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

Run Code Online (Sandbox Code Playgroud)

归档时间：	2 年，5 月前
查看次数：	82 次
最近记录：	2 年，2 月前