Создать последовательность обратного комплемента на основе AWK

Question

Уважаемые пользователи stackoverflow,

У меня есть данные TAB sep, например:

head -4 input.tsv

seq A      C change
seq T      A ok
seq C      C change
seq AC   CCT change

И мне нужно создать функцию обратного дополнения в awk, которая делает что-то вроде этого

head -4 output.tsv

seq T      G change
seq T      A ok
seq G      G change
seq GT   AGG change

Итак, если 4-й столбец отмечен как «изменение», мне нужно создать обратную последовательность дополнений.

СОВЕТ - то же самое, например tr в bash - Bash один вкладыш для этой задачи:

echo "ACCGA" | rev | tr "ATGC" "TACG"

Я пробовал что-то вроде этого

awk 'BEGIN {c["A"] = "T"; c["C"] = "G"; c["G"] = "C"; c["T"] = "A" }{OFS="\t"}
function revcomp(   i, o) {
    o = ""
    for(i = length; i > 0; i--)
        o = o c[substr($0, i, 1)]
    return(o)
}
{

if($4 == "change"){$2 = revcom(); $3 = revcom()} print $0; else print $0}' input

Значение обратной биологической последовательности:

A => T
C => G
G => C
T => A

и обратное дополнение означает:

ACCATG => CATGGT

Отредактировано: Также любой просто для образования может поделиться этим решением в python.

Answer 1

С небольшими переделками вашей попытки вы можете сделать что-то вроде ниже.

function revcomp(arg) {
    o = ""
    for(i = length(arg); i > 0; i--)
        o = o c[substr(arg, i, 1)]
    return(o)
}

BEGIN {c["A"] = "T"; c["C"] = "G"; c["G"] = "C"; c["T"] = "A" ; OFS="\t"}

{
    if($4 == "change") {
        $2 = revcomp($2); 
        $3 = revcomp($3)
    } 
}1

Ключевым моментом здесь было использование функции revcomp для принятия аргументов в качестве значений столбца и работы с ними путем итерации с конца. Раньше вы делали по всей строке $0, т.е. substr($0, i, 1), что вызывало бы множество необычных поисков в массиве c.

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

Если вы намереваетесь использовать вышеуказанное в части более крупного скрипта, я бы рекомендовал поместить весь код, как указано выше, в файл сценария, установите интерпретатор she-bang на #!/usr/bin/awk -f и запустите сценарий как awk -f script.awk input.tsv

Грубая версия bash, реализованная в awk, как показано ниже. Обратите внимание, что это не чисто и не рекомендуется 1021 * подход. Подробнее см. AllAboutGetline

Как и раньше, вызовите функцию как $2 = revcomp_bash($2) и $3 = revcomp_bash($3)

function revcomp_bash(arg) {
    o = ""
    cmd = "printf \"%s\" " arg "| rev | tr \"ATGC\" \"TACG\""
    while ( ( cmd | getline o ) > 0 ) { 

    }
    close(cmd);
    return(o)
}

Весь ваш код говорит о GNU awk -ism, поэтому не заботился о преобразовании его в совместимый с POSIX. Вы можете использовать split() с пустым де-ограничителем вместо length(), но в спецификации POSIX с радостью сказано, что «эффект пустой строки как значения fs не определен».

Answer 2

Не могли бы вы попробовать следующие, написанные и протестированные с показанными образцами (в GNU awk).

awk '
BEGIN{
  label["A"]="T"
  label["C"]="G"
  label["G"]="C"
  label["T"]="A"
}
function cVal(field){
  delete array
  num=split($field,array,"")
  for(k=1;k<=num;k++){
   if(array[k] in label){
     val=label[array[k]] val
   }
  }
  $field=val
  val=""
}
$NF=="change"{
  for(i=2;i<=(NF-1);i++){
    cVal(i)
  }
}
1
'  Input_file | column -t

Пояснение: Добавление подробного объяснения для код выше.

awk '                                ##Starting awk program from here.
BEGIN{                               ##Starting BEGIN section of this code here.
  label["A"]="T"                     ##Creating array label with index A and value T.
  label["C"]="G"                     ##Creating array label with index C and value G.
  label["G"]="C"                     ##Creating array label with index G and value C.
  label["T"]="A"                     ##Creating array label with index T and value A.
}
function cVal(field){                ##Creating function named cVal here with passing value field in it.
  delete array                       ##Deleting array here.
  num=split($field,array,"")         ##Splitting current field value passed to it and creating array.
  for(k=1;k<=num;k++){               ##Running for loop fromk=1 to till value of num here.
   if(array[k] in label){            ##Checking condition if array value with index k is present in label array then do following.
     val=label[array[k]] val         ##Creating val which has label value with index array with index k and keep concatenating its value to it.
   }
  }
  $field=val                         ##Setting current field value to val here.
  val=""                             ##Nullifying val here.
}
$NF=="change"{                       ##Checking condition if last field is change then do following.
  for(i=2;i<=(NF-1);i++){            ##Running for loop from 2nd field to 2nd last field.
    cVal(i)                          ##Calling function with passing current field number to it.
  }
}
1                                    ##1 will print current line here.
' Input_file | column -t             ##Mentioning Input_file name here.

Answer 3

Вроде неэффективно для этого конкретного приложения, поскольку он создает массив сопоставления при каждом вызове tr() и делает то же самое l oop в tr(), а затем снова в rev(), но решил, что покажу, как писать автономные функции tr() и rev(), и в любом случае, вероятно, будет достаточно быстро для ваших нужд:

$ cat tst.awk
BEGIN { FS=OFS="\t" }
$4 == "change" {
    for ( i=2; i<=3; i++) {
        $i = rev(tr($i,"ACGT","TGCA"))
    }
}
{ print }

function tr(instr,old,new,      outstr,pos,map) {
    for (pos=1; pos<=length(old); pos++) {
        map[substr(old,pos,1)] = substr(new,pos,1)
    }
    for (pos=1; pos<=length(instr); pos++) {
        outstr = outstr map[substr(instr,pos,1)]
    }
    return outstr
}

function rev(instr,     outstr,pos) {
    for (pos=1; pos<=length(instr); pos++) {
        outstr = substr(instr,pos,1) outstr
    }
    return outstr
}

.

$ awk -f tst.awk file
seq     T       G       change
seq     T       A       ok
seq     G       G       change
seq     GT      AGG     change

Answer 4

Если вас устраивает perl:

$ perl -F'\t' -lane 'if($F[3] eq "change") {
                     $F[1] = (reverse $F[1] =~ tr/ATGC/TACG/r);
                     $F[2] = (reverse $F[2] =~ tr/ATGC/TACG/r) }
                     print join "\t", @F' ip.txt
seq T   G   change
seq T   A   ok
seq G   G   change
seq GT  AGG change

Можно также использовать, но не c для столбцов, изменит любую последовательность ATCG символов:

perl -lpe 's/\t\K[ATCG]++(?=.*\tchange$)/reverse $&=~tr|ATGC|TACG|r/ge'