如何从fasta文件中移除所有N序列条目

Question

我想要删除一个fasta文件中所有碱基都是N的条目，但保留那些包含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

Answer 1

使用 shell 命令 egrep（就是 grep 加上正则表达式）

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 的标题行。

Answer 2

在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
#...

希望这对你有帮助。

Answer 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)

Answer 4

使用GNU awk

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

Answer 5

在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个字符）。