art-illumina模拟测序

最新推荐文章于 2022-08-29 22:29:57 发布

li_python

最新推荐文章于 2022-08-29 22:29:57 发布

阅读量1.8k

点赞数

本文链接：https://blog.csdn.net/li_python/article/details/101462498

版权

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