基因组回文序列挖掘——我的第一个perl脚本！

        这周病毒所的生物信息学课上老师简单讲了一下perl语言编程。之前主要用R语言编程，有时也用python编程，从来没有用过perl，所以那堂课我很认真地听了一下，并了解了perl编程的语法。巧妙的是这周的基因组学课上老师给我们布置了一个读文献的作业，搜集不同物种基因组的信息。考虑到短时间内读文献并不能搜集到足够的物种基因组信息，我突发奇想，想通过刚刚学的perl语言编写一个脚本挖掘基因组的回文序列。首先理解一下回文序列概念：双链DNA中含有的两个结构相同、方向相反的序列称为回文结构，每条单链以任一方向阅读时都与另一条链以相同方向阅读时的序列是一致的，例如5’GGTACC3′ 3’CCATGG5’。回文结构在某些因素作用下，可形成茎环式的十字形结构发夹结构，这些结构形式都可以互相转变。回文结构普遍存在于细胞的DNA分子中，它们与遗传信息表达的调控和基因转移有关。通常情况下短的回文结构可能是一种特别的信号，如二型限制性内切酶的识别位点。较长的回文结构容易转化成发夹结构。

        理解了回文序列的概念，我们再来了解一下perl编程，和前面我们所发的R语言编程的语法不同，perl编程环境下通常包含3种数据结构，标量$，数组@和哈希%，基础的内容读者可以百度perl编程入门或者看相关数据（羊驼书）了解（在我研一课程结束后我也会发表相关编程语言的基础入门学习教程，因为现在课业和分析任务繁忙，望读者们见谅）。和R语言一样，perl语言在分析生物学数据的时候也需要加载相关第三包，一般Bioperl使用较频繁。安装Linux系统安装Bioperl包通常用以下命令：

perl -MCPAN -e shell    

install Bundle::CPAN

install Module::Build    

o conf prefer_installer MB    

o conf commit　＃升级CPAN    

install Bio::SeqIO

        Windows系统直接从Bioperl官网下载即可，本次分析以流感病毒的pb-1片段为例。

use strict;

use warnings;

use Bio::Seq;

use Bio::SeqIO;

$in = Bio::SeqIO->new(-file => “C:\\Users\\admin\\Desktop\\华大拼接序列\\1-pb1.fasta”, -format => ‘Fasta’);

$totalseq = $in -> next_seq;

$longseq = $totalseq -> seq;

$nameseq = $totalseq -> display_name;

$longseq = uc($longseq)

        对于寻找回文序列大致思路是将这段序列的互补序列找到，然后再比较互补序列与原序列在某些区域上是否存在反向关系。那么要分析回文序列，第一步是找到自己序列的互补序列，根据碱基配对原则，可以编写一个函数来完成互补序列的创造。为了简化函数的编写，要将读取的字符串型的fasta序列转化为数组格式。

undef @longseq;

for($i=0; $i < length($longseq); $i = $i + 1,){

 @longseq[$i] = substr($longseq, $i, 1);

}

        接下来是函数的构造，在perl语言中构建函数为sub name{}形式构建。下面的代码中@_表示传递我们输入的数组变量，返回的字符串型数据$cmcha储存了原始序列的互补序列

sub comp{

 undef $cmcha;

 for($i=0; $i < scalar(@_); $i = $i + 1, ){

if(@_[$i] eq “A”){$cmcha = $cmcha . “T”}

elsif(@_[$i] eq “T”){$cmcha = $cmcha . “A”}

elsif(@_[$i] eq “C”){$cmcha = $cmcha . “G”}

elsif(@_[$i] eq “G”){$cmcha = $cmcha . “C”}

else{$cmcha = $cmcha};

 }; return $cmcha;

}

$longseq2 = comp(@longseq)

         对于一般情况而言，回文序列通常如5’GGTACC3′ 3’CCATGG5’这样，显而易见，他们其中一条序列反向后便与另外一条序列一致。所以，基于此，可以判断自己序列和互补序列取反向是否一致来判断是否存在简单回文序列，对于字符串而言，Perl语言中的reverse函数就可以让字符串反向输出，由于不清楚回文序列的具体长度，我们需要将原始序列从第一个位点开始，一个区域一个区域往后读，并逐渐加长区域范围。用代码展示为如下

undef @doubleseq

for($i=0; $i < length($longseq); $i = $i + 1, ){

 for($j=6; $j < length($longseq); $j = $j + 1, ){

  $seq1 = substr($longseq, $i, $j);

  $seq2 = substr($longseq2, $i, $j);

  $A = $seq1;

  $B = reverse($seq2);

  if(substr($A, 0, length($A)) eq substr($B, 0, length($B))){@doubleseq = (@doubleseq, $seq1)}

  else{@doubleseq = @doubleseq};

 };

}

        上述代码中substr函数为从字符串中提取子集（这里的子集大小与全集一样），eq代表字符串与字符串的对比，可以理解为等号的意思，即如果自己序列的某一区段与互补序列取反向后的对应区段一致，那么就是回文序列。运行此代码得到以下结果

        算法与代码的优化：是否所有回文序列都是与互补序列取反向一致？显然不是，再我们的免疫学课程上就学过一些不符合这个情况的回文序列，在抗体的基因重排过程中，抗体基因片段V基因J基因和D基因片段之间存在非编码的RSS（Recombination Signal Sequence）序列，这些序列中存在7bp和9bp的回文序列，5’CACAGTG3’和ACAAAAACC3’，用于RAG-1和RAG-2识别。以CACAGTG为例，它的互补序列为GTGTCAC,取反后中间的T并不与原序列的A一致，而是互补。所以，上面的算法是存在漏洞的算法。对于这样的情况而言，如何全面挖掘回文序列？通过观察发现，如果一个序列为回文序列，那么以序列的中心为对称轴，在自己序列和互补序列的对称轴两端是反向一致的，用图表示为

        所以，我们只需要判断序列和互补的中心两侧是否呈反向轴对称就可以判断序列是否为回文序列。那么改进优化的算法如下：

undef @oddseq

$longseq2 = comp(@longseq);

for($i=0; $i < length($longseq); $i = $i + 1, ){

 for($j=7; $j < length($longseq) – $i; $j = $j + 1, ){

  $seq1 = substr($longseq, $i, $j);

  $seq2 = substr($longseq2, $i, $j);

  $A = $seq1;

  $B = reverse($seq2);

  $k = 0; 

  $m = (length($seq1) – 1)/2 – 1;

  $l = (length($seq2) – 1)/2 + 1;

  $n = length($seq2);

  if(substr($A, $k, $m) eq substr($B, $k, $m) && substr($A, $l, $n) eq substr($B, $l, $n)){@oddseq = (@oddseq, $A)}

  else{@doubleseq = @doubleseq};

 };

}

        运行该代码得到如下结果

        从图片中可以看出，找到的回文序列也包含了之前找到的回文序列，这样的结果是全面的。

        完整的优化代码如下（其实还可以更优化）

        总体来说，这个脚本从产生想法到完成花费了我一下午加晚上几个小时，由于是我第一次写perl脚本，从我接触perl到写脚本总共花在上面的时间也就23小时左右，我个人认为perl语言入门还是非常简单的。只需要稍微适应语法即可。因为最近任务和课业都很多，已经两天没有更新推文了。希望能赶紧上完课考完试，专心把精力花在课题和编程上。