按序列大小对 fasta 进行排序

Question

我目前想按序列大小对 hudge fasta 文件（+10**8 行和序列）进行排序。 fasta 是生物学中用于存储序列（遗传或蛋白质）的明确定义格式：

>id1

sequence 1 # 可以在几行

>id2

序列 2

...

我有运行一个工具，可以 tsv 格式：

标识符、长度和标识符的字节位置。

现在我正在做的是按长度列对这个文件进行排序，然后我解析这个文件并使用 seek 来检索相应的序列，然后将它追加到一个新文件中。

# this fonction will get the sequence using seek
def get_seq(file, bites):  

    with open(file) as f_:
        f_.seek(bites, 0) # go to the line of interest
        line = f_.readline().strip() # this line is the begin of the 
                                     #sequence
        to_return = "" # init the string which will contains the sequence

        while not line.startswith('>') or not line:  # while we do not 
                                                     # encounter another identifiant
        to_return += line
        line = f_.readline().strip()

    return to_return
# simply append to a file the id and the sequence
def write_seq(out_file, id_, sequence):

    with open(out_file, 'a') as out_file:
        out_file.write('>{}\n{}\n'.format(id_.strip(), sequence))

# main loop will parse the index file and call the function defined below
with open(args.fai) as ref:

    indice = 0

    for line in ref:

        spt = line.split()
        id_ = spt[0]
        seq = get_seq(args.i, int(spt[2]))
        write_seq(out_file=args.out, id_=id_, sequence=seq)

我的问题是下面真的很慢这正常吗（需要好几天）？我有另一种方法吗？我不是一个纯粹的信息学家，所以我可能会漏掉一些要点，但我相信索引文件并使用搜索是实现这一目标的最佳方式，我错了吗？

Answer 1

似乎为每个序列打开两个文件可能对运行时间有很大贡献。您可以将文件句柄而不是文件名传递给 get/write 函数，但我建议使用已建立的 fasta parser/indexer，例如 biopython 或 samtools。这是一个使用 samtools 的（未经测试的）解决方案：

subprocess.call(["samtools", "faidx", args.i])
with open(args.fai) as ref:

    for line in ref:

        spt = line.split()
        id_ = spt[0]
        subprocess.call(["samtools", "faidx", args.i, id_, ">>", args.out], shell=True)

Answer 2

bash 和一些基本的 unix 命令（csplit 是线索）呢？我写了这个简单的脚本，但你可以 customize/improve 它。它没有高度优化，也不使用索引文件，但可能运行更快。

csplit -z -f tmp_fasta_file_  '/>/' '{*}'

for file in tmp_fasta_file_*
do
  TMP_FASTA_WC=$(wc -l < $file | tr -d ' ')
  FASTA_WC+=$(echo "$file $TMP_FASTA_WC\n")
done

for filename in $(echo -e $FASTA_WC | sort -k2 -r -n | awk -F" " '{print }')
do
  cat "$filename" >> 
done

rm tmp_fasta_file*

第一个位置参数是您的 fasta 文件的文件路径，第二个是输出的文件路径，即 ./script.sh input.fasta output.fasta

Answer 3

使用修改后的 fastq-sort 版本（目前可在 https://github.com/blaiseli/fastq-tools), we can convert the file to fastq format using bioawk 上使用，使用我添加的 -L 选项排序，然后转换回 fasta:

cat test.fasta \
    | tee >(wc -l > nb_lines_fasta.txt) \
    | bioawk -c fastx '{l = length($seq); printf "@"$name"\n"$seq"\n+\n%.*s\n", l, "IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII"}' \
    | tee >(wc -l > nb_lines_fastq.txt) \
    | fastq-sort -L \
    | tee >(wc -l > nb_lines_fastq_sorted.txt) \
    | bioawk -c fastx '{print ">"$name"\n"$seq}' \
    | tee >(wc -l > nb_lines_fasta_sorted.txt) \
    > test_sorted.fasta

fasta -> fastq 转换步骤非常难看。我们需要生成与序列长度相同的虚拟 fastq 质量。我发现没有比这个基于 https://www.gnu.org/software/gawk/manual/html_node/Format-Modifiers.html#Format-Modifiers.

末尾提到的 "dynamic width" 事情更好的方法来使用 (bio)awk 了

IIIII...字符串应该比最长的输入序列长，否则会得到无效的fastq，当转换回fasta时，bioawk似乎默默地跳过了这些无效的读取。

在上面的例子中，我添加了计算行数的步骤。如果行号不连贯，可能是因为 IIIII... 字符串太短。

生成的 fasta 文件将首先具有较短的序列。要在文件顶部获得最长的序列，请将 -r 选项添加到 fastq-sort.

请注意 fastq-sort 将中间文件写入 /tmp。如果由于某种原因在擦除它们之前中断，您可能需要手动清理 /tmp 而不是等待下一次重启。

编辑

我实际上找到了一种更好的方法来生成与序列长度相同的虚拟质量：简单地使用序列本身：

cat test.fasta \
    | bioawk -c fastx '{print "@"$name"\n"$seq"\n+\n"$seq}' \
    | fastq-sort -L \
    | bioawk -c fastx '{print ">"$name"\n"$seq}' \
    > test_sorted.fasta

这个解决方案更干净（并且稍微快一些），但我保留上面的原始版本，因为 printf 的 "dynamic width" 功能和 tee 用于检查中间数据长度的用法可能很有趣。

Answer 4

你也可以很方便的用awk，查看下面的代码：

awk '/^>/ {printf("%s%s\t",(N>0?"\n":""),[=10=]);N++;next;} {printf("%s",[=10=]);} END {printf("\n");}'  input.fasta  |\
awk -F '\t' '{printf("%d\t%s\n",length(),[=10=]);}' |\
sort -k1,1n | cut -f 2- | tr "\t" "\n"

此方法和其他方法已在 Biostars 中发布（例如使用 BBMap 的 sortbyname.sh 脚本），我强烈推荐此社区解决此类问题。

按序列大小对 fasta 进行排序

sort fasta by sequence size

sorting

bioinformatics

fasta

python-3.x

编辑