Usando Uniq -c com uma expressão regular ou contando o número de linhas removidas

Question

Eu tenho um arquivo delimitado por tabulações com informações sobre material genético. Algumas informações são cortadas em um arquivo de guia menor com algumas das colunas extraídas e o uniq é usado para garantir que não haja duplicatas. Uma contagem é armazenada, o que é importante mais tarde no pipeline. Eu quero alterar a base da função uniq de simplesmente basear-se na seqüência em um campo, para um regex de dentro de outro campo. O campo do qual quero extrair está neste formato:

14:50065421-50065521:12397472_t14_w100_x1

mas o bit após o segundo cólon muda dependendo da entrada do arquivo. Eu quero usar o uniq com base no primeiro semestre, então:

14:50065421-50065521

Eu testei o regex '((^ [0-9] {0,2} | x | y | MT): [0-9] {0,9} - [0-9] {0, 9} :) 'e funciona em uma pequena amostra dos dados. Eu encontrei algumas maneiras usando grep e um script perl que ambos trabalham para remover linhas com base no regex, mas nenhum deles fornece a contagem (é por isso que o uniq é muito mais ideal). Existe alguma maneira de usar um regex com uniq? Ou existe uma maneira melhor que também armazene a contagem daqueles removidos?

O código atual é este:

cat ${TAB_FILE} | \
    sed -e '1,2d' | \
    cut -f3,4 | \
    sort -k1 -u | \
    sort -k2 | \
    uniq -cf1| \
    sort -rn > t1

cat ${TAB_FILE} | \
    sed -e '1,2d' | \
    awk {'print $3"\t "$6"\t" $7"\t "$4'} | \
    sort -k4 > t2

awk 'FNR==NR {C[$2]=$1;next}FNR==1 \
    {print "Count Chromosome:Positions:QNAME Sequence Exon Transcript_ID"; next}$1 in C \
    {print C[$1], $1, $4, $3, $2}' t1 t2 > t3

cat t3 | awk '{print "//NODECLASS\t\"" $2"_"$1  "\"\t\"Exon " $4 "\"\t\"" $5 "\""}'

Na primeira etapa, é a coluna 1 do recorte que gostaria de basear minha regex em vez de usar a coluna 2. Qualquer ajuda seria muito apreciada, por favor, sinta-se à vontade para perguntar se você precisa de mim para esclarecer alguma coisa.

Exemplo de arquivo de guia:

queryHits   subjectHits readname    readSeq geneid  transcriptid    exonnumber  genename    biotype
350851  1   14:50065421-50065521:12397472_t14_w100_x1   CGCTGCCAGCTGCGCGCTCGGGGGAAAAGACGTTGCGCCCCCGCCGACTGCCGGTTTCCCGGGCGCGAGCCCGGATCCAGGTGGTCAGTCCCGGTACGCA    ENSG00000165501 ENST00000298288 1   LRR1    protein_coding
350851  5   14:50065421-50065521:12397472_t14_w100_x1   CGCTGCCAGCTGCGCGCTCGGGGGAAAAGACGTTGCGCCCCCGCCGACTGCCGGTTTCCCGGGCGCGAGCCCGGATCCAGGTGGTCAGTCCCGGTACGCA    ENSG00000165501 ENST00000318317 1   LRR1    protein_coding
350851  8   14:50065421-50065521:12397472_t14_w100_x1   CGCTGCCAGCTGCGCGCTCGGGGGAAAAGACGTTGCGCCCCCGCCGACTGCCGGTTTCCCGGGCGCGAGCCCGGATCCAGGTGGTCAGTCCCGGTACGCA    ENSG00000165501 ENST00000554869 1   LRR1    protein_coding
350852  1   14:50065461-50065561:12655987_t14_w100_x1   CCGCCGACTGCCGGTTTCCCGGGCGCGAGCCCGGATCCAGGTGGTCAGTCCCGGTACGCAACCACGGCGAGAACCCGGCCCTGCTAAGGGAGAAGGGAAG    ENSG00000165501 ENST00000298288 1   LRR1    protein_coding
350852  5   14:50065461-50065561:12655987_t14_w100_x1   CCGCCGACTGCCGGTTTCCCGGGCGCGAGCCCGGATCCAGGTGGTCAGTCCCGGTACGCAACCACGGCGAGAACCCGGCCCTGCTAAGGGAGAAGGGAAG    ENSG00000165501 ENST00000318317 1   LRR1    protein_coding
350852  8   14:50065461-50065561:12655987_t14_w100_x1   CCGCCGACTGCCGGTTTCCCGGGCGCGAGCCCGGATCCAGGTGGTCAGTCCCGGTACGCAACCACGGCGAGAACCCGGCCCTGCTAAGGGAGAAGGGAAG    ENSG00000165501 ENST00000554869 1   LRR1    protein_coding
350853  1   14:50065471-50065571:22679947_t13_w100_x1   CCGGTTTCCCGGGCGCGAGCCCGGATCCAGGTGGTCAGTCCCGGTACGCAACCACGGCGAGAACCCGGCCCTGCTAAGGGAGAAGGGAAGCCGTTTCCCG    ENSG00000165501 ENST00000298288 1   LRR1    protein_coding
350853  5   14:50065471-50065571:22679947_t13_w100_x1   CCGGTTTCCCGGGCGCGAGCCCGGATCCAGGTGGTCAGTCCCGGTACGCAACCACGGCGAGAACCCGGCCCTGCTAAGGGAGAAGGGAAGCCGTTTCCCG    ENSG00000165501 ENST00000318317 1   LRR1    protein_coding
350853  8   14:50065471-50065571:22679947_t13_w100_x1   CCGGTTTCCCGGGCGCGAGCCCGGATCCAGGTGGTCAGTCCCGGTACGCAACCACGGCGAGAACCCGGCCCTGCTAAGGGAGAAGGGAAGCCGTTTCCCG    ENSG00000165501 ENST00000554869 1   LRR1    protein_coding

Exemplo de arquivo de saída em que o uniq não foi usado, desejo remover as duplicatas, como aquelas nas linhas 4, 5 e 6 ou 8, 9 e 10:

//NODECLASS "Chromosome:Positions:QNAME_Count"  "Exon Exon" "Transcript_ID"

//NODECLASS "14:50067283-50067383:20149917_t14_w100_x1_1"   "Exon 1"    "ENST00000557531"
//NODECLASS "14:50067284-50067366:14257122_t14_w100_x1_2"   "Exon 1"    "ENST00000557531"
//NODECLASS "14:50067285-50067385:2072777_t12_w100_x1_1"    "Exon 1"    "ENST00000557531"
//NODECLASS "14:50074262-50074362:4355312_t12_w100_x1_1"    "Exon 3"    "ENST00000298288"
//NODECLASS "14:50074262-50074362:4355312_t12_w100_x1_1"    "Exon 4"    "ENST00000540712"
//NODECLASS "14:50074262-50074362:4355312_t12_w100_x1_1"    "Exon 4"    "ENST00000554869"
//NODECLASS "14:50067286-50067386:15839225_t12_w100_x1_3"   "Exon 1"    "ENST00000557531"
//NODECLASS "14:50074263-50074363:8914169_t11_w100_x1_1"    "Exon 3"    "ENST00000298288"
//NODECLASS "14:50074263-50074363:8914169_t11_w100_x1_1"    "Exon 4"    "ENST00000540712"
//NODECLASS "14:50074263-50074363:8914169_t11_w100_x1_1"    "Exon 4"    "ENST00000554869"
//NODECLASS "14:50067287-50067387:5439923_t13_w100_x1_1"    "Exon 1"    "ENST00000557531"
//NODECLASS "14:50067287-50067387:14106336_t12_w100_x1_3"   "Exon 1"    "ENST00000557531"
//NODECLASS "14:50074404-50074504:15492363_t11_w100_x1_1"   "Exon 3"    "ENST00000298288"
//NODECLASS "14:50074404-50074504:15492363_t11_w100_x1_1"   "Exon 4"    "ENST00000540712"
//NODECLASS "14:50074404-50074504:15492363_t11_w100_x1_1"   "Exon 4"    "ENST00000554869"
//NODECLASS "14:50074135-50074235:11346262_t11_w100_x1_2"   "Exon 3"    "ENST00000298288"

    

Answer 1

Eu não tenho certeza do que você quer exatamente, já que você parece excessivamente obcecado em seu tipo / coisas uniq.

No entanto, se tudo que você queria fosse a remoção de linhas repetitivas, por exemplo, 4,5,6 e 8,9,10, você poderia fazer isso no último arquivo que está mostrando, o qual tem // linhas NODECLASS:

perl -F\" -lane '
   print,next if $. < 3;
   print if ! $h{($F[1] =~ /:(.*?):/)[0]}++;
' NODE_CLASS_file

onde as duas primeiras linhas são ignoradas. Enquanto para o resto, nós olhamos para o 2º campo o número entre: /:(.*?):/ irá lhe fornecer o número então você precisa entrar no contexto escalar colocando-o em (...)[0] e passá-lo como a chave para o hash %código%. Imprima a linha atual somente se essa chave já tiver sido vista.

    

Answer 2

Consegui encontrar uma solução bruta em que adiciono as posições de início e fim no arquivo da guia como uma coluna separada. Então, ao criar o arquivo Nodeclass, também extraio essa coluna e uso sort | uniq -c baseado nisso em vez da coluna de seqüência. Parece estar fazendo o que eu queria, mas é um pouco mais lento porque o script que cria o arquivo de guia precisa produzir uma coluna extra inteira! Obrigado :)