按GC含量进行入库序列读取

Question

我想"bin“(分割成单独的文件)多fasta核苷酸序列文件(例如，Roche-454运行约500,000个读取平均读取长度为250bp)。我想要的垃圾桶基于GC内容的每一个阅读。结果输出将是8个多fasta文件：

GC含量<20%

GC含量21-30%

GC含量31-40%

GC含量41-50%

GC含量51-60%

GC含量61-70%

GC含量71-80%

80% GC含量

有没有人知道已经有脚本或程序这样做了？如果没有，有人可以建议如何根据GC内容对多fasta文件进行排序(然后我可以将其拆分到相关的存储箱中)？

Answer 1

在R/ Bioconductor中，任务是(a)加载适当的文库(b)读取fasta文件(c)计算核苷酸使用和gc % (d)将数据剪切到bin中，(e)将原始数据输出到单独的文件中。沿着下面的思路

## load
library(ShortRead)
## input
fa = readFasta("/path/to/fasta.fasta")
## gc content. 'abc' is a matrix, [, c("G", "C")] selects two columns
abc = alphabetFrequency(sread(fa), baseOnly=TRUE)
gc = rowSums(abc[,c("G", "C")]) / rowSums(abc)
## cut gc content into bins, with breaks at seq(0, 1, .2)
bin = cut(gc, seq(0, 1, .2))
## output, [bin==lvl] selects the reads whose 'bin' value is lvl
for (lvl in levels(bin)) {
    fileName = sprintf("%s.fasta", lvl)
    writeFasta(fa[bin==lvl], file=fileName)
}

要开始使用R/ Bioconductor，请参阅http://bioconductor.org/install。所示大小的454个数据的内存需求并不太差，并且这里的脚本将相当快(例如，260k读取需要7秒)。

Answer 2

我建议使用Python和Biopython或者Perl和Bioperl来读入FASTA文件。有一个脚本可以计算Bioperl here中序列的C含量，而Biopython有一个function for it。然后，只需将每个序列的GC内容存储在字典或散列中，并遍历每个序列，然后根据GC内容的高度将它们写入文件中。

你需要更具体的帮助吗？

Answer 3

如果我没弄错的话，你需要如下代码(Python)：

def GC(seq): # determine the GC content
    s = seq.upper()
    return 100.0 * (s.count('G') + s.count('C')) / len(s)

def bin(gc): # get the number of the 'bin' for this value of GC content
    if gc < 20: return 1
    elif gc > 80: return 8
    else:
        return int(gc/10)

然后，您只需从文件中读取条目，计算GC内容，找到正确的bin，并将条目写入相应的文件。下面的示例使用我们在实验室中使用的Python package实现了这一点：

from pyteomics import fasta

def split_to_bin_files(multifile):
"""Reads a file and writes the entries to appropriate 'bin' files.
`multifile` must be a filename (str)"""

    for entry in fasta.read(multifile):
        fasta.write((entry,), (multifile+'_bin_'+
                    str(bin(GC(entry[1])))))

然后你就像split_to_bin_files('mybig.fasta')一样叫它。