Perl語(yǔ)言統(tǒng)計(jì)核酸GC含量程序開(kāi)發(fā)

張婷婷 龐雪原 晏佳蕓

摘 ?要：生命科學(xué)領(lǐng)域中無(wú)論是PCR還是測(cè)序，都需要了解目的核酸的GC含量。由于目前缺乏獨(dú)立統(tǒng)計(jì)核酸GC含量的程序。該文采用Perl語(yǔ)言編程，設(shè)計(jì)處理了核酸序列統(tǒng)計(jì)GC含量的程序。程序采用了BioPerl環(huán)境，利用字符串拆分、遍歷算法和正則表達(dá)式對(duì)核酸序列進(jìn)行GC堿基的挖掘。最后該腳本可以應(yīng)用于Windows和Linux系統(tǒng)。適合用于多環(huán)境。

關(guān)鍵詞：GC含量 ?生命科學(xué) ?Perl語(yǔ)言編程 ?字符串拆分 ?正則表達(dá)式

中圖分類(lèi)號(hào)：TP311.1 ? 文獻(xiàn)標(biāo)識(shí)碼：A ? ? ? ? ? 文章編號(hào)：1672-3791（2019）05（a）-0022-02

在生命科學(xué)領(lǐng)域中無(wú)論是設(shè)計(jì)引物還是分析高通量測(cè)序數(shù)據(jù)質(zhì)量都非常關(guān)注核酸的GC含量[1]。這是因?yàn)樵贒NA雙鏈中G堿基與C堿基是互補(bǔ)配對(duì)的，而且兩個(gè)配對(duì)堿基之間有3個(gè)氫鍵連接（A堿基與T堿基之間是2個(gè)氫鍵連接）。所以GC含量的多少對(duì)于核酸穩(wěn)定性來(lái)說(shuō)是十分重要的一個(gè)因素，GC含量越高核酸的退火溫度越高，核酸也越相對(duì)更穩(wěn)定。在高通量數(shù)據(jù)質(zhì)控中，GC含量也是非常重要的因素。通過(guò)檢測(cè)GC含量可以判斷測(cè)序儀測(cè)序是否測(cè)序足夠隨機(jī)，是否有其他DNA污染。但是，目前沒(méi)有一個(gè)能夠獨(dú)立統(tǒng)計(jì)核酸GC含量的程序。雖然在Primer軟件中有統(tǒng)計(jì)GC含量的插件，但是不可能把所有的序列都放入Primer中進(jìn)行逐一統(tǒng)一序列的GC含量，高通量分析軟件fastqc也有統(tǒng)計(jì)GC含量插件，但是該軟件僅僅識(shí)別fq文件，對(duì)于fa文件則不識(shí)別。絕大多數(shù)時(shí)候我們所用的文件是fa（FASTA）文件。

Perl是一種穩(wěn)定的跨平臺(tái)編程語(yǔ)言，是報(bào)表提取和報(bào)告語(yǔ)言（Practical Extraction and Report Language）的縮寫(xiě)[2]。它適合用于處理字符串類(lèi)型的數(shù)據(jù)。Perl語(yǔ)言是一種開(kāi)源軟件，根據(jù)GUN通用公共許可證（GPL）獲得許可[3]。根據(jù)Perl語(yǔ)言的特點(diǎn)，可以設(shè)計(jì)相關(guān)程序完成對(duì)核酸序列的GC含量統(tǒng)計(jì)。由于DNA序列是由ATCG四種字符組成的字符串，RNA是由AUCG四種字符組成的字符串，所以它們也是一種字符串的數(shù)據(jù)。利用Perl語(yǔ)言設(shè)計(jì)程序統(tǒng)計(jì)GC含量是可行的。

1 ?分析方法

Perl語(yǔ)言中有專(zhuān)門(mén)的叫BioPerl的子程序可以進(jìn)行核酸文件的讀取。Bioperl可以安裝在普通Windows系統(tǒng)上。統(tǒng)計(jì)GC含量的基本原則是C堿基和G堿基總共的數(shù)量占整條核酸序列堿基數(shù)量的百分比。對(duì)于核酸的GC含量分析首先需要對(duì)核酸字符串進(jìn)行分割，拆分成一個(gè)一個(gè)堿基。拆分后的順序也需要與原序列保持一致。然后使用遍歷算法，逐一進(jìn)行正則表達(dá)式匹配，如果匹配上C或者G就記錄[4]，否則不記錄繼續(xù)進(jìn)行分析。具體代碼為：

use Bio：：SeqIO;#加載BioPerl子程序

my $input_mutiseq_file=$ARGV[0];#讀取輸入的核酸序列文件

my $Muti_seq=Bio：：SeqIO->new（#利用BioPerl中的函數(shù)打開(kāi)核酸序列文件

-file=>$input_mutiseq_file，

-format=>'fasta'，）;

my @GC_count;#先定義GC含量的變量，數(shù)據(jù)格式為數(shù)組格式，可儲(chǔ)存多條核酸GC含量

my$j=1;

while （my$seq=$Muti_seq->next_seq（））{#使用遍歷算法拆分讀取堿基

my $id=$seq->id;

my $seq_seq=$seq->seq;

my $seq_length=length$seq_seq;

my @split_seq = split//，$seq_seq;

my $GC_count;

for （my$i=0;$i<@split_seq;$i=$i+1，）{

if（@split_seq[$i]=～"G|C"）{$GC_count=$GC_count+1} #使用正則表達(dá)式匹配C或G

else{$GC_count = $GC_count}};

$GC_count = $GC_count/$seq_length;

@GC_count[$j]=$id."＼t"."GCcount=".$GC_count."＼n";

$j=$j+1;};

print @GC_count;#輸出結(jié)果

以不同HN分析的流感病毒的4號(hào)核酸片段為例，運(yùn)行該程序得到不同流感病毒4號(hào)核酸片段的GC含量結(jié)果（見(jiàn)圖1）。

2 ?結(jié)語(yǔ)

利用Perl語(yǔ)言開(kāi)發(fā)統(tǒng)計(jì)核酸序列的GC含量，對(duì)于研究核酸穩(wěn)定性、測(cè)序的隨機(jī)性和判斷測(cè)序是否有污染具有指導(dǎo)意義。缺陷，由于該腳本的正確運(yùn)行需要計(jì)算機(jī)同時(shí)擁有Perl語(yǔ)言和BioPerl環(huán)境，而整個(gè)Perl語(yǔ)言及其環(huán)境配置相對(duì)繁瑣。

參考文獻(xiàn)

[1] 陳永孜.基因表達(dá)譜芯片及核酸測(cè)序技術(shù)在癌癥研究中的應(yīng)用現(xiàn)狀[EB/OL].[2019-04-04].http：//kns.cnki.net/kcms/detail/23.1513.Q.20190308.1256.002.html.

[2] 劉磊，朱敏.Perl&R在語(yǔ)料庫(kù)語(yǔ)言學(xué)中的應(yīng)用[J].軟件導(dǎo)刊，2018，17（1）：53-55.

[3] 顧武雄.Linux代理程序安裝[J].網(wǎng)絡(luò)安全和信息化，2018（11）：87.

[4] 武明，孫壘，高磊，等.Perl腳本在試驗(yàn)數(shù)據(jù)解析中的應(yīng)用[J].遙測(cè)遙控，2016，37（3）：54-56.