☰

从 fasta 文件中提取序列

tag-icon • tag-icon text-processing
从 fasta 文件中提取序列

我有一个 fasta 文件（格式不正确），其中包含数十万个不同长度的 DNA 序列，如下所示：

>NODE_213384_length_62_cov_8686_ID_2134025ATCGAATGGAATCATCGAATGGACTCGAATGGAATAATCATTGAACGGAATCGAATGG>NODE_213385_length_62_cov_7933_ID_2134027ATCATCATCGAATGGAATCGAATGGAATCATCGAATGGACTCGAATGGAATAATCATTGAAC>NODE_213386_length_62_cov_7184_ID_2134029AATGATTATTCCATTCGAGTCCATTCGATGATTCCATTCGATTCCATTCGATGATGATTGCA>NODE_213387_length_62_cov_8639_ID_2134031CAGAGCAGACTTGAAACACTCTTTTTGTGGAATTTGCAAGTGGAGATTTCAGCCGCTTTGAG>NODE_213388_length_62_cov_6833_ID_2134033AGACTTGAAACACTCTTTTTGTGGAATTTGCAAGTGGAGATTTCAGCCGCTTTGAGGTCAAT

我想使用一个简单的 Linux 命令来提取那些长度超过 1000bp 的序列，并以正确的 fasta 格式输出，如下所示：

>NODE_213384_length_62_cov_8686_ID_2134025
ATCGAATGGAATCATCGAATGGACTCGAATGGAATAATCATTGAACGGAATCGAATGG
>NODE_213385_length_62_cov_7933_ID_2134027
ATCATCATCGAATGGAATCGAATGGAATCATCGAATGGACTCGAATGGAATAATCATTGAAC
>NODE_213386_length_62_cov_7184_ID_2134029
AATGATTATTCCATTCGAGTCCATTCGATGATTCCATTCGATTCCATTCGATGATGATTGCA
>NODE_213387_length_62_cov_8639_ID_2134031
CAGAGCAGACTTGAAACACTCTTTTTGTGGAATTTGCAAGTGGAGATTTCAGCCGCTTTGAG

感谢任何能提供帮助的人。

答案1

看起来你的队伍很长。您可能会遇到 sed 和 awk 的问题，因为它们会耗尽可用内存来处理此任务（当然取决于您的文件大小/行长度）。因此，采用两阶段方法，但由于内存限制，您可以使用tr上面的 awk、perl 或 sed 解决方案之一。

head -20 inputfile | 
tr '>' '\n'  > stage1
perl -ne 'print ">$1 $2\n" if /^(.*?)([ACGTU]+)$/ && length($2)>1000' < stage1 > output

一旦您对前 20 行的效果感到满意，就可以在单个管道中真正执行此操作：

tr '>' '\n' inputfile | 
perl -ne 'print ">$1 $2\n" if /^(.*?)([ACGTU]+)$/ && length($2)>1000' > output

我的 Perl 脚本可能不如其他脚本高效，但应该可以完成工作。为了清楚起见，我写了它。当且仅当存在包含超过 1000 个碱基对的行时，打印出标签，后跟一个空格，然后是相关的碱基对和换行符。

答案2

“简单的 Linux 命令”可能类似于：

sed 's/\(>[^ACGT]*_[0-9]\+\)\([ACGT]\+\)/\1\n\2\n/g' yourdnafile |egrep -B1 '^[ACGT]{1000}'

sed 部分将每组分成 2 行，grep 显示并匹配超过 1000 行及其前面的行 (-B1)。

或者这可能更简单：

sed 's/\(>[^>ACGT]*\)\([ACGT]\{1000\}[ACGT]*\)/\1\n\2\n/g;s/>[^>ACGT\n]*[ACGT]\+//g' yourdnafile

答案3

您可以分两步完成；首先分割成目标格式，然后打印DNA序列足够长的线对。例如，假设1000bp指的是 DNA 序列字符长度（而不是“长度”一词后面的值），它可能是以下内容

cat inputfile | sed -e 's/>/\n>/g' -e 's/\(ID_[0-9]*\)/\n\1/g' |\
awk -F_ '$1=="NODE" { n=$0; next; } length($0)>1000 { print n; print ">" $0;}' \
> outputfile

相关内容