基因组测序模拟

杨利霞

基因组测序模拟
9 y- g- ?; Y) g/ [: U; a基因组测序模拟

  T3 n! v+ \3 j  F6 J+ Y一、摘要

通过熟悉已有的基因组测序模拟和评估程序，加深全基因组鸟枪法测序原理的理解，并且能够编写程序模拟全基因组鸟枪法测序，理解覆盖度、测序深度、拷贝数等概念，设置测序相关参数，生成单端/双端测序结果文件
' C& f  q. W* P
二、材料和方法

# v( s0 c0 B0 P& W% b5 D1、硬件平台

$ B6 j7 |4 @8 `- F) F处理器:Intel（R） Core（TM）i7-4710MQ CPU @ 2.50GHz
# o, x, |; B( q; K- q+ p% q7 I安装内存(RAM)：16.0GB

8 N9 @( k1 L7 R# ?5 q2 D; O7 V2、系统平台
! m5 |# i$ Z  K2 p; W! F: [Windows 8.1，Ubuntu
) ]+ L9 e4 Q$ L/ O
5 O1 P5 j2 a5 y* x/ }. U  F3、软件平台
: n& j- N: B5 q
5 s6 i/ X9 v9 G4 Z0 F6 Bart_454
2 Y* v0 ]& \: [  t( K; GGenomeABC http://crdd.osdd.net/raghava/genomeabc/
' [; I7 |) d# m# J  U9 c. j# o+ v2 \4 WPython3.5
Biopython
4、数据库资源

NCBI数据库：https://www.ncbi.nlm.nih.gov/

  e! P) p2 y- P5、研究对象
6 F3 v( S" H8 o( H; H' b9 t0 T
酵母基因组Saccharomyces cerevisiae S288c (assembly R64)
  R7 z' B: ]% M2 ~9 |ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_genomic.fna.gz
& S+ z: e) ]8 O  H5 L
" o  g; k* g5 Z" {6、方法

6 N3 F* [; B# cart_454的使用
5 T9 Z$ x+ x' T首先至art系列软件的官网，下载软件，在ubuntu系统安装，然后阅读相关参数设置的帮助文档，运行程序。
) |. _" p1 B2 i. ^1 }. d4 eGenomeABC
# A1 m$ p2 n, @- F进入GenomeABC(http://crdd.osdd.net/raghava/genomeabc/)，输入参数，获得模拟测序结果。
; G1 d; V" d1 h1 C& o7 }5 V$ r编程模拟测序
下载安装python，并且安装biopython扩展模块，编写程序，模拟单端/双端测序。
. m  R3 q+ ~6 a三、结果

1、art_454的运行结果
7 l) O5 S1 o, g% C0 Z
无参数art_454运行，阅读帮助文档

图表 1无参数art_454运行
6 J; R" ]5 l6 Q/ t2 G- C对酵母基因组进行基因组单端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20.
下图为模拟单端测序，程序运行过程及结果
+ \) t8 x2 I4 z( H
4 n/ H' L' Q: ?, e图表 2 art454单端测序

图表 3 art454单端模拟结果
双端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20；MEAN_FRAG_LEN设为1500，即平均片段长度为1500；STD_DEV设为20，即长度的标准差为20
7 @. ?; \- Z: X# \下图为模拟双端测序，程序运行过程及结果

图表 4 art454双端测序
6 f  p" ~: w" z$ ~1 w: }7 `" Z
图表 5 art454双端模拟结果
; w8 Q7 M- v" Q, q2、GenomeABC
下图为设置参数页面
7 n% W1 {3 |, p. v* G. ^" c9 e2 |
- N$ |; M' A- z6 J8 Y4 `3 J" F/ k# a下图为结果下载页面

图表 6 结果下载页面
3、编程模拟测序结果
拷贝数是这里的N值；覆盖度是m，测序深度是宏观的量，在这里与覆盖度意思相同，就是测序仪10X，20X。
单端测序
" y: k# a* ~* x" W3 [+ a
" ~) Q, M; p* x8 P1 z7 l图表 7 程序模拟单端测序
; B& @$ U: I% t( z0 k0 R双端测序
, C* Z* Y  A, k0 z; `% W8 P; U
图表 8 程序模拟双端测序
6 I7 f. h7 X  ]: U& R测序结果
6 O) z! D4 u' Z3 R
图表 9 结果文件
: n9 i, X8 E8 h" X3 O! n+ o
- z( j% ^7 ~  ]3 W: `1 m# O' o因为期望片段长度是600bp，在片段长度区间200-1000bp内，所以大部分的片段都没有删除。
  t* S3 ^2 @# I# E) l7 ?3 w测序结果统计表
$ ?6 A2 G) Y) Y2 p7 L# l0 v
6 J0 f8 k( k1 W+ h- V7 G测序方式        基因组大小(bp)        片段长度区间 (bp)        N值        期望片段长度        克隆保留率        片段数量        Reads长度范围（bp）        Reads总数量        Reads总长度        覆盖度(m值)        理论丢失率(e-m)        覆盖率(1-e-m)
# P6 K$ H& B/ T: _8 Z5 c单端        12157kb        200-1000        10        600        0.95        107378        50-100        101968        7645.541kb        0.62889        0.53318        0.46682
3 g9 k4 h( K4 B7 Z  b; @单端        12157kb        200-1000        20        600        0.95        213722        50-100        202996        15227.882kb        1.25259        0.28576        0.71424
4 e- _" e  ~8 H, Q/ D双端        12157kb        200-1000        10        600        0.95        106704        50-100        202770        15212.662kb        1.25134        0.28612        0.71388
双端        12157kb        200-1000        20        600        0.95        214212        50-100        407186        30534.265kb        2.51164        0.08114        0.91886
四、讨论和结论
/ s* Q: S% v- o
程序运行方法
  `) U4 ]% P6 `; [- A4 @: r
, g$ z* Y. ?7 f1 H在类的构造方法init()中，调整参数。
% Q6 L8 \: d* {Averagefragmentlength为片段平均的长度；
minfragmentlength和maxfragmentlength是保留片段的范围；
7 y, U8 l6 l, gcloneRetainprobability是克隆的保留率；
minreadslength和maxreadslength是测序reads的长度范围

模拟测序的诸多方法都封装成了Sequencing类，只需要创建类，并调用singlereadsequencing()和pairreadsequencing()方法，传入文件名的参数即可。

附录
. w  r( P- w& J( s
from Bio import SeqIO
3 \* w, F0 ^# E8 F# x1 Z/ O8 P. \  _from math import exp
; p, h0 j# U4 l9 Fimport random
2 N! X# p; b5 T# M4 O& C, ^/ ^
class Sequencing:
    # N代表拷贝份数
' j4 O; b- ~! E+ n3 `    def __init__(self)
        self.fragmentList = []
1 C  g# S9 i7 B& h1 }9 }* b, {        self.readsID = 1
        self.readsList = []
        self.averagefragmentlength = 650
% n. u8 n# t6 g+ z. S        self.minfragmentlength = 500
        self.maxfragmentlength = 800
( Y  @* y! F1 d' C/ C        self.cloneRetainprobability = 1
5 b0 S8 g3 t1 y) x! S4 w        self.minreadslength = 50
" d+ Y; r2 t$ a4 x/ y" l! `  u3 L        self.maxreadslength = 150
        self.N = 10
        self.genomeLength = 0
        self.allreadslength = 0
; p: N; Q. W5 K4 A
    # 生成断裂点
' {: E& k! w/ Z! Q    def generatebreakpoint(self, seqlen, averageLength):
        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
        breakpoint.append(seqlen)
7 ]9 b( l; E" `, ?2 U        breakpoint.append(0)
        # 把随机断裂点从小到大排序
        breakpoint.sort()
6 x, F7 w% t8 q# Q- I! v        return breakpoint

* u! ?6 J( o0 g8 \    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
) u5 u: B' u1 ~) L* [! F6 d5 @) ?  w. H! _    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
9 r+ D+ L/ Z$ y' J        for i in range(len(breakpoint) - 1):
            fragment = seq[breakpoint:breakpoint[i + 1]]
- o, k& V4 I7 v5 u            if maxfragmentlength > len(fragment) > minfragmentlength:
                self.fragmentList.append(fragment)
        return self.fragmentList
  x- [& m! |3 D( \
/ q  ^# N) Z+ g) J( Y' @1 A    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
+ R8 \, F2 ~6 V    def clonefragment(self, fragmentList, cloneRetainprobability):
        clonedfragmentList = []
3 z7 w& P1 X% X$ }% t1 G+ d        Lossprobability = [random.random() for _ in range(len(fragmentList))]
        for i in range(len(fragmentList)):
            if Lossprobability <= cloneRetainprobability:
                clonedfragmentList.append(fragmentList)
        return clonedfragmentList

5 ~  Y' I9 L# e    # 模拟单端测序,并修改reads的ID号
    def singleread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
            fragment.id = ""
            fragment.name = ""
            fragment.description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + fragment.description
0 L, k8 ?* `2 x1 d            self.readsID += 1
            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
( \" o, \* B: R% I$ x8 t            self.readsList.append(fragment[:readslength])
2 j$ f0 B9 A) f! D( v9 A( Z4 W( c/ y
5 A- V; b9 C$ D5 ?3 A8 o8 l    def singlereadsequencing(self, genomedata, sequencingResult):
        for seq_record in SeqIO.parse(genomedata, "fasta"):
  {1 E' i/ H. c, f" e, \1 f            seqlen = len(seq_record)
: `' U" |, R6 u; ?* m            self.genomeLength += seqlen
            for i in range(self.N):
; s0 ~7 [5 t+ P/ b( H& c) i                # 生成断裂点
% ], z" B$ @; x                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
8 v6 ^1 }& g+ v; G+ m3 q                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
# m8 Y. v2 e6 ]. i        # 模拟克隆时的随机丢失情况
2 Q4 o# a/ z& h9 S1 b: o        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟单端测序
7 y. y9 J0 E% m' e5 I2 r        self.singleread(clonedfragmentList)
3 P) _( _5 K/ w& p2 Y$ C        SeqIO.write(self.readsList, sequencingResult, "fasta")
- u: G/ |5 W& S
    def pairread(self, clonedfragmentList):
  {7 H) k" ~9 e$ ~( ]        for fragment in clonedfragmentList:
            fragment.id = ""
            fragment.name = ""
! G. X5 Q1 z8 e2 r2 G            description = fragment.description[12:].split(",")[0]
" O2 v/ w6 l7 m4 F; p            fragment.description = str(self.readsID) + "." + description
) u( G9 R" Q. Y5 x7 x" K            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])

            readslength = random.randint(self.minreadslength, self.maxreadslength)
/ `: h. E- V: V) y2 m            self.allreadslength += readslength
' f* P& L( m5 H2 P, k7 F
            fragmentcomplement = fragment.reverse_complement()
% w! f+ q4 x8 U3 q" r) s            fragmentcomplement.id = ""
/ T9 n1 J! S9 ^$ e5 [            fragmentcomplement.name = ""
            fragmentcomplement.description = str(self.readsID) + "." + description
/ F& B! W4 l' B! L+ j" f) r            self.readsList.append(fragmentcomplement[:readslength])

9 p3 I3 l  X' P3 x1 x            self.readsID += 1

, q0 S/ B' _  P6 E! H    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
% j; Q( p% C" P9 F        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
0 e& c7 {2 ]$ h  F2 s& J            self.genomeLength += seqlen
            for i in range(self.N):
                # 生成断裂点
' }2 W3 ~; B8 y                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
' I6 N3 s3 Q, K$ n  l6 {                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
! h1 K8 r( l( W) _& V+ e7 S: H        # 模拟克隆时的随机丢失情况
, U4 u# f9 L  ~        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟双端测序
        self.pairread(clonedfragmentList)
- v& e' W) D' W        readsList_1 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 0]
* P7 @; l! n" E1 X$ B" P: t6 ^. R        readsList_2 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 1]
        SeqIO.write(readsList_1, sequencingResult_1, "fasta")
$ i9 q" O% U4 q! g$ o( j9 q. C        SeqIO.write(readsList_2, sequencingResult_2, "fasta")

/ P0 C) K7 G: L/ H3 ?    def resultsummary(self):
        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
        print("N值:" + str(self.N))
+ a9 ~5 a) d" _6 D; c9 a        print("期望片段长度:" + str(self.averagefragmentlength))
) X+ _7 J4 Q% F  p* ?% O: S8 r- ^% |        print("克隆保留率:" + str(self.cloneRetainprobability))
        print("片段数量:" + str(len(self.fragmentList)))
/ v) ~4 T* p- u) ?+ X( y1 h. V        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
# k% I* x# ~  J5 N/ G2 m) C        print("reads总数量:" + str(len(self.readsList)))
        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
        m = self.allreadslength / self.genomeLength
        print("覆盖度(m值):" + str(round(m, 5)))
        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
1 h) w; M9 @# z3 ?3 y0 `! ]0 Z        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
- K5 j* X) [  d4 r0 H1 F; s# -------------------------------------------主程序-------------------------------------------
# 模拟单端测序
sequencingObj = Sequencing()
, F& w6 `  [# x. W, L* V6 x% TsequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
" k: s, Q9 B1 x( G( |sequencingObj.resultsummary()
6 e6 w* B/ w3 n. f
* s& z2 Q8 H+ a, ~0 X# 模拟双端测序
  _  x# @& ^; M& t& L) usequencingObj = Sequencing()
9 G: a7 x2 a* QsequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
2 s- s7 f+ M  YsequencingObj.resultsummary()
from Bio import SeqIO
from math import exp
import random
4 ?0 c5 S( Q7 X, o, f6 Y
& q/ O1 X* m2 y) _& C/ Lclass Sequencing:
    # N代表拷贝份数
/ ~3 |. i. B2 S% U    def __init__(self):
( z9 ?/ L4 `! }        self.fragmentList = []
        self.readsID = 1
        self.readsList = []
* t& O( B9 x* c4 y* r        self.averagefragmentlength = 650
: y* t3 e  u7 @: l- k$ T& \        self.minfragmentlength = 500
        self.maxfragmentlength = 800
& R. R, T6 H& z8 ~        self.cloneRetainprobability = 1
8 X3 N$ G7 s4 k        self.minreadslength = 50
        self.maxreadslength = 150
2 w: z+ e% k' i  s' U        self.N = 10
, }; M% s' ^% u0 c' @) [$ `        self.genomeLength = 0
        self.allreadslength = 0
0 p/ V7 R" s9 _3 C( `* f% [
    # 生成断裂点
    def generatebreakpoint(self, seqlen, averageLength):
        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
        breakpoint.append(seqlen)
% {6 C" y* K: o" h        breakpoint.append(0)
8 S) @) J( g3 ]4 k7 i# p        # 把随机断裂点从小到大排序
        breakpoint.sort()
: J% T  M( O! L; S9 [7 X        return breakpoint

    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
; w, v6 Q0 A% V5 U        for i in range(len(breakpoint) - 1):
, g: i7 Z  b* Z9 K! O            fragment = seq[breakpoint:breakpoint[i + 1]]
            if maxfragmentlength > len(fragment) > minfragmentlength:
                self.fragmentList.append(fragment)
        return self.fragmentList

; D" W5 n7 }! E/ U" j* @" L    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
    def clonefragment(self, fragmentList, cloneRetainprobability):
        clonedfragmentList = []
: R( [  i' E! j4 ]        Lossprobability = [random.random() for _ in range(len(fragmentList))]
        for i in range(len(fragmentList)):
. Y1 v# [0 u: Q, V. }$ g            if Lossprobability <= cloneRetainprobability:
                clonedfragmentList.append(fragmentList)
        return clonedfragmentList

2 ]; L- ~8 ~; y/ @6 a    # 模拟单端测序,并修改reads的ID号
    def singleread(self, clonedfragmentList):
1 r  q0 E5 |' M) b0 }9 S        for fragment in clonedfragmentList:
            fragment.id = ""
            fragment.name = ""
            fragment.description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + fragment.description
6 \  o+ q1 D" F; j' D: ^0 V            self.readsID += 1
, \1 ]9 m7 Y/ |7 k! J! `            readslength = random.randint(self.minreadslength, self.maxreadslength)
0 Q4 Q8 I5 p7 A            self.allreadslength += readslength
# b) i  n) }" h6 N% S            self.readsList.append(fragment[:readslength])

% }6 H' |4 b1 s5 A. i! z    def singlereadsequencing(self, genomedata, sequencingResult):
1 V" W. W: s0 q3 ?0 k        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
            self.genomeLength += seqlen
6 ?8 U; R/ p, u9 r. x2 e: q            for i in range(self.N):
6 \# W8 d' {: @  D3 r1 }                # 生成断裂点
! G# ^" ]4 O) c! d. U                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
6 y' t8 [* ^, A1 U        # 模拟克隆时的随机丢失情况
3 h& F, ]" h6 N: `- h4 K        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
. h0 g, v- g2 j- b. a- ~        # 模拟单端测序
        self.singleread(clonedfragmentList)
        SeqIO.write(self.readsList, sequencingResult, "fasta")
4 a3 J# H* P2 A! l  k
    def pairread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
            fragment.id = ""
4 H8 T% X( |: B/ K( c            fragment.name = ""
! s2 F/ M- A" w( L2 N- c* T, R9 _) m            description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + description
            readslength = random.randint(self.minreadslength, self.maxreadslength)
6 U+ J4 e( U0 A. B$ u$ ^) }            self.allreadslength += readslength
( Q" n5 }. e% [            self.readsList.append(fragment[:readslength])

            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength

            fragmentcomplement = fragment.reverse_complement()
            fragmentcomplement.id = ""
/ ?# [" }4 P' c& }            fragmentcomplement.name = ""
% z5 D- q" f" C4 k% Y% w: k0 n            fragmentcomplement.description = str(self.readsID) + "." + description
- j9 ^" d9 E6 P. _, `' c# M            self.readsList.append(fragmentcomplement[:readslength])

            self.readsID += 1
: r3 y7 C4 Y0 |' B8 _
( {- u. E5 |5 ]! m# w4 M    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
- f1 t# b3 O7 }4 e5 W& W" C        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
0 i( \8 x* I' \. B6 l9 ^: \- l) p8 E5 X            self.genomeLength += seqlen
- c. x# f) o- S8 H            for i in range(self.N):
# o$ b. R, q8 ?' x7 k& f. t                # 生成断裂点
& D5 O( ~- s' J3 `7 r: C1 |1 D                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
9 [" u* s+ `3 g, r) V3 L                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
1 `. D" s7 V$ E9 ?, X        # 模拟双端测序
3 P7 z1 J3 K. T+ y6 `: k8 |2 B$ \( K8 M        self.pairread(clonedfragmentList)
        readsList_1 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 0]
- V/ A& c8 k( s% W        readsList_2 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 1]
3 N- c0 q4 j5 p0 D) S4 M! u& B2 f        SeqIO.write(readsList_1, sequencingResult_1, "fasta")
: x$ N$ R3 }# ]7 _2 J) C+ l        SeqIO.write(readsList_2, sequencingResult_2, "fasta")
5 D" o0 r( _  P+ ?: O
! C- d+ {1 S* y8 _    def resultsummary(self):
( X; O( N7 Z0 }, _5 \        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
  e0 M/ R( a  _' g4 M# I        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
        print("N值:" + str(self.N))
) ^( A+ I( x; Z' S* _        print("期望片段长度:" + str(self.averagefragmentlength))
4 M: z3 C& ~; y. q        print("克隆保留率:" + str(self.cloneRetainprobability))
        print("片段数量:" + str(len(self.fragmentList)))
        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
8 v0 g/ v  j# [5 x2 i        print("reads总数量:" + str(len(self.readsList)))
        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
        m = self.allreadslength / self.genomeLength
" @% p% q% C8 X8 a7 g+ u        print("覆盖度(m值):" + str(round(m, 5)))
8 R" F8 U# @) R; ]        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
& ], [: y/ ^  E1 Z        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
# -------------------------------------------主程序-------------------------------------------
# 模拟单端测序
sequencingObj = Sequencing()
sequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
sequencingObj.resultsummary()
8 H5 a9 l2 X, r% E. H
# 模拟双端测序
8 j0 W2 e# C" {' D9 ksequencingObj = Sequencing()
sequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
sequencingObj.resultsummary()
7 }9 {. I3 E- `4 [7 v9 ]1 u/ P. n
! A. W% ]& h9 e2 ~( b

8 ]' Q+ S5 {: M, [& c

3477959497

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