如何重新格式化 fasta 文件头的字段并折叠序列？答案

【问题标题】：How can I Re-Format the fields of a fasta file header and fold the sequence?如何重新格式化 fasta 文件头的字段并折叠序列？
【发布时间】：2020-12-22 02:24:45
【问题描述】：

生物信息学中最常用的文件之一是fasta file format。
Fasta 文件很简单：它们包含开头的“标题”记录带有“>”，后跟“Sequence”记录，即后面的所有内容标题但在下一个记录分隔符之前（即“>”）。

>ENSP00000488314.1 pep chromosome:GRCh38:X:143884071:143885255:1 gene:ENSG00000276380.2 transcript:ENST00000618570.1 gene_biotype:polymorphic_pseudogene transcript_biotype:polymorphic_pseudogene gene_symbol:UBE2NL description:ubiquitin conjugating enzyme E2 N like (gene/pseudogene) [Source:HGNC Symbol;Acc:HGNC:31710]
MAELPHRIIKETQRLLAEPVPGIKAEPDESNARYFHVVIAGESKDSPFEGGTFKRELLLA
EEYPMAAPKVRFMTKIYHPNVDKLERIS*DILKDKWSPALQIRTVLLSIQALLNAPNPDD
PLANDVVEQWKTNEAQAIETARAWTRLYAMNSI
> next record...
> another one...

标头可以非常简单（例如，">ENSP00000488314.1"）或复杂。
复杂的标题重要但可变的信息。
对于上面的示例序列（来自ENSEMBL），标题记录由以下部分组成：

Field 01: ENSP00000488314.1                         <=Protein ID
Field 02: pep                                       <=Peptide record
Field 03: chromosome:GRCh38:X:143884071:143885255:1 <=Chromosome and chromosomal coordinates
Field 04: gene:ENSG00000276380.2                    <=Gene ID
Field 05: transcript:ENST00000618570.1              <=Transcript ID
Field 06: gene_biotype:polymorphic_pseudogene       <=Gene Biotype
Field 07: transcript_biotype:polymorphic_pseudogene <=Transcript Biotype
Field 08: gene_symbol:UBE2NL                        <=Gene Symbol
Up to here the fields are all nicely separated by spaces, and then...Field 09 (Variable)
Field 09: description:ubiquitin conjugating enzyme E2 N like (gene/pseudogene)
Field 10: [Source:HGNC Symbol;Acc:HGNC:31710]       <=Predictable

很多时候 long 标头不会被其他生物信息学应用程序很好地接受，因此需要缩短标头。
如果能以一种聪明的方式做到这一点，那就太好了。因此，使用 AWK 并使用下面的示例序列，我想：

第一：控制表头记录的打印如下：

始终保留第一个字段：

>ENSP00000488314.1

但随后可以省略和/或包含字段。示例：

>ENSP00000488314.1 gene:ENSG00000276380.2 transcript:ENST00000618570.1
Field: 01          04                     05
>ENSP00000488314.1 pep chromosome:GRCh38:X:143884071:143885255:1 [Source:HGNC Symbol;Acc:HGNC:31710]
Field: 01          02  03                                        10

为简单起见，完全忽略字段 09 是完全可以接受的，但能够使用字段 10 会很好

然后可以将序列“折叠”到用户指定的数字。例如，每 60 个字符折叠一次的记录：

>ENSP00000441696.1 pep chromosome:GRCh38:14:21868839:21869365:1 gene:ENSG00000211788.2 transcript:ENST00000390436.2 gene_biotype:TR_V_gene transcript_biotype:TR_V_gene gene_symbol:TRAV13-1 description:T cell receptor alpha variable 13-1 [Source:HGNC Symbol;Acc:HGNC:12108]
MTSIRAVFIFLWLQLDLVNGENVEQHPSTLSVQEGDSAVIKCTYSDSASNYFPWYKQELG
KGPQLIIDIRSNVGEKKDQRIAVTLNKTAKHFSLHITETQPEDSAVYFCAAS
>ENSP00000488314.1 pep chromosome:GRCh38:X:143884071:143885255:1 gene:ENSG00000276380.2 transcript:ENST00000618570.1 gene_biotype:polymorphic_pseudogene transcript_biotype:polymorphic_pseudogene gene_symbol:UBE2NL description:ubiquitin conjugating enzyme E2 N like (gene/pseudogene) [Source:HGNC Symbol;Acc:HGNC:31710]
MAELPHRIIKETQRLLAEPVPGIKAEPDESNARYFHVVIAGESKDSPFEGGTFKRELLLA
EEYPMAAPKVRFMTKIYHPNVDKLERIS*DILKDKWSPALQIRTVLLSIQALLNAPNPDD
PLANDVVEQWKTNEAQAIETARAWTRLYAMNSI
>ENSP00000437680.2 pep chromosome:GRCh38:22:42140203:42141924:-1 gene:ENSG00000205702.11 transcript:ENST00000435101.1 gene_biotype:polymorphic_pseudogene transcript_biotype:nonsense_mediated_decay gene_symbol:CYP2D7 description:cytochrome P450 family 2 subfamily D member 7 (gene/pseudogene) [Source:HGNC Symbol;Acc:HGNC:2624]
DPAQPPRDLTEAFLAKKEKAKGSPESSFNDENLRIVSVSNRRSTT

可以变成（序列每 120 个字符折叠一次）：

>ENSP00000441696.1 gene:ENSG00000211788.2 transcript:ENST00000390436.2
MTSIRAVFIFLWLQLDLVNGENVEQHPSTLSVQEGDSAVIKCTYSDSASNYFPWYKQELGKGPQLIIDIRSNVGEKKDQRIAVTLNKTAKHFSLHITETQPEDSAVYFCAAS
>ENSP00000488314.1 gene:ENSG00000276380.2 transcript:ENST00000618570.1
MAELPHRIIKETQRLLAEPVPGIKAEPDESNARYFHVVIAGESKDSPFEGGTFKRELLLAEEYPMAAPKVRFMTKIYHPNVDKLERIS*DILKDKWSPALQIRTVLLSIQALLNAPNPDD
PLANDVVEQWKTNEAQAIETARAWTRLYAMNSI
>ENSP00000437680.2 gene:ENSG00000205702.11 transcript:ENST00000435101.1
DPAQPPRDLTEAFLAKKEKAKGSPESSFNDENLRIVSVSNRRSTT

到目前为止，我能做的最好的事情就是调用包含以下代码的脚本：

awk -v w=60 -f script.awk fasta_file.fa

#!/usr/bin/env gawk
## Script.awk

/^>/ {
    if (seq != "") print seq; print $1,$4,$5; seq = ""; next
}
    {
        seq = seq $1
        while (length(seq) > w) {
            print substr(seq, 1,w)
            seq = substr(seq, 1+w)
        }
    }
    END { if (seq != "") print seq }

上面代码的问题在于 $1、$4 和 $5 字段是硬编码。
针对类似问题提出了一个优雅的解决方案 Ed Morton，但是，它需要我理解 \s/\S gawk 扩展和 AWK 数组，这是我正在努力做的事情。
任何关于如何使用 AWK（不是 Perl/Python）改进上述代码的想法将不胜感激

【问题讨论】：

不要使用 shebang 在 shell 脚本中调用 awk 解释器，只需调用 awk 'script' "${@:--}"，因为 shebang 没有任何有价值的好处，但会剥夺你分离 shell 擅长的东西的机会从事情来看，awk 在你的 shell 脚本中做得很好。 \s 和 \S 分别表示“任何空格字符”和“任何非空格字符”，就像括号表达式 [[:space:]] 和 [^[:space:]] 中的 POSIX 字符类一样。 awk 数组只是普通的旧关联数组（将字符串映射到值），就像在任何语言中一样。

标签： awk text

【解决方案1】：

这不仅展示了如何使用 awk 做你想做的事，还展示了如何正确构建 shell 脚本以在解析参数后调用 awk（如果你用 shebang 调用 awk 则不能这样做 - 不要那样做）。它使用 GNU awk 进行 gensub() 和第三个参数来匹配()：

$ cat tst.sh
#!/usr/bin/env bash

while getopts ":w:f:" opt; do
    case "$opt" in
        w)  wid=${OPTARG}
            ;;
        f)  flds=${OPTARG}
            ;;
        *)  printf 'bad argument "%s"\n' "$opt" >&2
            exit 1
            ;;
    esac
done
shift "$((OPTIND-1))"

awk -v wid="$wid" -v flds="$flds" '
BEGIN {
    wid=(wid ? wid : 120)
    flds=(flds ? flds : "protein gene transcript")
    numTags = split(flds,tags)
}
sub(/^>/,"") {
    if (NR > 1) {
        prt()
    }
    match($0,/(description:.*\S)\s+\[([^]]+)/,a)
    $0 = substr($0,1,RSTART-1)
    f["description"] = a[1]
    f["predictable"] = a[2]
    f["protein"] = $1
    f["peptide"] = $2
    for (i=3; i<=NF; i++) {
        tag = gensub(/:.*/,"",1,$i)
        f[tag] = $i
    }
    next
}
{ f["sequence"] = f["sequence"] $0 }
END { prt() }

function prt(   tagNr, tag) {
    printf ">"
    for (tagNr=1; tagNr<=numTags; tagNr++) {
        tag = tags[tagNr]
        printf "%s%s", f[tag], (tagNr<numTags ? OFS : ORS)
    }
    print gensub(".{"wid"}","&"RS,"g",f["sequence"])
    delete f
}
' "${@:--}"

$ ./tst.sh file
>ENSP00000441696.1 gene:ENSG00000211788.2 transcript:ENST00000390436.2
MTSIRAVFIFLWLQLDLVNGENVEQHPSTLSVQEGDSAVIKCTYSDSASNYFPWYKQELGKGPQLIIDIRSNVGEKKDQRIAVTLNKTAKHFSLHITETQPEDSAVYFCAAS
>ENSP00000488314.1 gene:ENSG00000276380.2 transcript:ENST00000618570.1
MAELPHRIIKETQRLLAEPVPGIKAEPDESNARYFHVVIAGESKDSPFEGGTFKRELLLAEEYPMAAPKVRFMTKIYHPNVDKLERIS*DILKDKWSPALQIRTVLLSIQALLNAPNPDD
PLANDVVEQWKTNEAQAIETARAWTRLYAMNSI
>ENSP00000437680.2 gene:ENSG00000205702.11 transcript:ENST00000435101.1
DPAQPPRDLTEAFLAKKEKAKGSPESSFNDENLRIVSVSNRRSTT

$ ./tst.sh -w 60 -f 'gene_symbol chromosome' file
>gene_symbol:TRAV13-1 chromosome:GRCh38:14:21868839:21869365:1
MTSIRAVFIFLWLQLDLVNGENVEQHPSTLSVQEGDSAVIKCTYSDSASNYFPWYKQELG
KGPQLIIDIRSNVGEKKDQRIAVTLNKTAKHFSLHITETQPEDSAVYFCAAS
>gene_symbol:UBE2NL chromosome:GRCh38:X:143884071:143885255:1
MAELPHRIIKETQRLLAEPVPGIKAEPDESNARYFHVVIAGESKDSPFEGGTFKRELLLA
EEYPMAAPKVRFMTKIYHPNVDKLERIS*DILKDKWSPALQIRTVLLSIQALLNAPNPDD
PLANDVVEQWKTNEAQAIETARAWTRLYAMNSI
>gene_symbol:CYP2D7 chromosome:GRCh38:22:42140203:42141924:-1
DPAQPPRDLTEAFLAKKEKAKGSPESSFNDENLRIVSVSNRRSTT

$ ./tst.sh -w 10000 -f 'description' file
>description:T cell receptor alpha variable 13-1
MTSIRAVFIFLWLQLDLVNGENVEQHPSTLSVQEGDSAVIKCTYSDSASNYFPWYKQELGKGPQLIIDIRSNVGEKKDQRIAVTLNKTAKHFSLHITETQPEDSAVYFCAAS
>description:ubiquitin conjugating enzyme E2 N like (gene/pseudogene)
MAELPHRIIKETQRLLAEPVPGIKAEPDESNARYFHVVIAGESKDSPFEGGTFKRELLLAEEYPMAAPKVRFMTKIYHPNVDKLERIS*DILKDKWSPALQIRTVLLSIQALLNAPNPDDPLANDVVEQWKTNEAQAIETARAWTRLYAMNSI
>description:cytochrome P450 family 2 subfamily D member 7 (gene/pseudogene)
DPAQPPRDLTEAFLAKKEKAKGSPESSFNDENLRIVSVSNRRSTT

$ ./tst.sh -w 10000 -f 'predictable' file
>Source:HGNC Symbol;Acc:HGNC:12108
MTSIRAVFIFLWLQLDLVNGENVEQHPSTLSVQEGDSAVIKCTYSDSASNYFPWYKQELGKGPQLIIDIRSNVGEKKDQRIAVTLNKTAKHFSLHITETQPEDSAVYFCAAS
>Source:HGNC Symbol;Acc:HGNC:31710
MAELPHRIIKETQRLLAEPVPGIKAEPDESNARYFHVVIAGESKDSPFEGGTFKRELLLAEEYPMAAPKVRFMTKIYHPNVDKLERIS*DILKDKWSPALQIRTVLLSIQALLNAPNPDDPLANDVVEQWKTNEAQAIETARAWTRLYAMNSI
>Source:HGNC Symbol;Acc:HGNC:2624
DPAQPPRDLTEAFLAKKEKAKGSPESSFNDENLRIVSVSNRRSTT

【讨论】：