☰

尝试编写一个脚本从 6 个不同的框架中查找 ORF

tag-icon • tag-icon command-line regular-expression bioinformatics
尝试编写一个脚本从 6 个不同的框架中查找 ORF

所以基本上，你有一个序列，比如说

AAAGCATATGCTAGCCCGTATAGCGATACTAGCTATACGATATATATGATCAATGCCCGTATAG

您需要找到序列中的所有“ORF”，ORF 是以 ATG 开头并以 TGA 或 TAA 或 TAG 结尾的任何区域。

因此，例如，在上面的序列中，ORF 是

AAAGCAT**ATGCTAG**CCCGTATAGCGATACTAGCTATACGATATAT**ATGATCAATGCCCGTATAG**

你会注意到，在第二个ORF中，ORF内部有一个ATG，当它结束之前没有TGA或TAA或TAG时就会发生这种情况。

是的，基本上这就是问题。我知道在 C++ 上有 2-3 种方法，但就 bash 的语法而言我不知道。

我不能使用任何库或 perl 或类似的东西，没有特殊的函数，只是 grep、awk、sed 和循环之类的东西。

答案1

您可以使用grep 如果您正在使用 GNU 版本的grep，它可以-P选择与 perl 兼容的正则表达式 (PCRE)。

或者，您可以使用作者的pcregrep（又名）pgrep聚合酶链式反应图书馆。现在没有太多理由使用它，除非您使用旧版本的 GNU grep 或非 GNU grep，它们不支持该-P选项并且无法升级或替换。

例如（假设序列位于名为 的文件中input.txt）：

$ grep -oP 'ATG.*?TA[AG]'  input.txt 
ATGCTAG
ATGATCAATGCCCGTATAG

该-o选项告诉 GNU grep 仅输出匹配的文本，而不是整行，并-P告诉它使用与 perl 兼容的正则表达式。

或者，如果您正在使用pcregrep：

$ pcregrep -o 'ATG.*?TA[AG]'  input.txt 
ATGCTAG
ATGATCAATGCCCGTATAG

?正则表达式 ( ) 中的非贪婪修饰符.*?确保它捕获全部匹配模式，而不仅仅是最长的模式。在正则表达式的上下文中，“贪婪”意味着“尝试尽可能多地匹配”（默认值），“非贪婪”意味着“尝试尽可能少地匹配”。

对此有一个很好的解释https://www.regular-expressions.info/repeat.html。顺便说一句，该网站的其余部分是学习正则表达式的好地方，有很多教程和示例。

请注意，大多数正则表达式库不要实现非贪婪匹配，它是一个 Perl 扩展，也已被 GNU grep 采用。以及与以下链接的程序聚合酶链式反应， 当然。

顺便说一句，这就是输出的样子没有非贪婪修饰符：

$ grep -oP 'ATG.*TA[AG]'  input.txt 
ATGCTAGCCCGTATAGCGATACTAGCTATACGATATATATGATCAATGCCCGTATAG

答案2

可以使用正则表达式进行搜索和匹配perl（我最喜欢的sed不支持所需的非贪婪正则表达式（p）匹配）：

# echo AAAGCATATGCTAGCCCGTATAGCGATACTAGCTATACGATATATATGATCAATGCCCGTATAG | perl -pe 's|(.*?)(ATG.*?TA[AG])(.*?)|\2\n|g'
ATGCTAG
ATGATCAATGCCCGTATAG
AA

其中perl命令打印您请求的序列 - 以及最后一行和该行的其余部分。为了不使正则表达式变得更复杂，可以手动或使用以下命令轻松删除它head：

echo AAAGCATATGCTAGCCCGTATAGCGATACTAGCTATACGATATATATGATCAATGCCCGTATAGAA | perl -pe 's|(.*?)(ATG.*?TA[AG])(.*?)|\2\n|g' | head -n -1
ATGCTAG
ATGATCAATGCCCGTATAG

这些echo命令显示序列的正确结果。如果您希望从一个文件读取内容，并将结果发送到第二个文件，请执行以下操作：

cat original_file | perl -pe 's|(.*?)(ATG.*?TA[AG])(.*?)|\2\n|g' | head -n -1 > new_file

original_file您的源文件和new_file带有过滤模式的目标文件在哪里。

上面使用head来自 GNU coreutils 的“-n -1”语法。如果这对您不起作用，请尝试

cat original_file | perl -pe 's|(.*?)(ATG.*?TA[AG])(.*?)|\2\n|g' | awk 'NR>1 {print prev} {prev=$0}' > new_file

答案3

这在 bash 中是可能的，但它确实不是一个好的工具：

#!/bin/bash 

# Read the sequence into the variable $seq
seq=$1

## Check all three frames
for ((frame=0; frame<=3; frame++)); do
  ## Read the sequence in groups of 3
  for ((i=$frame;i<=${#seq};i+=3)); do
    ## The codon: three nucleotides starting from the current position.
    codon=${seq:i:3}
    ## set isORF to 1 if this is an ATG
    if [[ ${seq:i:3} = "ATG" ]]; then 
      isORF=1
    fi
    ## If we're in an ORF
    if [[ $isORF = 1 ]]; then
      ## Add this codon to the ORF's sequence
      ORF="${ORF}${codon}"
      ## Is this a STOP?
      if [[ ${seq:i:3} = "TGA" ||
            ${seq:i:3} = "TAA" ||
            ${seq:i:3} = "TAG" ]];
      then
        ## Print the ORF
        echo "ORF: $ORF"
        ## Set isORF to 0 and empty the ORF variable to repeat the process
        isORF=0
        ORF=""
    fi
  fi
  done
done

将其另存为foo.sh，使其可执行 ( chmod a+x foo.sh) 并像这样运行它：

/path/to/foo.sh AAAGCATATGCTAGCCCGTATAGCGATACTAGCTATACGATATATATGATCAATGCCCGTATAG

相关内容