美文网首页生信工具生信
art-illumina模拟测序

art-illumina模拟测序

作者: javaLi | 来源:发表于2019-04-13 21:58 被阅读16次

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')

相关文章

  • 序列组装

    1.利用fastqc对模拟测序的序列进行质控分析 1.1 使用art-illumina模拟测序,生成高通量数据(a...

  • art-illumina模拟测序

    1.安装软件及下载基因组数据 1.1 下载art-illumina测序软件 链接 1.2 下载基因组数据 从gen...

  • 54.《Bioinformatics Data Skills》之

    今天学习GenomicRanges包使用最后一小节,GRanges对象覆盖度计算。 测序简单模拟 测序深度计算公式...

  • perl模拟Illumina测序

    走过了一路的坑,终于走完了这条路,写下学习总结。 总体想法 要模拟Illumina测序,首先要了解的就必须是Ill...

  • perl模拟癌症细胞测序

    整体思路 与前面的测序文章相同,忘记了可以向前翻。可以说唯一的区别就是在于变异有很大不同。只要将癌细胞的变异理解,...

  • 微生物软件流程测试之数据模拟

    测序模拟器的系统综述:计算工具,功能和观点原文来自《Systematic review of next-gener...

  • 2020-06-15 `Splatter`是一个用于简单模拟sc

    splatter clp 15 June, 2020 Splatter是一个用于简单模拟单细胞RNA测序数据的R包...

  • biostar handbook|如何模拟测序结果

    为了评价一个工具的性能,通常我们都需要先模拟一批数据。这样相当于有了参考答案,才能检查工具的实际表现情况。因此对于...

  • 高通量测序原理

    测序类型 ROCHE/454 测序 illnumina 测序 Pacbio 测序 nanopore 测序 主流的...

  • RNA-seq测序数据模拟

    在评估不同软件性能的时候,我们会需要模拟一些数据。由于模拟数据当中的情况是已知的,例如差异表达基因的数目。因此,通...

网友评论

    本文标题:art-illumina模拟测序

    本文链接:https://www.haomeiwen.com/subject/yejrwqtx.html