Блог

У меня file1.gz это выглядит так

 NC_044998.1     3756    0       AA
NC_044998.1     3757    0       AA
NC_044998.1     3758    0       AT
NC_044998.1     3759    0       GG
NC_044998.1     3760    0       GG
NC_044998.1     3761    0       GC
NC_044998.1     3763    0       TT
NC_044998.1     3764    0       GG
  

Мне нужно аннотировать file1.gz в соответствии с этим col4 . If col4 — это 2x та же буква (что и в GG ) col5 homozygous , что и «else print « heterozygous «. Вот так:

 NC_044998.1     3756    0       A    A       homozygous
NC_044998.1     3757    0       A    A       homozygous
NC_044998.1     3758    0       A    T       heterozygous
NC_044998.1     3759    0       G    G       homozygous
NC_044998.1     3760    0       G    G       homozygous
NC_044998.1     3761    0       G    C       heterozygous
NC_044998.1     3763    0       T    T       homozygous
NC_044998.1     3764    0       G    G       homozygous
  

Решил это, выполнив (за исключением того, что разбивается col4 на 2).

 gunzip -c file1.gz | awk -v OFS="t" 'sub(/./,"amp;t",$4)1;' | awk -F $'t' 'BEGIN {OFS = FS} {if ($4==$5) print $0,"homozygous"; else print $0,"heterozygous";}' > file1.tmp
  

Затем мне нужно сравнить этот результат с file2.gz приведенным ниже

 chromo  position        major   minor   ref     anc     knownEM nInd
NC_044998.1     3756    A       C       A       A       0.000000        1
NC_044998.1     3757    A       C       T       T       0.000000        1
NC_044998.1     3758    A       C       A       A       0.000000        1
NC_044998.1     3759    G       A       G       G       0.000000        1
NC_044998.1     3760    G       A       G       G       0.000000        1
NC_044998.1     3761    G       A       G       G       0.000000        1
NC_044998.1     3762    A       C       A       A       0.000000        1
NC_044998.1     3764    G       A       C       C       0.000000        1
  

Добавление аннотации, чтобы, если col1, col2, col4 file1.gz == col1, col2, col5 file2.gz я добавлю « ref » в свою аннотацию и col1, col2, col5 file1.gz == col1, col2, col5 file2.gz добавлю « ref «. В противном случае добавьте « dev » в обоих случаях.

таким образом, конечный результат будет:

 NC_044998.1     3756    0       A    A       homozygous    ref    ref
NC_044998.1     3757    0       A    A       homozygous    dev    dev
NC_044998.1     3758    0       A    T       heterozygous  ref    dev
NC_044998.1     3759    0       G    G       homozygous    ref    ref
NC_044998.1     3760    0       G    G       homozygous    ref    ref
NC_044998.1     3761    0       G    C       heterozygous  ref    dev
NC_044998.1     3763    0       T    T       homozygous 
NC_044998.1     3764    0       G    G       homozygous    dev    dev
  

Я знаю, как это сделать отдельно

 gunzip -c  file2.gz | awk -F $'t' 'BEGIN {OFS = FS} NR==FNR{c[$1,$2,$5];next} {if (($1,$2,$4) in c ) print $0,"ref"; else print $0,"dev";}' - file1.tmp > file2.tmp

gunzip -c  file2.gz |awk -F $'t' 'BEGIN {OFS = FS} NR==FNR{c[$1,$2,$5];next} {if (($1,$2,$5) in c ) print $0,"ref"; else print $0,"dev";}' - file2.tmp > output
  

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

Вот awk для объединения последних 2 шагов в один:

 awk 'BEGIN {
   FS=OFS="t"
}
FNR==NR {
   a[$1,$2,$5]
   next
}
FNR==1 {
   next
}
{
   print $0, (($1,$2,$4) in a ? "ref":"dev"), (($1,$2,$5) in a ? "ref":"dev")
} file2 file1'
  

кстати, ваш первый сценарий awk также может быть выполнен за один шаг, как это:

 awk 'BEGIN {
   FS=OFS="t"
}
length($NF)==2 {
   p=substr($NF,1,1)
   s=substr($NF,2,1)
   $NF = p OFS s OFS (p==s ? "homozygous" : "heterozygous")
} 1' file