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()       #储存readsself.readsAllLen=0     #统计所有reads的长度self.readsLen=100      #默认reads长度为100bpself.avergeFrangementLen=600   #默认平均片段长为1000self.minFrangementLen=200         #默认最短片段self.maxFrangementLen=1000          #默认最长片段self.std=200self.dnaLen=0               #统计测序总长self.coverage=10            #默认测序深度为10Xself.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 breakListdef 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 PCRFrangementdef singleReads(self,PCRFrangement):for fragement in PCRFrangement:fragement_decs=fragement.description[:10]self.readsID[fragement_decs]+=1fragement.decsription='@'+fragement_decs+'-'+str(self.readsID[fragement_decs])self.reads.append(fragement[:self.readsLen])self.readsAllLen+=self.readsLendef 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]+=1fragement.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]+=1fragement_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')

art-illumina模拟测序相关推荐

基因组生物信息学实验（三）：基因组模拟测序（1）
基因组模拟测序(1):主线的内容 step1:方法通过 NCBI 的子库 Genome 获得 YJM1386 菌株的基因组测序数据. 使用 art 系列软件中的 art_illumina 程序,对下 ...
illumina 二代测序原理及过程
● 参考资料: illumina 双端测序二代测序中 barcodes index 的介绍 illumina 测序原理-百度文库 illumina 测序原理-丁香园 DNA 文库构建和 Illumi ...
图形化开放式生信分析系统开发 - 9 Illumina测序仪测序数据自动拆分
前文链接: 图形化开放式生信分析系统开发 - 1 需求分析及技术实现图形化开放式生信分析系统开发 - 2 样本信息处理图形化开放式生信分析系统开发 - 3 生信分析流程的进化图形化开放式生信分析 ...
Illumina测序原理
Illumina测序原理目录 Illumina测序原理 1. 文库制备 2. 成簇反应 3. 测序阶段 4. 数据分析 illumina二代测序平台特点:基于可逆终止的.荧光标记dNTP,实现边合成 ...
sc-RNA seq与Illumina测序
课程笔记 10X genomics 10X genomics技术:文献在10X官网中的publications.对单个细胞做RNA表达情况的定量分析.基于油包水乳浊液酶反应原理的分子生物学分析系统.可 ...
illumina测序的一些注意事项
在二代测序领域,illumina测序平台以其核心的三项技术(分别是碱基末端可逆修饰.桥式PCR和边合成变测序策略),独步天下.充分发挥了其碱基读取精确.通量提高.测序长度也稳步提高至双末端各300bp ...
ART基因序列生成器，究竟是做什么的？
ART是一款比较流行的模拟数据软件.可以模拟生成三大二代测序平台Illumina's Solexa, Roche's 454和Applied Biosystems' SOLiD的single-end, ...
Third-generation sequencing and the future of genomics 第三代测序和基因组学的未来
摘要第三代远程DNA测序和映射技术正在创造高质量基因组测序的复兴.第二代测序只能产生几百对碱基对的短序列,而第三代单分子测序技术可以产生超过10000个碱基对的短序列,或者映射超过100000个碱基 ...
MER：综述高通量测序应用于病原体和害虫诊断
高通量测序应用于病原体和害虫诊断--综述与实用性建议 High‐throughput identification and diagnostics of pathogens and pests: Ov ...

art-illumina模拟测序

1.安装软件及下载基因组数据

1.1 下载art-illumina测序软件链接

1.2 下载基因组数据

2.模拟测序

3.创建模型

python模拟测序过程（有些过程没有想明白，还存在错误，后面有时间再改）

art-illumina模拟测序相关推荐

最新文章

热门文章

art-illumina模拟测序

1.安装软件及下载基因组数据

1.1 下载art-illumina测序软件 链接

1.2 下载基因组数据

2.模拟测序

3.创建模型

python模拟测序过程（有些过程没有想明白，还存在错误，后面有时间再改）

art-illumina模拟测序相关推荐

最新文章

热门文章

1.1 下载art-illumina测序软件链接