☰

将带有脚手架的 fasta 文件分成相同长度的文件，并尊重脚手架 ID 和序列

tag-icon • tag-icon awk grep split csplit
将带有脚手架的 fasta 文件分成相同长度的文件，并尊重脚手架 ID 和序列

我目前正在处理一个大的 fasta 文件（3.7GB），其中有脚手架。每个支架都有一个唯一的标识符，从>第一行开始，在连续的行上，它具有如下所示的 DNA 序列：

>9999992:0-108
AAAGAATTGTATTCCCTCCAGGTAGGGGGGATAGTTGAGGGGATACATAG
TGGGAAGGCTTTTCATGCGGAGGGACTAGAATGTGCTCCCGACTGACAAA
GCAGCTTG
>9999993:0-118
AGGGACTAGAAATGAGATTAAAAAGAGTAAAAGCACTGATACAAGTACAA
AAACAAATTGCTTCACCTCCAAAACCCCAGAAACTGCCCCACTTGGCTCC
CATTTAACCTACCTTCAA
>9999994:0-113
CCATCCTCATCCTTTCCTCCCCATATCTTCCTCTGACCCCAAAGCTCAGG
TTTCCTGTCTTGTTTCCCAGAATCTGTACCTCATGGTAGTTAAACCTTCC
CCTCTGGCAGCCA
>9999997:0-87
AACATCCCTGTGGCCTGAGAGACTGCCAGCCACAGCGGTGACAGTCCCTG
CGAGAGGCTGCTGCAAAAAGACTGGAGAGAAAGCAGA
>9999998:0-100
AAACATCAGCGCCAAGTCCCCGAAACCAGCAGGGTCACTGGGCGGCCGGC
CTGAAATACCCCAGCAGGCCAGCAGTGCCGGGTGCCTGGGGAGGTGTCCT
>9999999:0-94
AAGAAACTTTTCCCTTAACCAATGAAGAGTTTTATGTAAAGGAAATTTAG
TAATTTTTTAAAAAATGGTAATGACAGATTTAAGTAATTTAATT

我想将文件分割成长度相同的小文件来处理它，但我需要同时尊重 ID 和序列，并获得如下所示的内容：

file1.fa
>9999992:0-108
AAAGAATTGTATTCCCTCCAGGTAGGGGGGATAGTTGAGGGGATACATAG
TGGGAAGGCTTTTCATGCGGAGGGACTAGAATGTGCTCCCGACTGACAAA
GCAGCTTG
>9999993:0-118
AGGGACTAGAAATGAGATTAAAAAGAGTAAAAGCACTGATACAAGTACAA
AAACAAATTGCTTCACCTCCAAAACCCCAGAAACTGCCCCACTTGGCTCC
CATTTAACCTACCTTCAA

file2.fasta
>9999994:0-113
CCATCCTCATCCTTTCCTCCCCATATCTTCCTCTGACCCCAAAGCTCAGG
TTTCCTGTCTTGTTTCCCAGAATCTGTACCTCATGGTAGTTAAACCTTCC
CCTCTGGCAGCCA
>9999997:0-87
AACATCCCTGTGGCCTGAGAGACTGCCAGCCACAGCGGTGACAGTCCCTG
CGAGAGGCTGCTGCAAAAAGACTGGAGAGAAAGCAGA

file3.fasta
>9999998:0-100
AAACATCAGCGCCAAGTCCCCGAAACCAGCAGGGTCACTGGGCGGCCGGC
CTGAAATACCCCAGCAGGCCAGCAGTGCCGGGTGCCTGGGGAGGTGTCCT
>9999999:0-94
AAGAAACTTTTCCCTTAACCAATGAAGAGTTTTATGTAAAGGAAATTTAG
TAATTTTTTAAAAAATGGTAATGACAGATTTAAGTAATTTAATT

请帮我。我尝试使用csplitgrep 但得到错误的输出。

答案1

下面的代码可能对您有帮助，但对于大文件来说可能会很慢（根据您的计算机规格）。

首先，您应该获得 的总数Unique Identifier。您可以通过使用来实现grep -c

total=$(grep -c "^>" largeFastaFile.txt)

上面的代码会将total以 开头的匹配行的计数分配给变量>。现在您应该得到Unique Identifier每个文件的数量。所以如果你想拥有10 个文件。你应该划分total/10：

max=$((total/10))
#If total has 3714529 then max will have 371,452.

最后，使用awk命令您可以将大文件分割为 10 个文件（实际上 11）大约有371,452 唯一标识符每个文件：

awk -v maxline=$max -v count=0\ 
 '{if(NR>1) { if( (NR-2)%maxline == 0 ) count++ ; print ">"$0 >("file"count".fasta")  } }'\
 RS='>' largeFastaFile.txt

该脚本应如下所示：

#!/bin/bash

total=$(grep -c "^>" largeFastaFile.txt)
max=$((total/10)) # where 10 is the number of files you will get

awk -v maxline=$max -v count=0 '{if(NR>1) { if( (NR-2)%maxline == 0 ) count++ ; print ">"$0 >("file"count".fasta")  } }' RS='>' largeFastaFile.txt

实际上你会得到 11 个文件，因为如果总数是371452910 个文件，并且371,452 Unique Identifier每个文件的数量相同，那么将有几行必须位于另一个文件中：

371,452 * 10 = 3,714,520

所以11号文件将有最后一个9 唯一标识符

答案2

我一直在使用一位老同事几年前编写的几个脚本来在 fasta 和 tbl（序列标识符，TAB，单行序列）格式之间进行转换。我已经发布了脚本这里。使用它们，您可以轻松地将文件转换为每行一个序列，然后分割指定每个文件中想要的行数，然后将分割文件转换回 fasta：

FastaToTbl file.fa | split -l 2 - splitFile

使用您的示例输入，创建这些文件：

$ for f in splitFile*; do printf "\n===== File: %s ====\n" "$f"; cat "$f"; done

===== File: splitFileaa ====
9999992:0-108   AAAGAATTGTATTCCCTCCAGGTAGGGGGGATAGTTGAGGGGATACATAGTGGGAAGGCTTTTCATGCGGAGGGACTAGAATGTGCTCCCGACTGACAAAGCAGCTTG
9999993:0-118   AGGGACTAGAAATGAGATTAAAAAGAGTAAAAGCACTGATACAAGTACAAAAACAAATTGCTTCACCTCCAAAACCCCAGAAACTGCCCCACTTGGCTCCCATTTAACCTACCTTCAA

===== File: splitFileab ====
9999994:0-113   CCATCCTCATCCTTTCCTCCCCATATCTTCCTCTGACCCCAAAGCTCAGGTTTCCTGTCTTGTTTCCCAGAATCTGTACCTCATGGTAGTTAAACCTTCCCCTCTGGCAGCCA
9999997:0-87    AACATCCCTGTGGCCTGAGAGACTGCCAGCCACAGCGGTGACAGTCCCTGCGAGAGGCTGCTGCAAAAAGACTGGAGAGAAAGCAGA

===== File: splitFileac ====
9999998:0-100   AAACATCAGCGCCAAGTCCCCGAAACCAGCAGGGTCACTGGGCGGCCGGCCTGAAATACCCCAGCAGGCCAGCAGTGCCGGGTGCCTGGGGAGGTGTCCT
9999999:0-94    AAGAAACTTTTCCCTTAACCAATGAAGAGTTTTATGTAAAGGAAATTTAGTAATTTTTTAAAAAATGGTAATGACAGATTTAAGTAATTTAATT

-l 2将每个文件更改为您想要的序列数。然后，遍历分割文件并将它们每个转换为 fasta：

c=0
for f in splitFile*; do 
    outFile=file.$((++c)).fa
    TblToFasta "$f" > "$outFile"
done

这会产生以下文件：

$ for f in file.*.fa; do printf "\n===== File: %s ====\n" "$f"; cat "$f"; done

===== File: file.1.fa ====
>9999992:0-108 
AAAGAATTGTATTCCCTCCAGGTAGGGGGGATAGTTGAGGGGATACATAGTGGGAAGGCT
TTTCATGCGGAGGGACTAGAATGTGCTCCCGACTGACAAAGCAGCTTG
>9999993:0-118 
AGGGACTAGAAATGAGATTAAAAAGAGTAAAAGCACTGATACAAGTACAAAAACAAATTG
CTTCACCTCCAAAACCCCAGAAACTGCCCCACTTGGCTCCCATTTAACCTACCTTCAA

===== File: file.2.fa ====
>9999994:0-113 
CCATCCTCATCCTTTCCTCCCCATATCTTCCTCTGACCCCAAAGCTCAGGTTTCCTGTCT
TGTTTCCCAGAATCTGTACCTCATGGTAGTTAAACCTTCCCCTCTGGCAGCCA
>9999997:0-87 
AACATCCCTGTGGCCTGAGAGACTGCCAGCCACAGCGGTGACAGTCCCTGCGAGAGGCTG
CTGCAAAAAGACTGGAGAGAAAGCAGA

===== File: file.3.fa ====
>9999998:0-100 
AAACATCAGCGCCAAGTCCCCGAAACCAGCAGGGTCACTGGGCGGCCGGCCTGAAATACC
CCAGCAGGCCAGCAGTGCCGGGTGCCTGGGGAGGTGTCCT
>9999999:0-94 
AAGAAACTTTTCCCTTAACCAATGAAGAGTTTTATGTAAAGGAAATTTAGTAATTTTTTA
AAAAATGGTAATGACAGATTTAAGTAATTTAATT

使用的两个脚本是：

  • 快速转Tbl

      #!/usr/bin/awk -f
  {
          if (substr($1,1,1)==">")
          if (NR>1)
                      printf "\n%s\t", substr($0,2,length($0)-1)
          else 
              printf "%s\t", substr($0,2,length($0)-1)
          else 
                  printf "%s", $0
  }END{printf "\n"}

  • TblToFasta

      #! /usr/bin/awk -f
  {
    sequence=$NF

    ls = length(sequence)
    is = 1
    fld  = 1
    while (fld < NF)
    {
       if (fld == 1){printf ">"}
       printf "%s " , $fld

       if (fld == NF-1)
        {
          printf "\n"
        }
        fld = fld+1
    }
    while (is <= ls)
    {
      printf "%s\n", substr(sequence,is,60)
      is=is+60
    }
  }

将它们保存在您的文件中$PATH并使它们可执行。

相关内容