走过了一路的坑，终于走完了这条路，写下学习总结。
总体想法
要模拟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 测序原理介绍
陈巍学基因
相关名词概念
SNP和Indel
SNP：基因上碱基的变异，原本是A，突变为G、C或者T了。
Indel：基因上碱基的增添、删除，原本是A，然后没有了。
SE和PE
SE：single end 单端读序，短于200的可以采用单端读序，因为可以一次测完。
PE：paired end 双端读序，因为很长，无法从一侧读完，所以采用双端读，延长读长。
coverage
coverage = 总碱基数/DNA长度
这里总碱基数是指的最后读完所有碱基的和，DNA长度是原始DNA未打断前的长度，coverage是人为设定的用于控制最后数据量大小的参数
5' 端和3' 端
DNA的复制是有方向的，必须从5' 端到3' 端进行复制，也就是说读取DNA时也要从5' 端到3' 端。
开始模拟基因测序，双端测序
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
网友评论