#在宏基因组/宏转录组数据进行组装后,常需要去除短片段,筛选出较长的片段以供后续分析#
在Linux系统中,您可以使用一些文本处理工具来提取长度大于n的序列。其中,常用的工具之一是awk
命令,它可以用于处理文本文件并提取符合条件的行。
以下是在Linux系统中使用awk
命令来提取长度大于1000的序列的示例命令:
awk '/^>/ {if (seqlen > 1000) {if (seqname != "") {print seqname; print seq} } seq=""; seqlen=0; seqname=$0; next} {seq = seq $0; seqlen += length($0)} END {if (seqlen > 1000) {if (seqname != "") {print seqname; print seq} } }' final_assembled.fasta > extracted_sequences.fasta
详细的:
awk '/^>/ {
if (seqlen > 1000) {
if (seqname != "") {
print seqname; # 打印上一个序列的名称
print seq; # 打印上一个序列的内容
}
}
seq = ""; # 重置当前序列内容为空
seqlen = 0; # 重置当前序列长度为0
seqname = $0; # 保存当前序列名称行,包括 ">" 符号
next # 跳过后续的处理并继续到下一行
}
{
seq = seq $0; # 将当前行的内容追加到当前序列内容中
seqlen += length($0); # 计算当前序列的长度
}
END {
if (seqlen > 1000) {
if (seqname != "") {
print seqname; # 打印最后一个序列的名称
print seq; # 打印最后一个序列的内容
}
}
}' final_contigs.fasta > extracted_sequences.fasta
代码的执行流程如下:
-
awk
是一个用于文本处理的命令行工具。我们用/^>/
来匹配以 ">" 开头的行,即序列名称行。 -
当遇到以 ">" 开头的行时,表示遇到了序列名称行。接下来执行大括号内的代码块。
-
在代码块中,首先执行
if (seqlen > 1000)
条件判断,检查上一个序列的长度是否大于1000。seqlen
是一个变量,用于保存当前序列的长度。 -
如果上一个序列的长度大于1000,则进入内层的
if (seqname != "")
条件判断。seqname
是一个变量,用于保存上一个序列的名称行。 -
如果上一个序列的名称行非空,则执行
print seqname; print seq;
,打印上一个序列的名称和内容。seq
是一个变量,用于保存上一个序列的内容。 -
然后,执行
seq = ""; seqlen = 0;
,将当前序列的内容和长度重置为空,为下一个序列做准备。 -
在序列名称行的处理完成后,使用
next
命令跳过后续的处理,并继续到下一行。 -
对于不是以 ">" 开头的行,表示是序列的内容行。程序将当前行的内容追加到
seq
变量中,同时增加当前序列的长度seqlen
。 -
当处理完所有行后,在
END
部分进行最后的处理。这部分代码在整个文件处理完毕后执行一次。 -
在
END
部分,首先执行if (seqlen > 1000)
条件判断,检查最后一个序列的长度是否大于1000。 -
如果最后一个序列的长度大于1000,则进入内层的
if (seqname != "")
条件判断。 -
如果最后一个序列的名称行非空,则执行
print seqname; print seq;
,打印最后一个序列的名称和内容。 -
最终,通过重定向操作符
>
,将提取的具有序列名称的序列保存到 "extracted_sequences.fasta" 文件中