0

我正在处理 .fasta 文件的标头（该文件在遗传学/生物信息学中普遍用于存储 DNA/RNA 序列数据）。Fasta 文件的标题以 > 符号（提供特定信息）开头，随后是标题描述的下一行的实际序列数据。序列数据无限延伸，直到下一个 \n 之后是下一个标头及其各自的序列。例如：

>scaffold1.1_size947603
ACGCTCGATCGTACCAGACTCAGCATGCATGACTGCATGCATGCATGCATCATCTGACTGATG....
>scaffold2.1_size747567.2.603063_605944
AGCTCTGATCGTCGAAATGCGCGCTCGCTAGCTCGATCGATCGATCGATCGACTCAGACCTCA....

等等...

所以，我对我正在使用的生物体的基因组的 fasta 标头有疑问。不幸的是，解决这个问题所需的 perl 专业知识似乎超出了我目前的技能水平：S 所以我希望这里有人可以告诉我如何做到这一点。

我的基因组由大约 25000 个 fasta 标头及其各自的序列组成，当前状态的标头给我尝试使用的序列比对器带来了很多麻烦，因此我必须大大简化它们。这是我的前几个标题的示例：

>scaffold1.1_size947603
>scaffold10.1_size550551
>scaffold100.1_size305125:1-38034
>scaffold100.1_size305125:38147-38987
>scaffold100.1_size305125:38995-44965
>scaffold100.1_size305125:76102-78738
>scaffold100.1_size305125:84171-87568
>scaffold100.1_size305125:87574-89457
>scaffold100.1_size305125:90495-305068
>scaffold1000.1_size94939

本质上，我想改进这些看起来像这样：

>scaffold1.1a
>scaffold10.1a
>scaffold100.1a
>scaffold100.1b
>scaffold100.1c
>scaffold100.1d
>scaffold100.1e
>scaffold100.1f
>scaffold100.1g
>scaffold1000.1a

或者甚至这个（但这似乎会更复杂）：

>scaffold1.1
>scaffold10.1
>scaffold100.1a
>scaffold100.1b
>scaffold100.1c
>scaffold100.1d
>scaffold100.1e
>scaffold100.1f
>scaffold100.1g
>scaffold1000.1

我在这里所做的是摆脱基因组每个支架的所有大小数据。对于碰巧碎片化的脚手架，我想用 a、b、c、d 等来表示它们。有一些脚手架有超过 26 个片段，所以也许我可以用 x、y、z、A 来表示它们， B，C，D ....等。

我正在考虑通过类似于以下的简单替换 foreach 循环来做到这一点：

#!/usr/bin/perl -w

### Open the files 
$gen = './Hc_genome/haemonchus_V1.fa';
open(FASTAFILE, $gen);
@lines = <FASTAFILE>;
#print @lines; 

###Add an @ symbol to the start of the label
my @refined;
foreach my $lines (@lines){ 
    chomp $lines;
    $lines =~ s/match everything after .1/replace it with a, b, c.. etc/g;
    push @refined, $lines;
}
#print @refined;


###Push the array on to a new fasta file
open FILE3, "> ./Hc_genome/modded_haemonchus_V1.fa" or die "Cannot open output.txt: $!";

foreach (@refined)
{
    print FILE3 "$_\n"; # Print each entry in our array to the file
}
close FILE3;

但我不知道必须在匹配和替换运算符中的 $1 和 \n 之间添加添加的字母标签。本质上是因为我不确定如何通过字母顺序为特定脚手架的每个片段执行此操作（我所能做的就是在每个片段的开头添加一个 a...）

如果您不介意，请告诉我如何实现这一目标！

非常感激！

安德鲁

4

1 回答 1

2

在 Perl 中，增量运算符++对字符串有“神奇”的作用。例如my $s = "a"; $a++递增$a到"b"。这一直持续到"z"，增量将产生"aa"等等。

您的文件的标题似乎已正确排序，因此我们可以遍历每个标题。从标题中，我们提取开始部分（包括 在内的所有内容.1）。如果这个起始部分与前一个标头的起始部分相同，我们增加我们的序列标识符。否则，我们将其设置为"a"：

use strict; use warnings;  # start every script with these

my $index = "a";
my $prev = "";

# iterate over all lines (rather than reading all 25E3 into memory at once)
while (<>) {

  # pass through non-header lines
  unless (/^>/) {
    print;  # comment this line to remove non-header lines
    next;
  }

  s/\.1\K.*//s;  # remove everything after ".1". Implies chomping

  # reset or increment $index
  if ($_ eq $prev) {
    $index++;
  } else {
    $index = "a";
  }

  # update the previous line
  $prev = $_;

  # output new header
  print "$_$index\n";
}

用法：$ perl script.pl <./Hc_genome/haemonchus_V1.fa >./Hc_genome/modded_haemonchus_V1.fa。

编写接受来自 STDIN 的输入并写入 STDOUT 的程序被认为是一种很好的风格，因为这提高了灵活性。与其在 perl 脚本中硬编码路径，不如保持脚本通用，并使用 shell 重定向操作符<来指定输入。这也为您省去了手动打开文件的麻烦。

示例输出：

>scaffold1.1a
>scaffold10.1a
>scaffold100.1a
>scaffold100.1b
>scaffold100.1c
>scaffold100.1d
>scaffold100.1e
>scaffold100.1f
>scaffold100.1g
>scaffold1000.1a
于 2013-07-03T00:45:30.817 回答