perl模拟Illumina测序

2019年6月13日 305次阅读来源: 破冰前行

走过了一路的坑，终于走完了这条路，写下学习总结。

总体想法

要模拟Illumina测序，首先要了解的就必须是Illumina的测序原理和流程，不能凭空造轮子。之后才是代码的编写，还有对结果的检测。所以首先了解下IIIumina测序原理。

IIIumina测序

对于IIIumina测序，概括起来可以分为以下几个简单的步骤，网上有很多我就稍微提下。

1、打断

打断是指将原本一条很长的DNA序列打断为一些较短的序列，因为IIIumina读长一般为100~200，不能对很长的DNA序列进行测序，所以要进行打断操作

2、PCR

基本上高中的时候都有了解了，就是将一个DNA序列进行大量复制，然后目的是进行大规模测序，获得大数据，降低误差等等原因

3、读取

最后一步就是对DNA进行读取，原理是用特殊处理的脱氧核糖核酸来进行基因的复制，然后这样每次就只能增加一个碱基，然后根据不同碱基测出的颜色不同来进行读取。

说的比较简略，了解下大概就行。
再留几个比较好的资料链接。
国外的一篇文章：http://thegenomefactory.blogspot.jp/2013/08/paired-end-read-confusion-library.html
20160405 illumina 测序原理介绍
 陈巍学基因

开始模拟基因测序，双端测序

1、随机模拟DNA序列

这个随便了，找真数据也无所谓，只是我这样做的，好处是可以随时获得指定长度的DNA，方便。
这个很简单就是使用随机数随机选取AGCT，生成伪随机序列。

my @DNA;
my @Nucleotide = ("A", "G", "C", "T");
for(my $i = 0;$i <$len; $i++)
{   
    my $rand = int(rand(4));
    $DNA[$i] = $Nucleotide[$rand];

}

2、在随机的DNA序列中模拟SNP和Indel

对于SNP，在代码中设定一个叫SNP_rate的变量用于设定每个位点发生SNP的概率，根据概率决定是否发生SNP，发生则替换对应位点的碱基。
Indel也是同样，只是将替换改为增加和删除，一般随机添加或删除1~3个碱基（听公司的前辈说的）。

# simulate the SNP and Indel problom
sub SNP_Indel{
    my ($ref, $sRate, $iRate) = @_;
    my $len = length $ref;
    my $sNum = int $sRate * $len;
    my $iNum = int $iRate * $len;
    my @array = ("A", "G", "C", "T");
    for(my $i=0; $i<$sNum; $i++) {
        substr($ref, int rand $sNum, 1) = $array[int rand 4];
    }
    for(my $i=0; $i<$iNum; $i++) {
        if(int rand 2 == 1) {
            substr($ref, int rand $iNum, int rand 3) = "";
        }
        else {
            my $position = int rand $iNum;
            for(my $i=0; $i<rand 3;$i++) {
                substr($ref, $position , 1) = @array[int rand 4];
            }
        }
    }
    return $ref;
}

3、模拟DNA打断

这一步就是将DNA序列进行一个分割操作随机提取一个长度，叫InsertSize，一般指定为500，但实际上不是100%的500，一般认为是根据正态分布，所以要有一个分布区间，可以设定为sd。
但，本人偷了个懒，没写那么复杂，直接随机了。
PCR则感觉在模拟过程中不是很必要，当然也可以做，可以模仿在PCR中可能产生的SNP和Indel情况，在对根据大数据还原原始序列，说着貌似很有道理啊。

# simulate the option of PCR
sub cutSeq{
    my ($ref, $mean, $range) = @_;
    my $len = int($mean + rand 2 * $range - $range);
    my $PCR = substr $ref, int rand(length($ref) -$len), $len;
    return $PCR;
}

4、PCR

（我省略了。。。。）

5、最后，开始读序

读序操作实际上有一点要注意的是，假设变量A读取的是DNA 5’端为始端的链，然后双端读序，变量B要读的就是末端的序列，但因为末端为3′ 端所以不能进行读取，只有去读5’端为末端的链，然后必须有一步的操作就是进行基因的互补反向，这个概念一定要清楚。
然后又一个我现在还不是非常理解但是能知道的就是，当打断的序列大于1000时，会形成环来进行读取，这是情况相反，变量A要进行互补反向操作。

# start read the DNA
sub getReads{
    my ($PCR, $rate, $read_len) = @_;
    my @array = ("A", "G", "C", "T");
    state $n = 0;
    my $len = length $PCR;
    my $read1 = substr($PCR, 0, $read_len);
    my $read2 = substr($PCR, -$read_len, $read_len);
    if($len > 1000){
        $read1 = &revcom($read1);
    }else{
        $read2 = &revcom($read2);
    }
    for(my $i=0; $i<$read_len; $i++) {
        if(rand(1) < $rate)
        {
            substr($read1, $i, 1) = @array[int rand 4];
        }
        if(rand(1)<$rate)
        {
            substr($read2, $i, 1) = @array[int rand 4];
        }
    }
    print READ1 ">reads_$n\_100 1\n$read1\n";
    print READ2 ">reads_$n\_100 2\n$read2\n";
    $n++;
}

你以为结束了吗？没呢

对结果进行检验

检验用的是kmer方法，统计长度为k的基因序列在所有读取序列中的出现次数。结果与coverage相关，如coverage设定为20，则kmer最后的峰值也应该在16、7左右。这样才正确。

介绍下kmer

kmer是指长度为k的所有可能的基因在所有读取文件中出现的次数，然后再统计出现相同次数的基因的数目。
举例：

我有1条reads:
>read1
ATGCAAA
>read2
AGTAGAA
K取6，则得到四个k-mer, 它们各自出现了一次
输出1：
ATGCAA 1
TGCAAA 1
AGTAGA 1
GTAGAA 1
则输出的统计文件结果如下：
输出2：
1 4

具体实现

首先是读入reads文件，按行读取，然后利用hash，用hash{key}++，来做，统计所有数据，最后得到一个hash，在利用相同的方法得到最后的出现x次的有y中kmer这个结果。

# get frequence that length is k
sub k_frequence
{
    my ($seq, $kmer,$modelHash) = @_;
    my $model;
    for(my $j=0; $j < length($seq)-$kmer; $j++)
    {
        $model = substr($seq, $j, $kmer);
        if(!defined($modelHash->{$model})){
            $modelHash->{$model} = 0;
        }
        $modelHash->{$model}++;
    }
}

# get kmer，input values
sub kmer
{
    my @modelArr = @_;
    my %kmers;
    for(my $i=0; $i<@modelArr; $i++)
    {
        if(!defined($kmers{$modelArr[$i]})){
            $kmers{$modelArr[$i]} = 0;
        }
        $kmers{$modelArr[$i]}++;
    }
    return %kmers;
}

验证结果

做完上述所有步骤后，现在只要打开你写入的文件，查看峰值出现的位置，你就可以知道最后做的结果如何了。

结束语

坑还是可以填满的，只要一直努力做下去。继续努力。
仓库地址：https://github.com/MyandMine/simulate_Illumina

    原文作者：破冰前行
    原文地址: https://www.jianshu.com/p/b04a04a8e0d4
    本文转自网络文章，转载此文章仅为分享知识，如有侵权，请联系博主进行删除。

算法网

perl模拟Illumina测序

总体想法

IIIumina测序

相关名词概念

开始模拟基因测序，双端测序

你以为结束了吗？没呢

结束语