1.安装软件及下载基因组数据
1.1 下载art-illumina测序软件 链接
tar zxvf artsrcmountrainier2016.06.05linux.tgz #解压软件
cd art_src_MountRainier_Linux/
./configure #检查依赖；结果报错，没有libgsl(c语言计算库)，没有报错的话不需要安装gsl,make就行了
#安装gsl
cd ../
wget http://ftp.club.cc.cmu.edu/pub/gnu/gsl/gsl-latest.tar.gz #下载最新版gsl库
tar zxvf gsl-latest.tar.gz
cd gsl-2.5/
./configure
make
sudo make install #安装，需要管理员权限
cd ../
#再次安装
cd art_src_MountRainier_Linux/
./configure
make
sudo make install
#安装没问题了，运行时报错，找不到libgsl.so，因为该库默认放在了/usr/local/lib，所以需指定环境变量
echo 'export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH' >>~/.bashrc #添加到环境变量
source ~/.bashrc
#现在可以运行了
1.2 下载基因组数据
从genome下载基因组序列，解压（酿酒酵母YJM1338）
2.模拟测序
使用art-illumina模拟测序
具体参数如下：
art_illumina是需要运行的程序
-ss illumina不同平台有不同的固定表示，HS20表示HiSeq 2000 (100bp)。
-sam 同时生成sam文件
-i 需要输入的参考基因组
-p 表示输出是paired-end数据，如果-m参数给出的值>=2000，则自动升级成mate-pair
-l 100 表示是100bp的双端数据
-f 表示输出数据的覆盖度，这里是10X
-m 表示paired-end的片段大小
-s 表示-m片段的偏差
-o 需要输出的数据，sequencing1是输出文件的前缀
-ef 加上-ef可以使输出的模拟数据没有错误值，加不加看自己的需求。
#使用nohup挂到后台，可以使用jobs查看
nohup art_illumina -ss HS20 -sam -i GCA_000977205.2_Sc_YJM1338_v1_genomic.fna -p -l 50 -f 1 -m 100 -s 10 -o sequencing1&
nohup art_illumina -ss HS20 -sam -i GCA_000977205.2_Sc_YJM1338_v1_genomic.fna -p -l 50 -f 1 -m 1000 -s 10 -o sequencing2&
nohup art_illumina -ss HS20 -sam -i GCA_000977205.2_Sc_YJM1338_v1_genomic.fna -p -l 100 -f 1 -m 1000 -s 10 -o sequencing3&
nohup art_illumina -ss HS20 -sam -i GCA_000977205.2_Sc_YJM1338_v1_genomic.fna -p -l 100 -f 1 -m 1000 -s 20 -o sequencing4&
nohup art_illumina -ss HS20 -sam -i GCA_000977205.2_Sc_YJM1338_v1_genomic.fna -p -l 100 -f 10 -m 1000 -s 10 -o sequencing5&
生成5个文件，2个fastq文件，2个aln文件，1个sam文件
使用awk统计碱基数
awk '{if (FNR%4==2) print $0}' sequencing11.fq sequencing12.fq | wc -m
awk '{if (FNR%4==2) print $0}' sequencing21.fq sequencing22.fq | wc -m
awk '{if (FNR%4==2) print $0}' sequencing31.fq sequencing32.fq | wc -m
awk '{if (FNR%4==2) print $0}' sequencing41.fq sequencing42.fq | wc -m
awk '{if (FNR%4==2) print $0}' sequencing51.fq sequencing52.fq | wc -m
使用表格统计不同参数下的覆盖度
| 序号 | 基因组大小 | -l(read_length) | -f | -m | -s | base number | 测序深度m(-f的实际值) | 理论丢失率(e-m) | 覆盖率(1-e-m) |
|---|---|---|---|---|---|---|---|---|---|
| 1 | 12372560 | 50 | 1 | 100 | 10 | 12454710 | 1.0066396929980537 | 36.57% | 63.43% |
| 2 | 12372560 | 50 | 1 | 1000 | 10 | 12450528 | 1.0063016869588832 | 36.56% | 63.44% |
| 3 | 12372560 | 100 | 1 | 1000 | 10 | 12319576 | 0.9957 | 36.95% | 63.05% |
| 4 | 12372560 | 100 | 1 | 1000 | 20 | 12322808 | 0.9959788 | 36.94% | 63.06% |
| 5 | 12372560 | 100 | 10 | 1000 | 10 | 123204850 | 9.95791 | 4.735e-5 | 99.995% |
从表中看出reads长度越小，片段越小，测序深度越大，偏差越大，则实际的测序深度越大，覆盖率越高。
3.创建模型
在Genbank的SRA子库中，搜索Saccharomyces cerevisiae YJM1338的测序结果文件。
sra数据
安装sratoolkit(https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/2.9.6/sratoolkit.2.9.6-ubuntu64.tar.gz)
使用sratoolkit的prefetch下载数据,使用vdb-validate进行数据校验，并使用fastq-dump提取文件。
#sratoolkit没有写到环境变量中，直接在软件目录下运行
./prefetch SRR800815 #下载数据，该软件会在/home目录下创建ncbi文件夹，下载的数据在该目录下
./vdb-validate /home/li/ncbi/public/sra/SRR800815.sra #校验数据是否下载完成
./fastq-dump --split-files /home/li/ncbi/public/sra/SRR800815.sra #将数据解压为fastq格式数据
创建文件夹SRA，将解压后的fastq数据保存到里面，使用art软件建模
art_profiler_illumina myHi2000.txt SRA fastq 1
python模拟测序过程（有些过程没有想明白，还存在错误，后面有时间再改）
from Bio import SeqIO
import random
from math import exp
from collections import defaultdict
class simulation:
def __init__(self):
self.fragement=list() #储存片段
self.reads=list() #储存reads
self.readsAllLen=0 #统计所有reads的长度
self.readsLen=100 #默认reads长度为100bp
self.avergeFrangementLen=600 #默认平均片段长为1000
self.minFrangementLen=200 #默认最短片段
self.maxFrangementLen=1000 #默认最长片段
self.std=200
self.dnaLen=0 #统计测序总长
self.coverage=10 #默认测序深度为10X
self.readsID=defaultdict(int) #统计每一个染色体的reads数
def breakSeq(self,seqLen,averageBreak): #计算打碎的片段区间
breakList=[random.randint(0,seqLen) for _ in range(seqLen//averageBreak)]
breakList.append(seqLen)
breakList.append(0)
breakList.sort()
return breakList
def breakGenome(self,breakList,minFrangementLen,maxFragementLen,seq): #去除掉过小或过长的片段，并打碎序列
for i in range(len(breakList)-1):
if (breakList[i+1]-breakList[i])>minFrangementLen and (breakList[i+1]-breakList[i])<maxFragementLen:
self.fragement.append(seq[breakList[i]:breakList[i+1]])
def PCRIncrease(self,probability): #模拟PCR中片段的丢失
PCRFrangement=list()
lossProbability=[random.random() for _ in range(len(self.fragement))]
for i in range(len(self.fragement)):
if lossProbability[i]<probability:
PCRFrangement.append(self.fragement[i])
return PCRFrangement
def singleReads(self,PCRFrangement):
for fragement in PCRFrangement:
fragement_decs=fragement.description[:10]
self.readsID[fragement_decs]+=1
fragement.decsription='@'+fragement_decs+'-'+str(self.readsID[fragement_decs])
self.reads.append(fragement[:self.readsLen])
self.readsAllLen+=self.readsLen
def singleSequencing(self,file,resultFile):
for seq in SeqIO.parse(file,'fasta'):
self.dnaLen+=len(seq)
seq_decs=seq.description[:10]
for i in range(self.coverage):
breakList=self.breakSeq(len(seq),self.avergeFrangementLen)
self.breakGenome(breakList,self.minFrangementLen,self.maxFrangementLen,seq)
PCRFrangement=self.PCRIncrease(0.95)
self.singleReads(PCRFrangement)
SeqIO.write(self.reads,resultFile,'fasta')
print('coverage:',self.readsAllLen/self.dnaLen)
def pairReads(self,PCRFrangement):
for fragement in PCRFrangement:
fragement_decs=fragement.description[:10]
self.readsID[fragement_decs]+=1
fragement.decsription='@'+fragement_decs+'-'+str(self.readsID[fragement_decs])+str(1)
fragement_rev_comp=fragement.reverse_complement()
fragement_rev_comp_desc=fragement_rev_comp.description[:10]
self.readsID[fragement_rev_comp_desc]+=1
fragement_rev_comp.description="@"+fragement_rev_comp_desc+"-"+str(self.readsID[fragement_rev_comp_desc])+str(2)
self.reads.append(fragement[:self.readsLen])
self.reads.append(fragement_rev_comp[:self.readsLen])
self.readsAllLen+=self.readsLen*2
def pairSequencing(self,file,resultFile1,resultFile2):
for seq in SeqIO.parse(file,'fasta'):
self.dnaLen+=len(seq)
seq_decs=seq.description[:10]
for i in range(self.coverage):
breakList=self.breakSeq(len(seq),self.avergeFrangementLen)
self.breakGenome(breakList,self.minFrangementLen,self.maxFrangementLen,seq)
PCRFrangement=self.PCRIncrease(0.95)
self.pairReads(PCRFrangement)
reads1=[self.reads[i] for i in range(len(self.reads)) if i%2==0]
reads2=[self.reads[i] for i in range(len(self.reads)) if i%2==1]
SeqIO.write(reads1,resultFile1,'fasta')
SeqIO.write(reads2,resultFile1,'fasta')
print('coverage:',self.readsAllLen/self.dnaLen)
situlate=simulation()
situlate.singleSequencing('GCA_000977205.2_Sc_YJM1338_v1_genomic.fna','result.fa')
situlatePair=simulation()
situlatePair.pairSequencing('GCA_000977205.2_Sc_YJM1338_v1_genomic.fna','result1.fa','result2.fa')
