删除小于某一长度的行及其上的行(删除FASTA文件中的短序列)

Question

我有一个包含以下文本的文件:

>seq1
GAAAT
>seq2
CATCTCGGGA
>seq3
GAC
>seq4
ATTCCGTGCC

如果一行不是以">"开头短于5个字符，我想删除它和它右边的一个

预期输出:

>seq2
CATCTCGGGA
>seq4
ATTCCGTGCC

我尝试过sed -r '/^.{,5}$/d'，但它也删除了带有">"的行。

Answer 1

使用sed

$ sed '/^>/N;/n[A-Z]{6,}$/!d' input_file
>seq2
CATCTCGGGA
>seq4
ATTCCGTGCC

Answer 2

对于您所显示的示例，对于awk，您可以尝试遵循awk代码。简单的解释是，如果一行从>开始，那么用当前行和next的值设置变量val将跳过这里的所有进一步的语句。然后，如果一行不是从>开始，如果当前行长度大于5，则打印val ORS和当前行。

awk '/^>/{val=$0;next} length($0)>5{print val ORS $0}' Input_file

Answer 3

对于GNU sed，您可以使用

sed -E '/>/N;/n[^>].{0,4}$/d'

细节:

• />/-查找包含>的行(如果必须在开头，则在>之前添加^)
• N-读取该行并将其附加到模式空间，并使用前导换行
• n[^>].{0,4}$-换行，>以外的字符(因为第一个字符不应该是>)，然后0到4个字符，直到字符串结束
• d删除模式空间中的值。

查看在线演示:

#!/bin/bash
s='>seq1
GAAAT
>seq2
CATCTCGGGA
>seq3
GAC
>seq4
ATTCCGTGCC'
sed -E '/>/N;/n[^>].{0,4}$/d' <<< "$s"

输出:

>seq2
CATCTCGGGA
>seq4
ATTCCGTGCC

Answer 4

不要重新发明轮子。使用常见的生物信息学工具，如seqtk或seqkit。除此之外，这些工具还可以处理多行FASTA序列。例子:

seqtk seq -L 5 in.fasta > out.fasta
seqkit seq -m 5 in.fa > out.fasta

安装这些工具，使用conda，特别是miniconda，例如:

conda create --channel bioconda --name seqtk seqtk
conda activate seqtk
# ... use seqtk here ...
conda deactivate

引用:

从FASTA文件中删除序列<300个碱基:https://www.biostars.org/p/329680/
seqtk: https://github.com/lh3/seqtk
seqkit: https://bioinf.shenwei.me/seqkit/
conda: https://docs.conda.io/projects/conda/en/latest/user-guide/install/index.html

Answer 5

非常不优雅的解决方案只是为了让它工作=(

{m,g}awk ' !_<NR && 5 < length($(NF-=_==$NF))
&&                           
($!NF =">seq" $!_ ORS $+NF)^_'  FS='[ tn]+' OFS='n' RS='n*>seq'

>seq2
CATCTCGGGA
>seq4
ATTCCGTGCC

Answer 6

我推荐Biopython。它使FASTA处理方便。

'''
Biopython version: 1.79
'''
from Bio import SeqIO
for seq_record in SeqIO.parse("input.fasta", "fasta"):
if len(seq_record.seq)>=5:
print(">",seq_record.id)
print(seq_record.seq)