Блог

Я пытаюсь обработать файл GTF / GFF, который я загружаю из ensemble. Усеченная версия файла выглядит следующим образом:

 1   ensembl gene    5273    10061   .   -   .   gene_id ENSGALG00000054818; gene_version 1; gene_source ensembl; gene_biotype protein_coding;
1   ensembl transcript  5273    10061   .   -   .   gene_id ENSGALG00000054818; gene_version 1; transcript_id ENSGALT00000098984; transcript_version 1; gene_source ensembl; gene_biotype protein_coding; transcript_source ensembl; transcript_biotype protein_coding;
1   ensembl gene    58427   58617   .       .   gene_id ENSGALG00000047594; gene_version 1; gene_name RF00004; gene_source ensembl; gene_biotype snRNA;
1   ensembl transcript  58427   58617   .       .   gene_id ENSGALG00000047594; gene_version 1; transcript_id ENSGALT00000094382; transcript_version 1; gene_name RF00004; gene_source ensembl; gene_biotype snRNA; transcript_name RF00004-201; transcript_source ensembl; transcript_biotype snRNA;
1   ensembl exon    58427   58617   .       .   gene_id ENSGALG00000047594; gene_version 1; transcript_id ENSGALT00000094382; transcript_version 1; exon_number 1; gene_name RF00004; gene_source ensembl; gene_biotype snRNA; transcript_name RF00004-201; transcript_source ensembl; transcript_biotype snRNA; exon_id ENSGALE00000460125; exon_version 1;
1   ensembl gene    63264   63454   .       .   gene_id ENSGALG00000049206; gene_version 1; gene_name RF00004; gene_source ensembl; gene_biotype snRNA;
1   ensembl transcript  63264   63454   .       .   gene_id ENSGALG00000049206; gene_version 1; transcript_id ENSGALT00000092780; transcript_version 1; gene_name RF00004; gene_source ensembl; gene_biotype snRNA; transcript_name RF00004-201; transcript_source ensembl; transcript_biotype snRNA;
1   ensembl exon    63264   63454   .       .   gene_id ENSGALG00000049206; gene_version 1; transcript_id ENSGALT00000092780; transcript_version 1; exon_number 1; gene_name RF00004; gene_source ensembl; gene_biotype snRNA; transcript_name RF00004-201; transcript_source ensembl; transcript_biotype snRNA; exon_id ENSGALE00000501941; exon_version 1;
  

(Девять столбцов, разделенных табуляцией.)

В некоторых строках отсутствуют атрибуты, такие как gene_name , transcript_id или transcript_name .

• Если gene_name отсутствует, я хотел заменить его на gene_id ,
• и если transcript_name отсутствует, я хотел заменить его на transcript_id (в случае отсутствия transcript_id оно заменяется на gene_id ).

Однако информация для transcript_id или, лучше сказать, положение этой информации неизвестно. Как бы я искал атрибут и, в случае его отсутствия, заменил его значением transcript_id с неизвестной информацией о местоположении

Я добился замены отсутствующего значения для gene_name на значение для gene_id следующим образом:

 awk '{if (!/gene_name/) print $0, "gene_name " $10; else print $0}' input.gtf > output.gtf
  

Это сработало довольно хорошо, но только потому, что в данном конкретном случае я знал позицию значения, которое я использовал в качестве замены. Я не мог понять, как я добьюсь этого, когда позиция совпадения неизвестна.

Я использовал следующий код для получения неизвестной информации о местоположении, но не смог интегрировать проверку на несоответствие, как в первом примере выше:

 awk '{for (i=1; i<=NF;   i) { if ($i ~ "transcript_name") print$0,"transcript_name ", $(i 1) } }' input.gtf > output.gtf
  

Условие состоит в том, что только если transcript_name оно еще не присутствует в строке, его следует заменить значением для transcript_id .

Я действительно был бы признателен за помощь в этом!

Использование скрипта awk;

script.awk:

 #!/usr/bin/awk -f
BEGIN {
  FS=OFS="t"
}
{
  gsub(/; *$/, "", $9)        # trim trailing `;'
  split($9, pairs, / *; */)   # split attributes into pairs
  for (i in pairs) {
    split(pairs[i], kv, / */) # split pair into key and value
    attr[kv[1]] = kv[2]       # add it to `attr'
  }
  # fill missing fields
  if (!("gene_name" in attr))
    attr["gene_name"] = attr["gene_id"]
  if (!("transcript_id" in attr))
    attr["transcript_id"] = attr["gene_id"]
  if (!("transcript_name" in attr))
    attr["transcript_name"] = attr["transcript_id"];
  # recreate the attributes field
  attr_all = sep = ""
  for (k in attr) {
    attr_all = attr_all sep k " " attr[k]
    sep = "; "
  }
  # update the record with new attributes
  $9 = attr_all 
}
1 # print record
  

Пример использования:

 awk -f script.awk inputfile
  

Онлайн-демонстрация.