标签: fasta

我有一个数据,总是以下列格式(称为FASTQ)以四块为单位:

@SRR018006.2016 GA2:6:1:20:650 length=36
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNGN
+SRR018006.2016 GA2:6:1:20:650 length=36
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!+!
@SRR018006.19405469 GA2:6:100:1793:611 length=36
ACCCGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
+SRR018006.19405469 GA2:6:100:1793:611 length=36
7);;).;);;/;*.2>/@@7;@77<..;)58)5/>/

是否有一种简单的sed/awk/bash方式将它们转换为这种格式(称为FASTA):

>SRR018006.2016 GA2:6:1:20:650 length=36
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNGN
>SRR018006.19405469 GA2:6:100:1793:611 length=36
ACCCGCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

原则上,我们想要在每个块中提取前两行并替换@为>.

上周我决定尝试Perl6并开始重新实现我的一个程序.我不得不说,Perl6对于对象编程来说非常简单,这在Perl5中对我来说非常痛苦.

我的程序必须读取和存储大文件,例如全基因组(高达3 Gb或更高,参见下面的示例1)或制表数据.

代码的第一个版本是通过逐行迭代("genome.fa".IO.lines)以Perl5方式制作的.对于正确的执行时间来说,它非常缓慢且无法确定.

my class fasta {
  has Str $.file is required;
  has %!seq;

  submethod TWEAK() {
    my $id;
    my $s;

    for $!file.IO.lines -> $line {
      if $line ~~ /^\>/ {
        say $id;
        if $id.defined {
          %!seq{$id} = sequence.new(id => $id, seq => $s);
        }
        my $l = $line;
        $l ~~ s:g/^\>//;
        $id = $l;
        $s = "";
      }
      else {
        $s ~= $line;
      }
    }
    %!seq{$id} = sequence.new(id => $id, seq => $s);
  }
}


sub MAIN()
{ …

我有一个小的DNA序列快速文件,看起来像这样:

>NM_000016 700 200 234
ACATATTGGAGGCCGAAACAATGAGGCGTGATCAACTCAGTATATCAC

>NM_000775 700 124 236
CTAACCTCTCCCAGTGTGGAACCTCTATCTCATGAGAAAGCTGGGATGAG

>NM_003820 700 111 222
ATTTCCTCCTGCTGCCCGGGAGGTAACACCCTGGACCCCTGGAGTCTGCA

问题:

1)如何将此fasta文件读入R作为数据帧,其中每一行是序列记录,第一列是refseqID,第二列是序列.

2)如何在(开始,结束)位置提取子序列？

NM_000016 1  3 #"ACA"
NM_000775 2  6 #"TAACC"
NM_003820 3  5 #"TTC"

我有以下FASTA文件:

>header1
CGCTCTCTCCATCTCTCTACCCTCTCCCTCTCTCTCGGATAGCTAGCTCTTCTTCCTCCT
TCCTCCGTTTGGATCAGACGAGAGGGTATGTAGTGGTGCACCACGAGTTGGTGAAGC
>header2
GGT
>header3
TTATGAT

我想要的输出:

>header1
117
>header2
3
>header3
7
# 3 sequences, total length 127.

这是我的代码:

awk '/^>/ {print; next; } { seqlen = length($0); print seqlen}' file.fa

我用这段代码得到的输出是:

>header1
60
57
>header2
3
>header3
7

我需要一个小修改来处理多个序列行.

我还需要一种方法来获得总序列和总长度.欢迎任何建议......请在bash或awk中.我知道在Perl/BioPerl中很容易做到这一点,实际上,我有一个脚本可以用这些方式来做.

我试图解析一个大的fasta文件,我遇到了内存错误.一些改进数据处理的建议将不胜感激.目前程序正确打印出名称,但部分通过文件我得到一个MemoryError

这是发电机

def readFastaEntry( fp ):
    name = ""
    seq = ""
    for line in fp:
        if line.startswith( ">" ):
            tmp = []
            tmp.append( name )
            tmp.append( seq )
            name = line
            seq = ""
            yield tmp
        else:
            seq = seq.join( line )

这部分工作后,这里是调用者存根

fp = open( sys.argv[1], 'r' )

for seq in readFastaEntry( fp ) :
    print seq[0]

对于那些与fasta格式不相似的人来说,这是一个例子

>1 (PB2)
AATATATTCAATATGGAGAGAATAAAAGAACTAAGAGATCTAATGTCACAGTCTCGCACTCGCGAGATAC
TCACCAAAACCACTGTGGACCACATGGCCATAATCAAAAAGTACACATCAGGAAGGCAAGAGAAGAACCC
TGCACTCAGGATGAAGTGGATGATG
>2 (PB1)
AACCATTTGAATGGATGTCAATCCGACTTTACTTTTCTTGAAAGTTCCAGCGCAAAATGCCATAAGCACC
ACATTTCCCTATACTGGAGACCCTCC

每个条目以">"开头,表示名称等,然后接下来的N行是数据.除了在开头有">"的下一行之外,没有定义的数据结尾.

我遇到的问题描述有点复杂,我会提供更完整的信息.对于不耐烦的人,这是我可以总结的最简单的方法:

抛出换行符时,将文本文件拆分为大小为N(绑定N,例如36)的所有(重叠)子串的最快(最短执行时间)方式是什么.

我正在编写一个模块,用于解析基于FASTA ascii的基因组格式的文件.这些文件包含所谓的'hg18'人类参考基因组,如果您愿意,可以从UCSC基因组浏览器(go slugs!)下载.

正如您将注意到的,基因组文件由chr [1..22] .fa和chr [XY] .fa组成,以及一组未在此模块中使用的其他小文件.

已经存在几个用于解析FASTA文件的模块,例如BioPython的SeqIO.(对不起,我发布了一个链接,但我还没有这么做.)不幸的是,我能找到的每个模块都没有做我想做的具体操作.

我的模块需要将基因组数据(例如,'CAGTACGTCAGACTATACGGAGCTA'可以是一条线)分成每个重叠的N长度子串.让我举一个例子,使用一个非常小的文件(实际的染色体文件长度在355到2000万个字符之间),N = 8

>>>import cStringIO
>>>example_file = cStringIO.StringIO("""\
>header
CAGTcag
TFgcACF
""")
>>>for read in parse(example_file):
...    print read
...
CAGTCAGTF
AGTCAGTFG
GTCAGTFGC
TCAGTFGCA
CAGTFGCAC
AGTFGCACF

我发现的功能从我能想到的方法中获得了绝对最佳的性能:

def parse(file):
  size = 8 # of course in my code this is a function argument
  file.readline() # skip past the header
  buffer = ''
  for line in file:
    buffer += line.rstrip().upper()
    while len(buffer) >= size:
      yield buffer[:size]
      buffer = …

我有一个fasta文件,用换行符分解序列.我想删除换行符.这是我的文件的一个例子:

>accession1
ATGGCCCATG
GGATCCTAGC
>accession2
GATATCCATG
AAACGGCTTA

我想把它转换成这个:

>accession1
ATGGCCCATGGGATCCTAGC
>accession2
GATATCCATGAAACGGCTTA

我在这个网站上找到了一个潜在的解决方案,如下所示:

cat input.fasta | awk '{if (substr($0,1,1)==">"){if (p){print "\n";} print $0} else printf("%s",$0);p++;}END{print "\n"}' > joinedlineoutput.fasta

但是,这会在每个条目之间放置一个额外的换行符,因此文件如下所示:

>accession1
ATGGCCCATGGGATCCTAGC

>accession2
GATATCCATGAAACGGCTTA

我是一个awk noob,但我开始修改命令.我的猜测是if (p){print "\n";}罪魁祸首......可能print "\n"会增加两个换行符.我无法弄清楚如何只添加一个换行符...这可能很简单,但就像我说的那样,我是一个菜鸟.这是我的(不成功)解决方案:

awk '{if (substr($0,1,1)==">"){print "\n"$0} else printf("%s",$0);p++;}END{print "\n"}' input.fasta > joinedoutput.fasta

但是,这会在文件开头添加一个空行,因为它在打印第一个入藏号之前始终打印一个新行:

{empty line} 
>accession1
ATGGCCCATGGGATCCTAGC
>accession2
GATATCCATGAAACGGCTTA

任何人都有解决方案来获取正确格式的文件？谢谢!

我已经尝试了mathematica代码,用于在这个地址中发布DNA序列的混沌游戏:http: //facstaff.unca.edu/mcmcclur/blog/GeneCGR.html

这是这样的:

genome = Import["c:\data\sequence.fasta", "Sequence"];
genome = StringReplace[ToString[genome], {"{" -> "", "}" -> ""}];
chars = StringCases[genome, "G" | "C" | "T" | "A"];
f[x_, "A"] := x/2;
f[x_, "T"] := x/2 + {1/2, 0};
f[x_, "G"] := x/2 + {1/2, 1/2};
f[x_, "C"] := x/2 + {0, 1/2};
pts = FoldList[f, {0.5, 0.5}, chars];
Graphics[{PointSize[Tiny], Point[pts]}]

我所拥有的fasta序列只是像AACCTTTGATCAAA这样的字母序列,要生成的图形如下:

代码适用于小序列,但是当我想要放置一个巨大的序列,例如几乎40Mb的染色体时,该程序需要花费大量时间并且只显示黑色方块,因此无法进行分析.是否有可能改进上述代码,以便显示它的方格会更大？,方式必须只是方形单位.感谢您的帮助

我试图理解FASTA算法在数据库中搜索查询序列的类似序列的基本步骤.这些是算法的步骤:

1. 识别I和J之间的常见k字
2. 使用k字匹配对角线进行评分,确定10个最佳对角线
3. 使用替换分数矩阵重新构造初始区域
4. 使用间隙加入初始区域,惩罚差距
5. 执行动态编程以查找最终对齐

我对使用PAM250评分矩阵的第3步和第4步以及如何"加入使用差距"感到困惑.

有人可以"尽可能具体地"为我解释这两个步骤.谢谢

我想从multifasta文件中提取与单独的ID列表给出的ID匹配的序列.

FASTA文件seq.fasta:

>7P58X:01332:11636
TTCAGCAAGCCGAGTCCTGCGTCGTTACTTCGCTT
CAAGTCCCTGTTCGGGCGCC
>7P58X:01334:11605
TTCAGCAAGCCGAGTCCTGCGTCGAGAGTTCAAGTC
CCTGTTCGGGCGCCACTGCTAG
>7P58X:01334:11613
ACGAGTGCGTCAGACCCTTTTAGTCAGTGTGGAAAC
>7P58X:01334:11635
TTCAGCAAGCCGAGTCCTGCGTCGAGAGATCGCTTT
CAAGTCCCTGTTCGGGCGCCACTGCGGGTCTGTGTC
GAGCG
>7P58X:01336:11621
ACGCTCGACACAGACCTTTAGTCAGTGTGGAAATCT
CTAGCAGTAGAGGAGATCTCCTCGACGCAGGACT

ID文件id.txt:

7P58X:01332:11636
7P58X:01334:11613

我想获取fasta文件,只包含与id.txt文件中的ID匹配的序列:

>7P58X:01332:11636
TTCAGCAAGCCGAGTCCTGCGTCGTTACTTCGCTTT
CAAGTCCCTGTTCGGGCGCC
>7P58X:01334:11613
ACGAGTGCGTCAGACCCTTTTAGTCAGTGTGGAAAC

我真的很喜欢我在这里和这里的答案中找到的awk方法,但是那里给出的代码仍然不能完美地用于我给出的例子.原因如下:

(1)

awk -v seq="7P58X:01332:11636" -v RS='>' '$1 == seq {print RS $0}' seq.fasta

此代码适用于多行序列,但ID必须单独插入代码.

(2)

awk 'NR==FNR{n[">"$0];next} f{print f ORS $0;f=""} $0 in n{f=$0}' id.txt seq.fasta

此代码可以从id.txt文件中获取ID,但只返回多行序列的第一行.

我想好的方法是修改代码中的RS变量(2),但到目前为止我的所有尝试都失败了.请问,有人帮我吗？