Início / unix / Perguntas / 427831
Accepted
Age87
Asked: 2018-03-03 18:38:11 +0800 CST

extrair palavras do arquivo usando informações de um subconjunto do arquivo (várias etapas)

  • 772

Eu tenho dois arquivos grandes que se parecem com o seguinte:

Arquivo 1:

NW_006502347.1 316684
NW_006527876.1 351
NW_006502151.1 27628
NW_006526579.1 232
NW_006525259.1 132
NW_006501641.1 437014
NW_006525259.1 378
NW_006523082.1 215
NW_006522424.1 153
NW_006522101.1 815
NW_006521985.1 505
NW_006521985.1 527
NW_006521722.1 920
NW_006521525.1 73
NW_006521432.1 258
NW_006521302.1 938
NW_006521272.1 585
NW_006521272.1 745
NW_006521038.1 202
NW_006519846.1 1528
NW_006519837.1 10215

Arquivo 2:

NW_006502347.1  Gnomon  CDS   305319  305340   .  +  0  ID=cds43608  Parent=rna48098  Dbxref=GeneID:102908761,Genbank:XP_006997436.2  Name=XP_006997436.2  gbkey=CDS  gene=LOC102908761  partial=true  product=histone deacetylase 4-like  protein_id=XP_006997436.2
NW_006501037.1  Gnomon  gene  6936    115174   .  -  .  ID=gene0                      Dbxref=GeneID:102922816  Name=Efl1  gbkey=Gene  gene=Efl1  gene_biotype=protein_coding  partial=true  start_range=.,6936
NW_006501037.1  Gnomon  mRNA  6936    115174   .  -  .  ID=rna0      Parent=gene0     Dbxref=GeneID:102922816,Genbank:XM_006970114.2  Name=XM_006970114.2  gbkey=mRNA  gene=Efl1  model_evidence=Supporting evidence includes similarity to: 5 mRNAs%2C 7 Proteins%2C and 99%25 coverage of the annotated genomic feature by RNAseq alignments%2C including 16 samples with support for all annotated introns  partial=true  product=elongation factor like GTPase 1  start_range=.,6936  transcript_id=XM_006970114.2
NW_006501037.1  Gnomon  exon  115095  115174   .  -  .  ID=id1       Parent=rna0      Dbxref=GeneID:102922816,Genbank:XM_006970114.2  gbkey=mRNA  gene=Efl1  partial=true  product=elongation factor like GTPase 1  transcript_id=XM_006970114.2
NW_006501037.1  Gnomon  exon  114246  114355   .  -  .  ID=id2       Parent=rna0      Dbxref=GeneID:102922816,Genbank:XM_006970114.2  gbkey=mRNA  gene=Efl1  partial=true  product=elongation factor like GTPase 1  transcript_id=XM_006970114.2
NW_006502347.1  Gnomon  mRNA  272091  477077   .  +  .  ID=rna48098  Parent=gene26399 Dbxref=GeneID:102908761,Genbank:XM_006997374.2  Name=XM_006997374.2  end_range=477077,.  gbkey=mRNA  gene=LOC102908761  model_evidence=Supporting evidence includes similarity to: 1 mRNA%2C and 90%25 coverage of the annotated genomic feature by RNAseq alignments  partial=true  product=histone deacetylase 4-like  transcript_id=XM_006997374.2
NW_006501037.1  Gnomon  exon  92339   92472    .  -  .  ID=id5       Parent=rna0      Dbxref=GeneID:102922816,Genbank:XM_006970114.2  gbkey=mRNA  gene=Efl1  partial=true  product=elongation factor like GTPase 1  transcript_id=XM_006970114.2
NW_006501037.1  Gnomon  exon  90969   91106    .  -  .  ID=id6       Parent=rna0      Dbxref=GeneID:102922816,Genbank:XM_006970114.2  gbkey=mRNA  gene=Efl1  partial=true  product=elongation factor like GTPase 1  transcript_id=XM_006970114.2
NW_006501037.1  Gnomon  exon  89261   89475    .  -  .  ID=id7       Parent=rna0      Dbxref=GeneID:102922816,Genbank:XM_006970114.2  gbkey=mRNA  gene=Efl1  partial=true  product=elongation factor like GTPase 1  transcript_id=XM_006970114.2
NW_006502151.1  Gnomon  exon  26099   27657    .  -  .  ID=id586652  Parent=rna47002  Dbxref=GeneID:102918658,Genbank:XM_006996663.2  gbkey=mRNA  gene=Rftn1  product=raftlin%2C lipid raft linker 1  transcript_id=XM_006996663.2
NW_006501641.1  Gnomon  mRNA  393496  438556   .  +  .  ID=rna40001  Parent=gene21212 Dbxref=GeneID:102913870,Genbank:XM_015986269.1  Name=XM_015986269.1  gbkey=mRNA  gene=LOC102913870  model_evidence=Supporting evidence includes similarity to: 9 mRNAs%2C 5 Proteins%2C and 81%25 coverage of the annotated genomic feature by RNAseq alignments  product=transmembrane protein 189  transcript_id=XM_015986269.1
NW_006501053.1  Gnomon  exon  5104713 5104872  .  +  .  ID=id45206   Parent=rna3590   Dbxref=GeneID:102916769,Genbank:XR_001580019.1  gbkey=misc_RNA  gene=Mycbpap  product=MYCBP associated protein%2C transcript variant X4  transcript_id=XR_001580019.1
NW_006501053.1  Gnomon  exon  5104959 5105062  .  +  .  ID=id45207   Parent=rna3590   Dbxref=GeneID:102916769,Genbank:XR_001580019.1  gbkey=misc_RNA  gene=Mycbpap  product=MYCBP associated protein%2C transcript variant X4  transcript_id=XR_001580019.1
NW_006501053.1  Gnomon  exon  5105698 5105881  .  +  .  ID=id45208   Parent=rna3590   Dbxref=GeneID:102916769,Genbank:XR_001580019.1  gbkey=misc_RNA  gene=Mycbpap  product=MYCBP associated protein%2C transcript variant X4  transcript_id=XR_001580019.1
NW_006501053.1  Gnomon  exon  5106131 5106246  .  +  .  ID=id45209   Parent=rna3590   Dbxref=GeneID:102916769,Genbank:XR_001580019.1  gbkey=misc_RNA  gene=Mycbpap  product=MYCBP associated protein%2C transcript variant X4  transcript_id=XR_001580019.1

Desejo usar as informações do arquivo 1 para extrair a palavra seguinte "gene=" na coluna 13 ou 14 (por exemplo, "Efl1"). Mais especificamente, quero:

Etapa 1) Compare os rótulos da coluna 1 do arquivo 1 (por exemplo, NW_006527876.1) com os rótulos da coluna 1 do arquivo 2, extraia todas as linhas para as quais a coluna 1 (arquivo 2) corresponde à coluna 1 do arquivo 1.

Como você pode ver, os rótulos da coluna 1 (arquivo 2) se repetem, portanto, haverá várias correspondências para cada rótulo do arquivo 1.

As colunas 4 e 5 do arquivo 2 representam um intervalo, sendo a coluna 4 o início e a coluna 5 o final do intervalo. A coluna 2 do arquivo 1 representa os números entre esses intervalos.

Etapa 2) Das linhas isoladas da etapa 1, extraia as linhas para as quais o número na coluna 2 (arquivo 1) está entre o intervalo indicado pela coluna 4 e 5 do arquivo 2.

Isso está muito além do que eu sei, mas abaixo está uma ideia de como os comandos podem ser.

awk '{ print $1 }' file 1 |  
awk `$4 *(file2)* < $2 *(file1)*' | awk '$5 *(file2)* > $2 *(file1)*' > output.tsv

A saída deve ter linhas com rótulos exclusivos na coluna 1.

Etapa 3) Do arquivo output.tsv criado acima, gostaria de extrair a palavra que segue o sinal de igual em "gene=" na coluna 13 ou 14 (veja abaixo), de forma que acabe com um arquivo apenas com as palavras seguindo o sinal de igual.

Arquivo de saída final (com base neste exemplo):

LOC102908761
Rftn1
LOC102913870

RESPONDA:

while read -r id pos; do awk -v id="$id" -v pos="$pos" '$1 == id && pos > $4 && pos < $5 { if (gensub(/.*gene=([A-Za-z0-9]*).*/, "\\1", 1) !~ /\s/) print gensub(/.*gene=([A-Za-z0-9]*).*/, "\\1", 1); }' < file2.txt; done <file1.txt > gene_hits.txt

Em seguida, para se livrar de replicações (mesmo locus em vários intervalos), mantendo diferentes loci mapeados para o mesmo gene, faça:

perl -ne 'print if ++$k{$_}==1' A_gene_hits.txt > A_genes.txt
awk sed

1 respostas

  1. Best Answer

    Supondo que você tenha espaço em branco como delimitador:

    $ while read -r id pos; do awk -v id="$id" -v pos="$pos" '$1 == id && pos > $4 && pos < $5 { print gensub(/.*gene=([A-Za-z0-9]*).*/, "\\1", 1) }' <file2; done <file1
LOC102908761
Rftn1
LOC102913870

    Explicação

    • while read -r id pos; do FOO; done <file1: isso lê file1linha por linha e coloca o primeiro campo (por exemplo NW_006502347.1) na variável do shell $ide o segundo campo (por exemplo 316684) na variável do shell $pos. Em seguida, ele é executado FOOpara cada linha.
    • awk -v id="$id" -v pos="$pos" 'BAR' <file2: para cada linha de file1, executaremos um awkcomando que será executado BAR. Isso irá procurar file2as peças correspondentes. Precisamos informar a este awkscript as duas variáveis ​​"externas" do shell. isto é, a variável awk idrecebe o mesmo valor que a variável shell $id, e o mesmo para a variável awk pose a variável shell $pos.
    • $1 == id && pos > $4 && pos < $5: esta é a parte "condicional" do awkscript. Se essas condições forem atendidas, os seguintes comandos serão executados. Aqui, estamos verificando se o primeiro campo $1de file2é o mesmo da idlinha atual de file1, E se posestá entre o 4º $4e o 5º $5campos de file2.
    • { print gensub(/.*gene=([A-Za-z0-9]*).*/, "\\1", 1) }: se as condições acima forem atendidas, esse código será executado. Queremos fazer uma substituição com gensubprimeiro. Isso procura gene=seguido por uma string alfanumérica de qualquer comprimento ([A-Za-z0-9]*). Essa string alfanumérica é (capturada )pelos parênteses. Também "pesquisaremos" todos os caracteres antes e depois .*da string completa gene=([A-Za-z0-9]*). Portanto, isso "procura" a linha inteira, que é substituída pelo (primeiro e único) grupo de captura "\\1", ou seja, a sequência alfanumérica após gene=. O final 1significa substituir a primeira ocorrência, embora isso não faça muito sentido, pois presumo que haverá apenas uma correspondência gene=por linha.

    Versão delimitada por tabulação

    Prefiro usar arquivos delimitados por tabulação em geral, especialmente para o que presumo ser um arquivo GFF/GTF. Isso permite diferenciar os espaços, principalmente no 9º campo.

    while IFS=$'\t' read -r id pos; do awk -F'\t' -v id="$id" -v pos="$pos" '$1 == id && pos > $4 && pos < $5 { print gensub(/.*gene=([A-Za-z0-9]*).*/, "\\1", 1) }' <file2.tsv ; done <file1.tsv

    As modificações no script estão dividindo explicitamente linhas de shell em guias com IFS=$'\t', e awklinhas com -F'\t'.

    • 2

relate perguntas