Блог

Я работаю над скриптом, который извлекает строки из файла GTF с помощью gawk, и он добавляет дополнительный пробел между последовательными вызовами «printf». Для тех, кто не знаком с файлами GTF, они представляют собой обычный геномный формат, состоящий из 9 полей, разделенных табуляцией, где в 9-м поле хранится список пар атрибутов ключ — значение, разделенных комбинированной точкой с запятой и пробелом «; «. Цель состоит в том, чтобы извлечь строки с определенным «gene_name», переданные в виде одного столбца текста во входном файле 1.

Все в сценарии обработки работает так, как ожидалось, ЗА исключением того, что каким-то образом вводится дополнительный пробел между последней итерацией printf во внутреннем цикле for и инструкцией printf, которая вставляет символ «новой строки».

Пример входного файла 1:

 (base) [user@host MouseEnsembl100]$ head gene_names.txt
Cryaa
Cryab
Crygc
  

Пример входного файла GTF (файл 2):

 (base) [user@host MouseEnsembl100]$ head example.gtf
17  ensembl_havana  gene    31677807    31681733    .       .   gene_id "ENSMUSG00000024041"; gene_version "10"; gene_name "Cryaa";
17  havana  transcript  31677807    31681733    .       .   gene_id "ENSMUSG00000024041"; gene_version "10"; transcript_id "ENSMUST00000228716"; gene_name "Cryaa";
17  havana  exon    31677807    31678189    .       .   gene_id "ENSMUSG00000024041"; gene_version "10"; transcript_id "ENSMUST00000228716"; gene_name "Cryaa";
17  havana  CDS 31678001    31678189    .       0   gene_id "ENSMUSG00000024041"; gene_version "10"; transcript_id "ENSMUST00000228716"; gene_name "Cryaa";
17  havana  start_codon 31678001    31678003    .       0   gene_id "ENSMUSG00000024041"; gene_version "10"; transcript_id "ENSMUST00000228716"; gene_name "Cryaa";
17  havana  exon    31679559    31679681    .       .   gene_id "ENSMUSG00000024041"; gene_version "10"; transcript_id "ENSMUST00000228716"; gene_name "Cryaa";
  

Сценарий обработки:

 #!/bin/bash
#SBATCH --job-name=make_rseqc_bed
#SBATCH --mem=32000
#SBATCH --ntasks=4

TARGETS=/work/abf/MouseEnsembl100/gene_names.txt
TGTLABL=gene_name
GTFPATH=/work/abf/MouseEnsembl100/example.gtf

if [ ! -z $TGTLABL ] amp;amp; [ ! -z $TARGETS ]
then
    gawk -v lbl=${TGTLABL}
     -v FS="t| |;"
     -v OFS=''
     -v ORS=''
     '(NR == FNR) {tgt[$1]; next}
      (NR != FNR) {gsub("; ",";")}
          (NR != FNR)
              {
                 for(i=0; i<=NF; i  ){                     
                     if($i == lbl){
                         gsub("42","",$(i 1))
                         if($(i 1) in tgt){
                             $(i 1)="42"$(i 1)"42"
                             for(j=1; j<=NF;j  ){
                                if(j < 9) {
                                    printf($j"t")
                                }
                                else if( (j % 2) == 1){
                                    printf($j" ")
                                }
                                else if( (j % 2) == 0 amp;amp; (j 1) < NF){
                                    printf($j"; ")
                                }
                                else if((j 1) == NF){
                                    printf($j";LAST_FIELD")
                                }
                             }
                             printf("%sn","NEXT LINE")
                         }
                  }   
                 } 
              }' $TARGETS $GTFPATH >> extracted_targets.gtf
fi
  

Пример вывода:

 17  ensembl_havana  gene    31677807    31681733    .       .   gene_id "ENSMUSG00000024041"; gene_version "10"; gene_name "Cryaa"; NEXT LINE
17  havana  transcript  31677807    31681733    .       .   gene_id "ENSMUSG00000024041"; gene_version "10"; transcript_id "ENSMUST00000228716"; gene_name "Cryaa"; NEXT LINE
17  havana  exon    31677807    31678189    .       .   gene_id "ENSMUSG00000024041"; gene_version "10"; transcript_id "ENSMUST00000228716"; gene_name "Cryaa"; NEXT LINE
17  havana  CDS 31678001    31678189    .       0   gene_id "ENSMUSG00000024041"; gene_version "10"; transcript_id "ENSMUST00000228716"; gene_name "Cryaa"; NEXT LINE
17  havana  start_codon 31678001    31678003    .       0   gene_id "ENSMUSG00000024041"; gene_version "10"; transcript_id "ENSMUST00000228716"; gene_name "Cryaa"; NEXT LINE
17  havana  exon    31679559    31679681    .       .   gene_id "ENSMUSG00000024041"; gene_version "10"; transcript_id "ENSMUST00000228716"; gene_name "Cryaa"; NEXT LINE
  

После нескольких часов попыток разных действий я наконец понял, что с помощью установленного разделителя полей FS="t| |;" Gawk генерирует пустое поле для терминала ; (В моем файле GTF каждая строка заканчивается точкой с запятой).

Когда я изменился:

 for(j=1; j<=NF;j  ){
  

Для:

   for(j=1; j<NF;j  ){
  

Проблема была решена. Надеюсь, кто-нибудь сочтет этот пост полезным.