Início / unix / Perguntas / 746286
Accepted
Anna1364
Asked: 2023-05-19 01:21:24 +0800 CST

filtrar linhas com base em alguns critérios

  • 772

Eu tenho alguns arquivos vcf e quero filtrar algumas variantes. Esta é apenas uma pequena parte do meu vcf: existem algumas linhas de cabeçalho no início do arquivo (começando com ##) e depois variantes (uma linha por variante).

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

Mantendo as linhas de cabeçalho, quero remover as linhas com SVLEN < 100 e aquelas com apenas um SVCALLERS incluído (são dois critérios que ambos devem atender, em outras palavras, quero manter apenas as linhas com SVLEN > 100 e pelo menos dois SVCALLERS ). Além disso, existem algumas linhas que ALT é BND e vcf não fornece nenhum SVLEN para esse tipo de variante, se a linha contiver BND, quero apenas mantê-la se for suportada por dois chamadores. Exemplos: desejo descartar esta variante porque o SVLEN é menor que 100 e apenas um SVCALLERS o detectou

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

Ou esta linha também, embora haja dois chamadores, mas SVLEN é menor que 100:

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

Existe uma maneira fácil de fazer isso? Obrigado

Meu arquivo final deve ficar assim:

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

  1. Best Answer

    Aqui está uma maneira 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

    Explicação

    • perl -F'\t' -lane: o -afaz o perl funcionar como awkse ele dividisse automaticamente cada linha de entrada no caractere fornecido por -F(espaço em branco, por padrão, como awk) e o salvasse no array @F. Portanto, como os arrays começam a contar em 0, o primeiro campo é $F[0[, o segundo $F[1]e assim por diante. Em seguida, -lremove as novas linhas à direita de cada linha de entrada e também adiciona um \na cada printchamada, o -nque significa "ler o arquivo de entrada linha por linha e aplicar o script fornecido por -ea cada linha.
    • if(/^#/){ print; next }; : se esta for uma linha de cabeçalho, basta imprimi-la e passar para a próxima linha.
    • $F[7] =~ /\bSVLEN=(\d+)/; $svlen=$1;: corresponde à sequência de números mais longa depois SVLEN=no 8º campo e salve-a como $svlen. Este será o comprimento. Isso \bgarante que correspondamos apenas nos limites das palavras, portanto, isso não falhará se você tiver algo parecido NOTSVLEN=em seu arquivo.
    • $F[7] =~ /\bSVCALLERS=([^;]+)/; @callers=split(/,/,$1);: agora procure a string no 8º campo SVCALLERS=, pegue o trecho mais longo de não ;caracteres depois dela e, em seguida, divida-o ,na matriz @callers. Isso agora tem a lista de chamadores CNV usados ​​para este CNV.
    • print if $svlen > 100 && scalar(@callers) > 1: agora imprimimos a linha se o comprimento for maior que 100 e o número de chamadores ( scalar(@array)fornece o número de elementos em uma matriz) for maior que 1.

    E aqui está a mesma coisa básica de uma forma mais concisa e menos clara se você preferir seus comandos levemente golfados:

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

    Se você também deseja manter as linhas sem SVLEN desde que tenham pelo menos dois chamadores, use isto:

    $ 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. Isso ignora as linhas de cabeçalho, apenas lida com as linhas de dados.

    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

    Atualização 1

    Para a pergunta ajustada

    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

relate perguntas