基因组测序模拟

杨利霞

基因组测序模拟
基因组测序模拟

一、摘要
" j' }( u5 F3 R  p, u8 `; f
$ k8 F# |0 W! M1 f. [通过熟悉已有的基因组测序模拟和评估程序，加深全基因组鸟枪法测序原理的理解，并且能够编写程序模拟全基因组鸟枪法测序，理解覆盖度、测序深度、拷贝数等概念，设置测序相关参数，生成单端/双端测序结果文件
5 W$ l% N/ l' Z9 ^7 Z, j+ E6 h- U& f
二、材料和方法

1、硬件平台
9 [& O2 t* _+ M# C: W& w7 z# _8 n
处理器:Intel（R） Core（TM）i7-4710MQ CPU @ 2.50GHz
安装内存(RAM)：16.0GB
+ h, l2 z/ S, u; a" g0 Z
3 {7 g+ @, z) }. K% K2、系统平台
: c. k4 y! Z% k: G; I+ hWindows 8.1，Ubuntu
" ?* \5 i( q; y& x
3、软件平台

8 U- z, L+ b1 ?. J& z, ]: S1 ^art_454
GenomeABC http://crdd.osdd.net/raghava/genomeabc/
Python3.5
Biopython
! B, h3 _# D! H; G$ e  d; n$ l& V4、数据库资源

9 m: G: R- c5 |$ O) \3 oNCBI数据库：https://www.ncbi.nlm.nih.gov/

5、研究对象

* d( J5 E/ G$ |0 m' W酵母基因组Saccharomyces cerevisiae S288c (assembly R64)
9 m# E" V1 L' ?& U( j/ fftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_genomic.fna.gz
$ J; k: ?5 ^) Y# o6 Y
6、方法
. C( Q9 I6 U/ J/ o) w) u
5 y2 o' v. o, s+ |# T5 Z. O  Aart_454的使用
! d0 |7 f9 M3 ~  u首先至art系列软件的官网，下载软件，在ubuntu系统安装，然后阅读相关参数设置的帮助文档，运行程序。
GenomeABC
进入GenomeABC(http://crdd.osdd.net/raghava/genomeabc/)，输入参数，获得模拟测序结果。
编程模拟测序
下载安装python，并且安装biopython扩展模块，编写程序，模拟单端/双端测序。
% H# p2 G2 \0 \* H) N- c三、结果

7 @: U# R- k. N( O3 N1、art_454的运行结果
; N; e1 \; |) P; ?, q7 A5 u
& j3 y) C0 \8 h* G: j无参数art_454运行，阅读帮助文档

9 W8 A- K0 w: m% J4 c* U0 l& x图表 1无参数art_454运行
对酵母基因组进行基因组单端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20.
下图为模拟单端测序，程序运行过程及结果
$ s0 S" o! w* ^3 ~7 J
图表 2 art454单端测序

图表 3 art454单端模拟结果
双端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20；MEAN_FRAG_LEN设为1500，即平均片段长度为1500；STD_DEV设为20，即长度的标准差为20
下图为模拟双端测序，程序运行过程及结果

图表 4 art454双端测序
1 f6 h: A2 l& ?
& i* [0 O  y9 \7 s图表 5 art454双端模拟结果
7 i( `4 |% d& k3 V7 k2、GenomeABC
/ F: t, U9 o  }下图为设置参数页面
/ y9 }6 i* L1 B! ~7 E+ s
下图为结果下载页面
8 ?6 h6 J& V2 S  X& ^* O. F0 R; Y
图表 6 结果下载页面
3、编程模拟测序结果
拷贝数是这里的N值；覆盖度是m，测序深度是宏观的量，在这里与覆盖度意思相同，就是测序仪10X，20X。
单端测序

1 Q5 f3 I* O8 ~# y  ^* h& |6 l0 I1 s图表 7 程序模拟单端测序
- F# s& m9 g, V. }& z2 T# o双端测序

图表 8 程序模拟双端测序
测序结果

图表 9 结果文件

因为期望片段长度是600bp，在片段长度区间200-1000bp内，所以大部分的片段都没有删除。
5 J2 V5 T/ y3 V测序结果统计表
9 ]; e4 J* i4 C% R) Y
6 F9 K7 W9 v* s测序方式        基因组大小(bp)        片段长度区间 (bp)        N值        期望片段长度        克隆保留率        片段数量        Reads长度范围（bp）        Reads总数量        Reads总长度        覆盖度(m值)        理论丢失率(e-m)        覆盖率(1-e-m)
单端        12157kb        200-1000        10        600        0.95        107378        50-100        101968        7645.541kb        0.62889        0.53318        0.46682
单端        12157kb        200-1000        20        600        0.95        213722        50-100        202996        15227.882kb        1.25259        0.28576        0.71424
双端        12157kb        200-1000        10        600        0.95        106704        50-100        202770        15212.662kb        1.25134        0.28612        0.71388
7 v  @4 }" D- M& L% P: b双端        12157kb        200-1000        20        600        0.95        214212        50-100        407186        30534.265kb        2.51164        0.08114        0.91886
3 Q) n8 {8 r' m! |四、讨论和结论
+ s* ^6 t& d( i* Z0 x3 |" }2 w. s7 H' @
; Q/ \' b7 ~' ^2 }程序运行方法
: W$ Y$ [( a3 C% o& G0 ]' L
. N: F9 G$ L; |在类的构造方法init()中，调整参数。
Averagefragmentlength为片段平均的长度；
& Z+ h. X( c. v# v; Sminfragmentlength和maxfragmentlength是保留片段的范围；
0 d) j: E  P2 |4 K0 j9 McloneRetainprobability是克隆的保留率；
minreadslength和maxreadslength是测序reads的长度范围

模拟测序的诸多方法都封装成了Sequencing类，只需要创建类，并调用singlereadsequencing()和pairreadsequencing()方法，传入文件名的参数即可。
# U" U; g0 u* f7 h; Q7 c
附录
9 a( Q3 K8 w6 U! Y+ K1 I( T$ U
% X$ M( D0 r9 m+ ]from Bio import SeqIO
from math import exp
import random

class Sequencing:
( Y9 d0 }6 J, `$ I7 I    # N代表拷贝份数
    def __init__(self)
7 `, k/ X, W, V# F) N: g        self.fragmentList = []
* k, c6 u& r, P        self.readsID = 1
2 S; h* T2 E& G9 e        self.readsList = []
7 N* C' ~" Q$ C4 |7 F* D9 ^- {        self.averagefragmentlength = 650
        self.minfragmentlength = 500
: W! b( E, Z! \( D! \: W- |/ Z        self.maxfragmentlength = 800
% Z- G: L, B5 e" H; c  L' X        self.cloneRetainprobability = 1
        self.minreadslength = 50
        self.maxreadslength = 150
; O* u3 s! |* {9 V/ T) l5 ^1 r        self.N = 10
        self.genomeLength = 0
        self.allreadslength = 0
  d  q: J! G5 F' ]( S; q
    # 生成断裂点
    def generatebreakpoint(self, seqlen, averageLength):
4 u% `8 l! E# z0 w* z0 j& }9 g        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
* r$ `* n' |! i6 Y2 F        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
6 y: v( P' t3 h& f- G        breakpoint.append(seqlen)
        breakpoint.append(0)
5 ~1 w3 y4 y9 ^# m, B; r, I3 {        # 把随机断裂点从小到大排序
6 G% d: J# N! W9 V( K        breakpoint.sort()
6 t5 P2 b( e( a        return breakpoint

7 z6 \' J, A/ H    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
3 |! y* I2 |( ^) T; X    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
        for i in range(len(breakpoint) - 1):
            fragment = seq[breakpoint:breakpoint[i + 1]]
            if maxfragmentlength > len(fragment) > minfragmentlength:
                self.fragmentList.append(fragment)
        return self.fragmentList

    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
    def clonefragment(self, fragmentList, cloneRetainprobability):
5 k  v5 k' A" g7 n" H2 \        clonedfragmentList = []
' T6 [2 j% m0 p# e  Y3 ?, K/ ?0 ~        Lossprobability = [random.random() for _ in range(len(fragmentList))]
7 }: B6 w4 O/ s: h- }% V3 f        for i in range(len(fragmentList)):
0 V7 a. i, f( U* P, \& S            if Lossprobability <= cloneRetainprobability:
                clonedfragmentList.append(fragmentList)
1 F" g. s1 d8 V$ {        return clonedfragmentList
! B7 w5 r4 J& H# x8 ^' F
    # 模拟单端测序,并修改reads的ID号
8 [$ P+ m! s" u3 u# T3 [8 B+ O    def singleread(self, clonedfragmentList):
' q$ s5 h3 c( [. c8 O        for fragment in clonedfragmentList:
, q3 B6 ?4 f* _            fragment.id = ""
            fragment.name = ""
  E, ?2 }5 H5 \9 d7 x1 r            fragment.description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + fragment.description
            self.readsID += 1
            readslength = random.randint(self.minreadslength, self.maxreadslength)
' b, s& E0 B# e' t9 M            self.allreadslength += readslength
( v4 ~' L- D3 H; P4 v1 d            self.readsList.append(fragment[:readslength])
  T: _5 a7 j  L7 a7 n" \
( G- ^  m! k3 O8 g' a2 A5 B    def singlereadsequencing(self, genomedata, sequencingResult):
        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
            self.genomeLength += seqlen
            for i in range(self.N):
3 O; o$ i$ X* @                # 生成断裂点
' m3 m' f5 g  ^& ?; A8 |                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
& u! U# t  B2 G6 ~" i0 j                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟单端测序
        self.singleread(clonedfragmentList)
5 ^' G- r& U1 G$ @% u3 `        SeqIO.write(self.readsList, sequencingResult, "fasta")
( h" E% z. y6 u! V% u* j& R. c
3 T" r# f: [4 c. F    def pairread(self, clonedfragmentList):
( ^" Z5 u! G7 o; j& O        for fragment in clonedfragmentList:
            fragment.id = ""
            fragment.name = ""
6 D  K! G* l7 N, J! o& S6 X. T' P            description = fragment.description[12:].split(",")[0]
! p& H! F+ ?7 D8 q, v) H$ ?            fragment.description = str(self.readsID) + "." + description
            readslength = random.randint(self.minreadslength, self.maxreadslength)
: Y5 p9 Q" u4 a. q7 H* R            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])
8 o1 d& Z& w. @+ M; ]0 ?
( |' h, {/ i+ _+ k- x1 E            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength

& w+ j% \( H; `" @: y3 e            fragmentcomplement = fragment.reverse_complement()
1 ]- H' U; y- q% V' Z, E0 C) s            fragmentcomplement.id = ""
( c) N+ y! J6 {4 c6 X  M) F2 x& c            fragmentcomplement.name = ""
            fragmentcomplement.description = str(self.readsID) + "." + description
* B' T6 C; `$ ]' ~. s! o: T5 Z            self.readsList.append(fragmentcomplement[:readslength])

            self.readsID += 1
* e2 }/ [) L: k9 @0 X
    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
/ N6 C7 _* Y+ P' `% l! y+ n. [) J. T" W        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
: ^% e4 P2 j, H3 q+ p' [  D7 U            self.genomeLength += seqlen
            for i in range(self.N):
$ s- P# }6 o4 Y1 M$ |/ \8 n8 ^* k                # 生成断裂点
7 K4 F; H$ p  |, ~. N4 X                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
0 E- P( l, U* I1 c* d. A& z                # 沿断裂点打断基因组
2 Z( N, ]! ~7 {! c: @                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
        # 模拟克隆时的随机丢失情况
, n0 {: H7 L( c0 v. Y        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
) ~- c# u" F- Z        # 模拟双端测序
9 p- _, e1 A4 f9 [) i        self.pairread(clonedfragmentList)
1 ]  y  o3 F' k7 k  u5 U9 I/ I/ D        readsList_1 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 0]
0 z9 {4 g5 @, ], U  @+ L        readsList_2 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 1]
' Z% f1 k) R1 r1 \        SeqIO.write(readsList_1, sequencingResult_1, "fasta")
        SeqIO.write(readsList_2, sequencingResult_2, "fasta")

    def resultsummary(self):
$ M* ]8 k1 v% r/ G) ?        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
+ J4 f' `# P# s9 n. `. v( v8 i. v        print("N值:" + str(self.N))
        print("期望片段长度:" + str(self.averagefragmentlength))
        print("克隆保留率:" + str(self.cloneRetainprobability))
4 {. L" Q  Q6 u' |5 ?! ~$ ]9 L        print("片段数量:" + str(len(self.fragmentList)))
        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
7 s* C5 u. |& @$ h; L! N3 B4 ~# e5 B        print("reads总数量:" + str(len(self.readsList)))
+ y! ?8 F; C* I" N; s- I/ Q        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
        m = self.allreadslength / self.genomeLength
        print("覆盖度(m值):" + str(round(m, 5)))
" b  n+ G) k4 o. \' h0 K# F3 }        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
# P, m/ e& R9 K0 z7 z" B  F; W        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
# -------------------------------------------主程序-------------------------------------------
# 模拟单端测序
/ g4 C) Y& r( s7 XsequencingObj = Sequencing()
sequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
sequencingObj.resultsummary()
4 M0 K0 w; a  A4 v, T
9 q+ {3 `: y$ l# 模拟双端测序
8 h6 a+ U2 z1 I1 v  |. LsequencingObj = Sequencing()
' @0 F2 {; [% c* H  z, E, GsequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
' @% A, s- O, d( lsequencingObj.resultsummary()
from Bio import SeqIO
# l5 z& v. l! z8 _from math import exp
2 Q- L$ u& N# O/ x6 Wimport random
6 {4 Q, f$ ?# P# V6 v$ O- j
class Sequencing:
    # N代表拷贝份数
    def __init__(self):
8 @8 W3 {, |5 I; L6 i. V        self.fragmentList = []
6 ?; Q% w8 p" H        self.readsID = 1
        self.readsList = []
        self.averagefragmentlength = 650
        self.minfragmentlength = 500
        self.maxfragmentlength = 800
        self.cloneRetainprobability = 1
        self.minreadslength = 50
        self.maxreadslength = 150
        self.N = 10
- N6 U: G( ^& S: ~- l        self.genomeLength = 0
        self.allreadslength = 0

1 h& j3 `5 }5 w, }' ]3 q    # 生成断裂点
    def generatebreakpoint(self, seqlen, averageLength):
3 @8 A, i2 P) }) G% o6 h        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
+ N3 ?3 M. H$ m- Z" T# `        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
( u* f' H/ g$ M- h& \        breakpoint.append(seqlen)
% G. q) F6 C% q: J0 p1 Z        breakpoint.append(0)
        # 把随机断裂点从小到大排序
7 n! R) D2 R, l        breakpoint.sort()
9 A  ^+ S6 ?" p+ g2 K        return breakpoint

    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
: V: r1 j& e, f    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
        for i in range(len(breakpoint) - 1):
            fragment = seq[breakpoint:breakpoint[i + 1]]
            if maxfragmentlength > len(fragment) > minfragmentlength:
" n- [! M3 U' D6 ~/ k                self.fragmentList.append(fragment)
$ W* m$ c) M( i' R$ V        return self.fragmentList

5 t* m2 W, B/ V' o( U( Q7 G    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
    def clonefragment(self, fragmentList, cloneRetainprobability):
( e6 |) f" r8 j( \$ p+ @& o  }        clonedfragmentList = []
5 s9 g+ C$ a+ z        Lossprobability = [random.random() for _ in range(len(fragmentList))]
9 {6 M. |; A# n$ z3 G8 t8 H* R        for i in range(len(fragmentList)):
            if Lossprobability <= cloneRetainprobability:
2 v, E. V& Q$ N                clonedfragmentList.append(fragmentList)
" M) x# T. T+ Y0 L        return clonedfragmentList

    # 模拟单端测序,并修改reads的ID号
    def singleread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
            fragment.id = ""
8 T; t& F2 @! c( B3 f9 N& K            fragment.name = ""
            fragment.description = fragment.description[12:].split(",")[0]
( t3 \: B# R( e( z3 s3 Y) Y- z            fragment.description = str(self.readsID) + "." + fragment.description
+ z& o  Q: o6 w/ t' u+ {0 i* M, i0 B            self.readsID += 1
            readslength = random.randint(self.minreadslength, self.maxreadslength)
* D+ J, ^% K; Z+ v            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])
) J& M; u" X9 P* o
6 z$ t) G% v; l; K    def singlereadsequencing(self, genomedata, sequencingResult):
8 ~+ O' Y3 v6 {9 a$ P        for seq_record in SeqIO.parse(genomedata, "fasta"):
6 r$ L4 ]% k9 Z# U. G8 ]$ E            seqlen = len(seq_record)
" i" x# f- Q9 Y' ~            self.genomeLength += seqlen
6 e4 m. `  ]1 d4 V) o            for i in range(self.N):
                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
# ^% F$ C) A7 W$ y) Y! c+ H                # 沿断裂点打断基因组
( `; R8 {/ \% `$ t: M9 Z: Q; g: S                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
; ^% I" K/ u& g7 `        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟单端测序
        self.singleread(clonedfragmentList)
        SeqIO.write(self.readsList, sequencingResult, "fasta")
% N+ ^/ ]0 c4 [9 d1 w) Z2 x! p
    def pairread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
            fragment.id = ""
            fragment.name = ""
            description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + description
            readslength = random.randint(self.minreadslength, self.maxreadslength)
. T$ l/ S: M6 }5 t# k            self.allreadslength += readslength
6 S7 s/ B6 l( r; H' d! D            self.readsList.append(fragment[:readslength])
; v4 K& Z6 B9 n+ W0 w  }
2 W/ j" T9 I/ \7 t+ Y            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
' R. ?/ W: v( ?( N" D* k' K
            fragmentcomplement = fragment.reverse_complement()
            fragmentcomplement.id = ""
3 h+ W, Y& k; r, f; u            fragmentcomplement.name = ""
2 Z  |5 ^4 A5 ~4 w5 d- G4 m' S3 l            fragmentcomplement.description = str(self.readsID) + "." + description
1 f- U9 T2 a, j) o            self.readsList.append(fragmentcomplement[:readslength])

            self.readsID += 1

    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
( {* `2 K6 ~7 ~6 p        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
8 N2 [" Z; [4 s6 V            self.genomeLength += seqlen
2 a1 O+ q7 S/ A% k* t4 t7 J7 k            for i in range(self.N):
- [  v9 h( L# K2 ?2 ?# z, p( v                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
+ H0 k7 n. W8 Y2 N, _) u                # 沿断裂点打断基因组
9 v( m* _; d9 z4 o3 ?1 L                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
        # 模拟克隆时的随机丢失情况
* b7 X# _% n; e: i; A: O        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟双端测序
        self.pairread(clonedfragmentList)
; A% B$ O4 I! p8 W  m- d        readsList_1 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 0]
        readsList_2 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 1]
3 G  M$ W4 a6 k; E9 K; n8 W! T) i4 t7 p        SeqIO.write(readsList_1, sequencingResult_1, "fasta")
2 I/ T8 U1 Y" l        SeqIO.write(readsList_2, sequencingResult_2, "fasta")
+ A! m! _8 S6 j9 m2 h0 A7 H+ O
    def resultsummary(self):
. i3 M6 {; ^2 A$ N        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
        print("N值:" + str(self.N))
        print("期望片段长度:" + str(self.averagefragmentlength))
- O: p" Y; l! r1 ^" r; c! W        print("克隆保留率:" + str(self.cloneRetainprobability))
        print("片段数量:" + str(len(self.fragmentList)))
% \" f! P" ^# l; r( J4 N" i        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
        print("reads总数量:" + str(len(self.readsList)))
        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
# N. O& d9 b3 `7 C, }$ V- Z* l6 @        m = self.allreadslength / self.genomeLength
2 U7 I: ~  z) A) y6 M        print("覆盖度(m值):" + str(round(m, 5)))
        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
- z3 Q2 V2 X. R* s$ {        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
# -------------------------------------------主程序-------------------------------------------
8 }" P3 n+ g) I* W* t0 D5 R# 模拟单端测序
sequencingObj = Sequencing()
3 k- y! b2 G3 x2 g( WsequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
sequencingObj.resultsummary()
& {9 z) C" i' S+ x- |) D
* m' ~: ]- d+ d9 H& D6 x$ @" u1 ^# 模拟双端测序
5 R0 S, K+ T3 X" F3 R9 gsequencingObj = Sequencing()
sequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
# k) t$ J2 q" h' c9 XsequencingObj.resultsummary()
9 Q4 ?  w' U9 @# Y2 F; i, k
1 n. k7 C0 _  i/ ^
3 z  z' X% `8 a
; }: H. Z" A: K+ J

3477959497

不错。。。。。。。。。。。。。。。。。