将类似GenBank的多行记录转换为新的文件格式

将类似GenBank的多行记录转换为新的文件格式

关注
发信
关注(28)粉丝(399)

将类似GenBank的多行记录转换为新的文件格式(fasta格式)

本文介绍了将类似GenBank的多行记录转换为新的文件格式(fasta格式)的处理方法，对大家解决问题具有一定的参考价值，需要的朋友们下面随着小编来一起学习吧！

问题描述

我有一个具有以下体系结构的文件:

I have a file with the following architecture:

 source          1..3566367
                 /organism="Laccaria bicolor S238N-H82"
                 /mol_type="genomic DNA"
                 /strain="S238N-H82"
                 /db_xref="taxon:486041"
 gene            <143..>783
                 /locus_tag="LACBIDRAFT_300552"
                 /db_xref="GeneID:6069681"
 mRNA            join(<143..224,274..309,357..470,524..>783)
                 /locus_tag="LACBIDRAFT_300552"
                 /product="helix-turn-helix transcription factor, AraC
                 type"
                 /note="Has EST support"
                 /transcript_id="XM_001873113.1"
                 /db_xref="GeneID:6069681"
 CDS             join(143..224,274..309,357..470,524..783)
                 /locus_tag="LACBIDRAFT_300552"
                 /note="Helix-turn-helix transcription factor; AraC type"
                 /codon_start=1
                 /product="helix-turn-helix transcription factor, AraC
                 type"
                 /protein_id="XP_001873148.1"
                 /db_xref="GeneID:6069681"
                 /translation="MHAKIFLTILFASAVSVYASPQLEERQILSGVVSAITSAGGAVA
                 SGVTSVAGQVTSVAGSIGGDITSEAGQVFQTVTSIGGKAVTIVTSVGGDAITLATSGA
                 GVATSKFGSVYTVATAAAASEASAATGKSSAALPIHGLQSSLIVGLVTVVGSALLGAA
                 ITL"
 gap             1104..3435
                 /estimated_length=2332
 gene            complement(<3738..>4636)
                 /locus_tag="LACBIDRAFT_242762"
                 /db_xref="GeneID:6069481"
 mRNA            complement(join(<3738..3852,3910..4045,4101..4244,
                 4296..4356,4409..4491,4540..4565,4620..>4636))
                 /locus_tag="LACBIDRAFT_242762"
                 /product="predicted protein"
                 /transcript_id="XM_001873722.1"
                 /db_xref="GeneID:6069481"
 CDS             complement(join(3738..3852,3910..4045,4101..4244,
                 4296..4356,4409..4491,4540..4565,4620..>4636))
                 /locus_tag="LACBIDRAFT_242762"
                 /note="Ribulose-5-phosphate 3-epimerase. RPE2. Fourth step
                 of pentose phosphate pathway; Ribulose-5-phosphate
                 3-epimerase. RPE2. TRuncated gene. No EST evidence"
                 /codon_start=1
                 /product="hypothetical protein"
                 /protein_id="XP_001873757.1"
                 /db_xref="InterPro:IPR000056"
                 /db_xref="GeneID:6069481"
                 /translation="LDVMDGHFVPNITMGAPILSCVHKGVPGIFMDCHMMVAKPEQWV
                 DDIADAGGSLYCFHIEATSDPVSLINTIHKRNMKAGVAISPDTPSTAITDEIANAADM
                 LLVMTVYPGRGGQKFIERCVPKVAELRARFPEKDIEVDGGVGPNTIGICADAGCNVIV
                 AGTAIFGSENPMEVIQRLKDTVNAAQAQSGAKY"
 gap             4637..5256
                 /estimated_length=620
 gap             7701..7750
                 /estimated_length=50

在此文件中，我感兴趣的记录是以 CDS 开头的记录:

Within this file, the records I am interested are the ones that start with the word CDS:

 CDS             join(143..224,274..309,357..470,524..783)
                 /locus_tag="LACBIDRAFT_300552"
                 /note="Helix-turn-helix transcription factor; AraC type"
                 /codon_start=1
                 /product="helix-turn-helix transcription factor, AraC
                 type"
                 /protein_id="XP_001873148.1"
                 /db_xref="GeneID:6069681"
                 /translation="MHAKIFLTILFASAVSVYASPQLEERQILSGVVSAITSAGGAVA
                 SGVTSVAGQVTSVAGSIGGDITSEAGQVFQTVTSIGGKAVTIVTSVGGDAITLATSGA
                 GVATSKFGSVYTVATAAAASEASAATGKSSAALPIHGLQSSLIVGLVTVVGSALLGAA
                 ITL"
 CDS             complement(join(3738..3852,3910..4045,4101..4244,
                 4296..4356,4409..4491,4540..4565,4620..>4636))
                 /locus_tag="LACBIDRAFT_242762"
                 /note="Ribulose-5-phosphate 3-epimerase. RPE2. Fourth step
                 of pentose phosphate pathway; Ribulose-5-phosphate
                 3-epimerase. RPE2. TRuncated gene. No EST evidence"
                 /codon_start=1
                 /product="hypothetical protein"
                 /protein_id="XP_001873757.1"
                 /db_xref="InterPro:IPR000056"
                 /db_xref="GeneID:6069481"
                 /translation="LDVMDGHFVPNITMGAPILSCVHKGVPGIFMDCHMMVAKPEQWV
                 DDIADAGGSLYCFHIEATSDPVSLINTIHKRNMKAGVAISPDTPSTAITDEIANAADM
                 LLVMTVYPGRGGQKFIERCVPKVAELRARFPEKDIEVDGGVGPNTIGICADAGCNVIV
                 AGTAIFGSENPMEVIQRLKDTVNAAQAQSGAKY"

从他们那里，我想将这些记录上的信息转换为以下格式:

From them, I would like to convert the information on these records into the following format:

>XP_001873148.1 GeneID:6069681 LACBIDRAFT_300552 helix-turn-helix transcription factor, AraC
MHAKIFLTILFASAVSVYASPQLEERQILSGVVSAITSAGGAVA
SGVTSVAGQVTSVAGSIGGDITSEAGQVFQTVTSIGGKAVTIVT
SVGGDAITLATSGAGVATSKFGSVYTVATAAAASEASAATGKSS
AALPIHGLQSSLIVGLVTVVGSALLGAAITL
>XP_001873757.1 GeneID:6069481 LACBIDRAFT_242762 hypothetical protein
LDVMDGHFVPNITMGAPILSCVHKGVPGIFMDCHMMVAKPEQWV
DDIADAGGSLYCFHIEATSDPVSLINTIHKRNMKAGVAISPDTP
STAITDEIANAADMLLVMTVYPGRGGQKFIERCVPKVAELRARF
PEKDIEVDGGVGPNTIGICADAGCNVIVAGTAIFGSENPMEVIQ
RLKDTVNAAQAQSGAKY

此新文件的记录应以>"开头.符号，并包含以下摘录自的信息:

The records of this new file should start with the ">" symbol and contain the following information extracted from:

XP_001873148.1 来自:/protein_id=XP_001873148.1"

GeneID:6069681 来自:/db_xref ="GeneID:6069681"

LACBIDRAFT_300552 来自:/locus_tag ="LACBIDRAFT_300552"

helix-turn-helix转录因子，AraC 来自:/product =" helix-turn-helix转录因子，AraC

最后，顺序:

MHAKIFLTILFASAVSVYASPQLEERQILSGVVSAITSAGGAVA
SGVTSVAGQVTSVAGSIGGDITSEAGQVFQTVTSIGGKAVTIVT
SVGGDAITLATSGAGVATSKFGSVYTVATAAAASEASAATGKSS
AALPIHGLQSSLIVGLVTVVGSALLGAAITL

我们可以折叠"到任意数字(在此示例中通常为80、44)

Can we 'folded' to any number (usually 80, 44 in this example)

如果有人能指出正确的方向来使用AWK(而不是Perl或Python)来完成此任务，我将不胜感激.我用来执行此任务的Python/Perl工具越来越难以维护/编译.我希望/相信老旧的AWK有潜力执行这些任务.

I would be very grateful if anyone could point me in the right direction to get this done using AWK, not Perl or Python. The Python/Perl tools I used to perform this task are increasingly becoming harder to maintain/compile. I hope/believe that good/old AWK has the potential to perform these tasks.

推荐答案

每当数据中有标记值对时，我发现最好首先构建一个包含该映射的数组( f [] )，然后只需按其标签/名称访问值即可.

Whenever you have tag-value pairs in your data I find it best to first build an array to contain that mapping (f[]) below and then you can just access the values by their tags/names.

使用GNU awk将第三个参数添加到 match()，然后利用 \ s/\ S gawk扩展:

Using GNU awk for the 3rd arg to match() and then taking advantage of the \s/\S gawk extensions:

$ cat tst.awk
BEGIN {
    wid=(wid ? wid : 44)
}
/^ ?\S/ {
    prt()
    key = $1
    sub(/\s*\S+/,"")
}
{
    gsub(/^\s+|\s+$/,"")
    rec = (rec == "" ? "" : rec " ") $0
}
END { prt() }

function prt(   f, tag, val) {
    if ( key == "CDS" ) {
        while ( match(rec,/\/([^=]+)=(\S+|"[^"]+")/,a) ) {
            tag = a[1]
            val = a[2]
            gsub(/^"|"$/,"",val)
            f[tag] = val
            rec = substr(rec,RSTART+RLENGTH)
        }

        gsub(/\s+/,"",f["translation"])
        gsub(".{"wid"}","&"RS,f["translation"])
        sub(RS"$","",f["translation"])

        print ">" f["protein_id"], f["db_xref"], f["locus_tag"], f["product"]
        print f["translation"]
    }
    rec = ""
}

.

$ awk -f tst.awk file
>XP_001873148.1 GeneID:6069681 LACBIDRAFT_300552 helix-turn-helix transcription factor, AraC type
MHAKIFLTILFASAVSVYASPQLEERQILSGVVSAITSAGGAVA
SGVTSVAGQVTSVAGSIGGDITSEAGQVFQTVTSIGGKAVTIVT
SVGGDAITLATSGAGVATSKFGSVYTVATAAAASEASAATGKSS
AALPIHGLQSSLIVGLVTVVGSALLGAAITL
>XP_001873757.1 GeneID:6069481 LACBIDRAFT_242762 hypothetical protein
LDVMDGHFVPNITMGAPILSCVHKGVPGIFMDCHMMVAKPEQWV
DDIADAGGSLYCFHIEATSDPVSLINTIHKRNMKAGVAISPDTP
STAITDEIANAADMLLVMTVYPGRGGQKFIERCVPKVAELRARF
PEKDIEVDGGVGPNTIGICADAGCNVIVAGTAIFGSENPMEVIQ
RLKDTVNAAQAQSGAKY

如果要为 translation 字段使用其他宽度，则可以更改代码或在命令行上指定代码:

If you want to use a different width for the translation field then you can either change the code or specify it on the command line:

$ awk -v wid=80 -f tst.awk file
>XP_001873148.1 GeneID:6069681 LACBIDRAFT_300552 helix-turn-helix transcription factor, AraC type
MHAKIFLTILFASAVSVYASPQLEERQILSGVVSAITSAGGAVASGVTSVAGQVTSVAGSIGGDITSEAGQVFQTVTSIG
GKAVTIVTSVGGDAITLATSGAGVATSKFGSVYTVATAAAASEASAATGKSSAALPIHGLQSSLIVGLVTVVGSALLGAA
ITL
>XP_001873757.1 GeneID:6069481 LACBIDRAFT_242762 hypothetical protein
LDVMDGHFVPNITMGAPILSCVHKGVPGIFMDCHMMVAKPEQWVDDIADAGGSLYCFHIEATSDPVSLINTIHKRNMKAG
VAISPDTPSTAITDEIANAADMLLVMTVYPGRGGQKFIERCVPKVAELRARFPEKDIEVDGGVGPNTIGICADAGCNVIV
AGTAIFGSENPMEVIQRLKDTVNAAQAQSGAKY

这篇关于将类似GenBank的多行记录转换为新的文件格式(fasta格式)的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！

08-20 09:35
Powered by LMLPHP ©2025 将类似GenBank的多行记录转换为新的文件格式 0.052443