extrai a palavra do arquivo usando informações de um subconjunto do arquivo (várias etapas)

Question

Eu tenho dois arquivos grandes que parecem com os seguintes:

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

Eu quero usar 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) aos 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, com a coluna 4 sendo o início e a coluna 5 o final do intervalo. A coluna 2 do arquivo 1 representa 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) fica entre o intervalo denotado pelas colunas 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 parecer.

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) No 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 tal forma que acabo com um arquivo apenas com as palavras após o sinal de igual.

Arquivo final de saída (baseado neste exemplo):

LOC102908761
Rftn1
LOC102913870

RESPOSTA:

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

Próximo para se livrar de réplicas (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

    

Answer 1

Supondo que você tenha espaço em branco como um 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 : lê file1 linha por linha e coloca o primeiro campo (por exemplo, NW_006502347.1 ) na variável de casca $id e o segundo campo (por exemplo, 316684 ) na variável de casca $pos . Em seguida, ele executa FOO para cada linha.
• awk -v id="$id" -v pos="$pos" 'BAR' <file2 : para cada linha de file1 , executaremos um comando awk que executará BAR . Isto irá procurar file2 para as partes correspondentes. Precisamos dizer a este script awk as duas variáveis "externas" do shell. ou seja, a variável awk id recebe o mesmo valor que a variável de shell $id e a mesma para a variável awk pos e a variável de shell $pos .
• $1 == id && pos > $4 && pos < $5 : esta é a parte "condicional" do script awk . Se essas condições forem atendidas, os seguintes comandos serão executados. Aqui, estamos verificando se o primeiro campo $1 de file2 é igual ao id da linha atual de file1 , E se pos está entre o 4º $4 e 5 $5 campos 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 gensub primeiro. Isso procura por gene= seguido por uma string alfanumérica de qualquer tamanho ([A-Za-z0-9]*) . Esta cadeia alfanumérica é ( capturada ) pelos parênteses. Também "procuraremos" por todos os caracteres antes e depois de .* da string completa gene=([A-Za-z0-9]*) . Portanto, isso "procura" pela 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 1 significa substituir a primeira ocorrência, embora isso não faça muito sentido, pois suponho que haverá apenas uma correspondência de gene= por linha.

Versão delimitada por tabulação

Eu prefiro usar arquivos delimitados por tabulação em geral, especialmente para o que eu suponho ser um arquivo GFF / GTF. Isso permite diferenciar espaços, especialmente 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 awk linhas com -F'\t' .