如何从 fasta 文件中删除所有 N 序列条目答案

【问题标题】：How to remove all-N sequence entries from fasta file(s)如何从 fasta 文件中删除所有 N 序列条目
【发布时间】：2014-05-18 04:16:13
【问题描述】：

我想删除所有核苷酸均为 N 的 fasta 文件的条目，但不删除包含 ACGT 和 N 核苷酸的条目。

来自输入文件内容的示例：

#>seq_1
TGCTAGCTAGCTGATCGTGTCGATCG
CACCACANNNNNCACGTGTCG
#>seq2
NNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNN
#>seq3
catgcatcgacgatgctgacgatc
#>seq4
cacacaccNNNNttgtgca
#...

希望输出文件内容为：

#>seq_1
TGCTAGCTAGCTGATCGTGTCGATCG
CACCACANNNNNCACGTGTCG
#>seq3
catgcatcgacgatgctgacgatc
#>seq4
cacacaccNNNNttgtgca
#...

在使用 awk、perl、python 或其他方法时有什么建议吗？谢谢！ FDS

【问题讨论】：

标签： python perl awk

【解决方案1】：

使用 GNU awk

awk -v RS='#>seq[[:digit:]]+' '!/^[N\n]+$/{printf "%s",term""$0}; {term=RT}' file
#>seq_1
TGCTAGCTAGCTGATCGTGTCGATCG
CACCACANNNNNCACGTGTCG
#>seq3
catgcatcgacgatgctgacgatc
#>seq4
cacacaccNNNNttgtgca

【讨论】：

【解决方案2】：

使用 python，使用 BioPython 模块：

import Bio

INPUT  = "bio_input.fas"
OUTPUT = "bio_output.fas"

def main():
    records = Bio.SeqIO.parse(INPUT, 'fasta')
    filtered = (rec for rec in records if any(ch != 'N' for ch in rec.seq))
    Bio.SeqIO.write(filtered, OUTPUT, 'fasta')

if __name__=="__main__":
    main()

但是请注意，FastA 规范说序列 ID 应该以 > 开头，而不是 #>！

对抗

>seq_1
TGCTAGCTAGCTGATCGTGTCGATCG
CACCACANNNNNCACGTGTCG
>seq2
NNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNN
>seq3
catgcatcgacgatgctgacgatc
>seq4
cacacaccNNNNttgtgca

这会产生

>seq_1
TGCTAGCTAGCTGATCGTGTCGATCGCACCACANNNNNCACGTGTCG
>seq3
catgcatcgacgatgctgacgatc
>seq4
cacacaccNNNNttgtgca

（注意，默认换行长度为 60 个字符）。

【讨论】：

【解决方案3】：

所以基本上块以 #> 标记开头，并且您希望删除没有任何行包含除 N 以外的任何内容的块。Python 中的一种方法：

#! /usr/bin/env python
import fileinput, sys, re

block=[]
nonN=re.compile('[^N\n]')
for line in fileinput.input():
    if line.startswith('#>'):
        if len(block)==1 or any(map(nonN.search, block[1:])):
            sys.stdout.writelines(block)
        block=[line]
    else:
        block.append(line)
if len(block)==1 or any(map(nonN.search, block[1:])):
    sys.stdout.writelines(block)

【讨论】：

【解决方案4】：

在 python 中使用正则表达式：

#!/usr/bin/env python
import re
ff = open('test', 'r')
data = ff.read()
ff.close()

m = re.compile(r'(#>seq\d+[N\n]+)$', re.M)
f = re.sub(m, '', data)

fo = open('out', 'w')
fo.write(f)
fo.close()

然后你会进入你的输出文件：

#>seq_1
TGCTAGCTAGCTGATCGTGTCGATCG
CACCACANNNNNCACGTGTCG

#>seq3
catgcatcgacgatgctgacgatc
#>seq4
cacacaccNNNNttgtgca
#...

希望这会有所帮助。

【讨论】：

【解决方案5】：

使用 shell 命令 egrep (grep + regrex)

egrep -B 1 "^[^NNNN && ^#seq]" your.fa >conert.fa
# -B 1 means print match row and the previous row
# ^[^NNNN && ^#seq] is regrex pattern means not match with (begin with NNNN and 
# #seq)

所以只匹配以普通 A/T/G/C 序列行开头的行及其前一行，即 fasta 标题

【讨论】：

虽然这段代码可能会解决问题，但一个好的答案还应该解释代码的作用以及它如何提供帮助。