いくつかのファイルがあり、.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行未満の行とSVCALLERS
1行だけを含む行を削除したいと思います。これは満たさなければならない2つの基準です。つまり、SVLEN
少なくとも2つ以上の行> 100だけ維持したいのですSVCALLERS
。
また、いくつかの行があり、ファイルはそのようなバリエーションを提供しないので、行にが含まれているALT
場合は両方の呼び出し元がサポートしている場合にのみ維持したいと思います。BND
SVLEN
BND
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
@F
0
$F[0[
$F[1]
-l
\n
print
-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
答え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