0

我有一个多 fasta 文件，我需要从中提取 100-200 范围内的碱基，包括它们相应的标题。我知道'cut -c 100-200'可以在没有相应标题的情况下做到这一点。有没有办法在 Perl 或 bash 中做到这一点？

示例文件：

8YS68_00009_00025 GAGTTTGATCCTGGCTCAGAGCGAACGCTGGCGGCAGGCTTAACACATGCAAGTCGAGCGGGCGTAGCAATACGTCAGCGGCAGACGGGTGAGTAACGCGTGGGAACATACCTTTTGGTTCGGAACAACACAGGGAAACTTGTGCTAATACCGGATAAGCTACGGGAAGATT 8YS68_00009_00027 GAGTTTGATCATGGCTCAGAGCGAACGCTGGCGGCAGGCCTAACACATGCAAGTCGAGCGCCGTAGCAATACGGAGCGGCAGACGGGTGAGTAACGCGTGGGAACGTACCTTTCGGTTCGGAATAACTCAGGGAAACTTGAGCTAATACCGAATACGTCCGTAAGGAGAAAGATTTATCGCCGAAAGATCGGCCCGCGTAAGATTAGCTAGTTGGTGAGGTAAGGCTCACCAAGCGACGATCGTTAGCTTGTC 8YS68_00012_00035 GAGTTTGATCATGGCTCAGAACGAACGTTGGCGGCGTGGATTAGGCATGCAAGTCGAACGAATCCCATCTGGGTAACTGGGTGGGGGAAGTGGCGAAAGGGGCAGTAATGCGTGGGTAACCTACCTGGGGACCGGGATAGCCTCCTAACGGATGGGTAATACCGGATACGACCTTCGGAGGCATCTCCTGAAGG

所需输出：seq id ------ATCGATCGATCG-----

seq id ------ATCGATCGATCG-----

seq id ------ATCGATCGATCG-----

这意味着，我想准确地提取每个序列的 100-200 之间的碱基，以及它们的标题。如果序列短于 100 bp，则忽略它。

4

4 回答 4

1

使用Bio::SeqIO，以下代码将从 100 提取到 200 并打印标题。

#!/usr/bin/perl 
use strict; 
use warnings;
use Bio::SeqIO;

my $in_file = "fasta_dat.txt"; 

my $in = Bio::SeqIO->new (-file=> $in_file, -format=>'fasta');
my $out = Bio::SeqIO->new( -file   => '>test.fasta',
                           -format => 'fasta');


while(my $seq = $in->next_seq() ) {
    my $subseq = $seq->trunc(100, 200);
    $out->write_seq($subseq);
}

更新：或者只是在这里采用 choroba 的解决方案

于 2013-05-13T15:17:27.380 回答
0

如果您想要的输出是另一个多 fasta 文件，那么您只需要一点awk. 只需子串你想要的。

awk '!/^>/ { print substr($0, 100, 100); next }1' file.fa

最后1返回 true，启用文件中所有行的默认打印。其余的应该是不言自明的。HTH。

一个推测：

awk '/^>/ { h = $0; getline; print h RS substr($0, 100, 100) }' file.fa

或没有getline：

awk '/^>/ { h = $0; next } h { print h RS substr($0, 100, 100); h = "" }' file.fa
于 2013-05-13T11:47:42.843 回答
0

也许您可以使用以下 python 脚本：

    import sys,re
    i,list1 =0,[]
    for line in open(sys.argv[1]):
      if re.match(r'^[>|;]',line):  print line,
      else:
        for x in line:
          if x != "\n": i+=1
          if 100 < i < 200: list1.append(x)
    print "".join(list1)
于 2013-05-13T12:28:47.247 回答
0

在查看了这些建议并为这个问题工作了一段时间后，我在 Perl 中找到了解决方案。这是我编写的在 Perl 中完成这项工作的重要“循环”。

my $seq  = '';
my $head ;

while (my $seq = <IN>) {
if ($seq =~ m/^>/){
    $head = $seq;
    }
    else{
    my $dna .=$seq;
    my $subseq = substr ($seq, 100, 100);
    my $size = length($subseq);
    if ($size > 99){
        print OUT "$head";
        print OUT "$subseq";
        } 
  }

}

感谢大家的帮助和支持。

于 2013-05-14T11:22:34.620 回答