走過了一路的坑，終于走完了這條路，寫下學習總結。
總體想法
要模擬Illumina測序，首先要了解的就必須是Illumina的測序原理和流程，不能憑空造輪子。之后才是代碼的編寫，還有對結果的檢測。所以首先了解下IIIumina測序原理。
IIIumina測序
對于IIIumina測序，概括起來可以分為以下幾個簡單的步驟，網(wǎng)上有很多我就稍微提下。
1、打斷
打斷是指將原本一條很長的DNA序列打斷為一些較短的序列，因為IIIumina讀長一般為100~200，不能對很長的DNA序列進行測序，所以要進行打斷操作
2、PCR
基本上高中的時候都有了解了，就是將一個DNA序列進行大量復制，然后目的是進行大規(guī)模測序，獲得大數(shù)據(jù)，降低誤差等等原因
3、讀取
最后一步就是對DNA進行讀取，原理是用特殊處理的脫氧核糖核酸來進行基因的復制，然后這樣每次就只能增加一個堿基，然后根據(jù)不同堿基測出的顏色不同來進行讀取。
說的比較簡略，了解下大概就行。
再留幾個比較好的資料鏈接。
國外的一篇文章：http://thegenomefactory.blogspot.jp/2013/08/paired-end-read-confusion-library.html
20160405 illumina 測序原理介紹
陳巍學基因
相關名詞概念
SNP和Indel
SNP：基因上堿基的變異，原本是A，突變?yōu)镚、C或者T了。
Indel：基因上堿基的增添、刪除，原本是A，然后沒有了。
SE和PE
SE：single end 單端讀序，短于200的可以采用單端讀序，因為可以一次測完。
PE：paired end 雙端讀序，因為很長，無法從一側讀完，所以采用雙端讀，延長讀長。
coverage
coverage = 總堿基數(shù)/DNA長度
這里總堿基數(shù)是指的最后讀完所有堿基的和，DNA長度是原始DNA未打斷前的長度，coverage是人為設定的用于控制最后數(shù)據(jù)量大小的參數(shù)
5' 端和3' 端
DNA的復制是有方向的，必須從5' 端到3' 端進行復制，也就是說讀取DNA時也要從5' 端到3' 端。
開始模擬基因測序，雙端測序
1、隨機模擬DNA序列
這個隨便了，找真數(shù)據(jù)也無所謂，只是我這樣做的，好處是可以隨時獲得指定長度的DNA，方便。
這個很簡單就是使用隨機數(shù)隨機選取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的變量用于設定每個位點發(fā)生SNP的概率，根據(jù)概率決定是否發(fā)生SNP，發(fā)生則替換對應位點的堿基。
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，一般認為是根據(jù)正態(tài)分布，所以要有一個分布區(qū)間，可以設定為sd。
但，本人偷了個懶，沒寫那么復雜，直接隨機了。
PCR則感覺在模擬過程中不是很必要，當然也可以做，可以模仿在PCR中可能產(chǎn)生的SNP和Indel情況，在對根據(jù)大數(shù)據(jù)還原原始序列，說著貌似很有道理啊。
# 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'端為末端的鏈，然后必須有一步的操作就是進行基因的互補反向，這個概念一定要清楚。
然后又一個我現(xiàn)在還不是非常理解但是能知道的就是，當打斷的序列大于1000時，會形成環(huán)來進行讀取，這是情況相反，變量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方法，統(tǒng)計長度為k的基因序列在所有讀取序列中的出現(xiàn)次數(shù)。結果與coverage相關，如coverage設定為20，則kmer最后的峰值也應該在16、7左右。這樣才正確。
介紹下kmer
kmer是指長度為k的所有可能的基因在所有讀取文件中出現(xiàn)的次數(shù)，然后再統(tǒng)計出現(xiàn)相同次數(shù)的基因的數(shù)目。
舉例：
- 我有1條reads:
>read1
ATGCAAA
>read2
AGTAGAA
K取6，則得到四個k-mer, 它們各自出現(xiàn)了一次
輸出1：
ATGCAA 1
TGCAAA 1
AGTAGA 1
GTAGAA 1
則輸出的統(tǒng)計文件結果如下：
輸出2：
1 4
具體實現(xiàn)
首先是讀入reads文件，按行讀取，然后利用hash，用hash{key}++，來做，統(tǒng)計所有數(shù)據(jù)，最后得到一個hash，在利用相同的方法得到最后的出現(xiàn)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;
}
驗證結果
做完上述所有步驟后，現(xiàn)在只要打開你寫入的文件，查看峰值出現(xiàn)的位置，你就可以知道最后做的結果如何了。
結束語
坑還是可以填滿的，只要一直努力做下去。繼續(xù)努力。
倉庫地址：https://github.com/MyandMine/simulate_Illumina
