☰

アップデート1

tag-icon • tag-icon text-processing command-line bioinformatics
アップデート1

いくつかのファイルがあり、.vcfいくつかのバリエーションをフィルタリングしたいと思います。これは私のファイルのほんの一部です。ファイルの先頭（##で始まる）にはいくつかのヘッダー行があり、それにはバリアント（各バリアントごとに1行）があります。

##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行未満の行とSVCALLERS1行だけを含む行を削除したいと思います。これは満たさなければならない2つの基準です。つまり、SVLEN少なくとも2つ以上の行> 100だけ維持したいのですSVCALLERS。

また、いくつかの行があり、ファイルはそのようなバリエーションを提供しないので、行にが含まれているALT場合は両方の呼び出し元がサポートしている場合にのみ維持したいと思います。BNDSVLENBND

SVLEN例: このバリアントが 100 個未満で 1 つしか検出SVCALLERSされないため、削除したいと思います。

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

または、次の行は呼び出し側が2人ですが、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：これは、Perlが（デフォルトではawkなどの空白）、指定された文字の各入力行を自動的に分割し、それを配列に保存するように-a機能します。したがって、配列は最初から計算されるため、最初のフィールドはで、2番目のフィールドはこのようになります。次に、各入力行から末尾の改行を削除し、各呼び出しにaを追加します。これは、「入力ファイルを1行ずつ読み込み、与えられたスクリプトを各行に適用することを意味します。awk-F@F0$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=最長の非文字セグメントを取得して配列に;分割します。これで、このCNVのCNV送信者のリストがあります。,@callers
  • 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なしで回線を維持するには、発信者が2人以上の場合は、次のコマンドを使用します。

$ 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

おすすめ記事