基因组测序模拟

杨利霞

基因组测序模拟
: n* h7 q- }. \1 s- W. w+ A基因组测序模拟

一、摘要
2 j5 a; `4 Y  ]  p% L( M
; r4 o0 i9 d$ S通过熟悉已有的基因组测序模拟和评估程序，加深全基因组鸟枪法测序原理的理解，并且能够编写程序模拟全基因组鸟枪法测序，理解覆盖度、测序深度、拷贝数等概念，设置测序相关参数，生成单端/双端测序结果文件
8 v4 ~! z6 L9 m  T1 X7 y/ D
; q. f3 ^$ u7 s9 G+ O2 U+ I1 Y/ k二、材料和方法
! r. Y3 `& j& B, l+ }
1、硬件平台
  S: i' f% a) p' r9 g4 c
处理器:Intel（R） Core（TM）i7-4710MQ CPU @ 2.50GHz
( t7 n1 s, w6 y6 a: G安装内存(RAM)：16.0GB

, u, t' V7 O+ V7 ^  c( d2、系统平台
+ H* \. y+ j- T! X' JWindows 8.1，Ubuntu

3、软件平台

art_454
GenomeABC http://crdd.osdd.net/raghava/genomeabc/
2 c$ f8 m: C+ M" OPython3.5
8 I  `7 }8 N% ]6 t0 j% v  DBiopython
4、数据库资源

# Z/ Z3 w, G2 E; K4 X. ZNCBI数据库：https://www.ncbi.nlm.nih.gov/

# n8 o9 U3 y" e5 d5、研究对象

- n7 N$ B* {" x( y. F酵母基因组Saccharomyces cerevisiae S288c (assembly R64)
ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_genomic.fna.gz

* E# P7 c! c2 s: _. z* r, n+ ]8 J4 E6、方法
9 q" F3 |. S2 n% M9 j1 \. B/ ]
art_454的使用
首先至art系列软件的官网，下载软件，在ubuntu系统安装，然后阅读相关参数设置的帮助文档，运行程序。
3 C+ M) S! P& ~2 y' XGenomeABC
3 z6 N* ]9 q4 F3 D  d6 ?5 |进入GenomeABC(http://crdd.osdd.net/raghava/genomeabc/)，输入参数，获得模拟测序结果。
编程模拟测序
! x0 {0 k2 l* _; y& I下载安装python，并且安装biopython扩展模块，编写程序，模拟单端/双端测序。
/ ~( X3 m( x* D: F, L三、结果

1、art_454的运行结果

无参数art_454运行，阅读帮助文档

1 I5 q% p2 n8 A/ D/ R7 Q% ?% u图表 1无参数art_454运行
对酵母基因组进行基因组单端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20.
1 p7 F$ T5 E0 n下图为模拟单端测序，程序运行过程及结果
/ E5 Y- w/ e/ F' K% l) G! c" o( j
8 f1 [- R9 {, G图表 2 art454单端测序

: z2 C, R7 w( s7 T& B2 U图表 3 art454单端模拟结果
: H6 x& }9 U- A8 m: w双端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20；MEAN_FRAG_LEN设为1500，即平均片段长度为1500；STD_DEV设为20，即长度的标准差为20
" ?6 V: r* s& C, I4 k下图为模拟双端测序，程序运行过程及结果

图表 4 art454双端测序
& t. G2 V0 Y2 e& T
图表 5 art454双端模拟结果
2、GenomeABC
下图为设置参数页面

2 I4 F5 B1 `  I下图为结果下载页面

" U- y0 }6 f( Q, P, R* I$ A/ v1 M- G. [图表 6 结果下载页面
( s9 s) R1 q3 X5 O# G7 c/ H0 q3、编程模拟测序结果
$ }8 ?2 @3 w+ F! `- ~  o拷贝数是这里的N值；覆盖度是m，测序深度是宏观的量，在这里与覆盖度意思相同，就是测序仪10X，20X。
单端测序

图表 7 程序模拟单端测序
, K! U+ n! E: X5 \双端测序

# p9 K. G( K+ z2 j/ `/ R图表 8 程序模拟双端测序
# y/ g. @: N4 }4 N0 \测序结果

图表 9 结果文件

因为期望片段长度是600bp，在片段长度区间200-1000bp内，所以大部分的片段都没有删除。
测序结果统计表
1 `$ o; j6 l! u" d- G( k0 A/ E
测序方式        基因组大小(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
双端        12157kb        200-1000        20        600        0.95        214212        50-100        407186        30534.265kb        2.51164        0.08114        0.91886
四、讨论和结论

程序运行方法
; d4 [1 ^; L6 Q3 D
- U4 w+ V* y9 ^  H# r: K, H在类的构造方法init()中，调整参数。
Averagefragmentlength为片段平均的长度；
minfragmentlength和maxfragmentlength是保留片段的范围；
cloneRetainprobability是克隆的保留率；
- @& b! I. @4 L, W" O# p8 Aminreadslength和maxreadslength是测序reads的长度范围

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

附录
6 o& |2 M  |% S$ f$ t7 J
3 l0 v, V6 [4 e: m/ wfrom Bio import SeqIO
% v+ n. N7 E. @  M) ofrom math import exp
2 R; r$ a5 Y5 Himport random
) |* |" O' ]4 r% x7 f# a; k
! R2 b) P9 l* _) w/ uclass Sequencing:
) d8 X% c% F# U5 R3 J    # N代表拷贝份数
6 u7 {) D! _1 C8 X* x7 G. N9 q$ u    def __init__(self)
        self.fragmentList = []
4 Y% c% y4 O7 V        self.readsID = 1
        self.readsList = []
9 `5 X# {" D/ r& R1 c        self.averagefragmentlength = 650
( p: ~6 u1 o4 \" }( i3 b/ Q/ E        self.minfragmentlength = 500
        self.maxfragmentlength = 800
        self.cloneRetainprobability = 1
& a4 i8 D' e' m& J0 N4 Z        self.minreadslength = 50
1 z/ R+ [: H/ l- V) ?, k1 J3 B        self.maxreadslength = 150
        self.N = 10
4 [& \+ ^; K; t3 `- H8 N5 n* S        self.genomeLength = 0
) W% F' @% G, d+ h6 |        self.allreadslength = 0
5 {" K9 @4 v% U
    # 生成断裂点
" V/ J' g1 Q9 Y3 P" j7 n6 _. P    def generatebreakpoint(self, seqlen, averageLength):
        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
  _5 \+ ?, U" P6 A1 K" x        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
/ w# u- }( `- ^3 b# p, P, i1 g        breakpoint.append(seqlen)
+ \, a7 T: i6 L, a5 l$ z# X. G$ b        breakpoint.append(0)
        # 把随机断裂点从小到大排序
        breakpoint.sort()
* I" w* r  ^$ E; w- a        return breakpoint

    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
        for i in range(len(breakpoint) - 1):
2 X9 T7 \2 |$ b1 e- d1 k            fragment = seq[breakpoint:breakpoint[i + 1]]
* F1 @+ b5 |1 A* f6 F0 R: L            if maxfragmentlength > len(fragment) > minfragmentlength:
                self.fragmentList.append(fragment)
) e5 Y$ B* _4 f6 |        return self.fragmentList

2 s# e2 ~4 z2 |  u+ v3 r$ k5 y! @    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
    def clonefragment(self, fragmentList, cloneRetainprobability):
5 y8 F# N8 g3 A* N# n; M        clonedfragmentList = []
        Lossprobability = [random.random() for _ in range(len(fragmentList))]
        for i in range(len(fragmentList)):
            if Lossprobability <= cloneRetainprobability:
* e& M* X$ m( Z' M5 b$ t                clonedfragmentList.append(fragmentList)
        return clonedfragmentList
/ T* x" F8 C! e/ h' Z$ j
: @. W+ J1 z) s# n! t    # 模拟单端测序,并修改reads的ID号
    def singleread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
2 O$ g1 \/ \2 o" b6 x            fragment.id = ""
: q1 r: A) u- A# `$ Z6 Y  L            fragment.name = ""
  C/ j  c: q8 T3 U' Y            fragment.description = fragment.description[12:].split(",")[0]
9 R/ d) F; o. C; p            fragment.description = str(self.readsID) + "." + fragment.description
6 K) s$ p" M) V  Q            self.readsID += 1
# c* h3 @( c  L            readslength = random.randint(self.minreadslength, self.maxreadslength)
. M$ e9 L1 _- }* X# A            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])

    def singlereadsequencing(self, genomedata, sequencingResult):
        for seq_record in SeqIO.parse(genomedata, "fasta"):
) G, G$ H/ ]# f# T: |            seqlen = len(seq_record)
            self.genomeLength += seqlen
            for i in range(self.N):
" T1 U0 h5 t2 h  J0 _+ V7 l" j                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
+ v* T3 P+ P+ T' l                # 沿断裂点打断基因组
* Y2 {6 c# a: A3 \( s# w3 |                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
        # 模拟克隆时的随机丢失情况
! J8 ^6 {' h* V% y  i        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
: V) T* a/ n9 p2 Y        # 模拟单端测序
        self.singleread(clonedfragmentList)
        SeqIO.write(self.readsList, sequencingResult, "fasta")

    def pairread(self, clonedfragmentList):
( W" c# S. L. T1 u# ?        for fragment in clonedfragmentList:
            fragment.id = ""
            fragment.name = ""
            description = fragment.description[12:].split(",")[0]
  H% H$ V, t! {6 ^+ T8 A& r            fragment.description = str(self.readsID) + "." + description
$ R2 r, M. k& M3 W            readslength = random.randint(self.minreadslength, self.maxreadslength)
, F# L6 v, r: d) ]6 x3 w8 M            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])
6 v7 P" |# [' l& j, V7 }% _# V
' a2 |. V/ }% H  I            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
5 ?5 }0 l% v6 N4 m! k3 |* V% F
            fragmentcomplement = fragment.reverse_complement()
, P  S  O$ g! Y9 R1 H            fragmentcomplement.id = ""
0 ?; N. ~1 f3 j7 C2 \8 Q            fragmentcomplement.name = ""
            fragmentcomplement.description = str(self.readsID) + "." + description
- c. m. h7 U+ \# D& x- I4 l            self.readsList.append(fragmentcomplement[:readslength])
9 l/ F- a7 S7 D# }9 J9 H
            self.readsID += 1
. U- z# r4 Z8 K; c# w
( W( Q6 ~2 p0 L3 Q2 y3 g- B    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
% N1 ]1 e1 e! Z$ V            self.genomeLength += seqlen
            for i in range(self.N):
6 l8 V) d- q$ X1 b7 f9 F1 Q                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
' h6 [- o7 W: s0 T/ a* ^        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟双端测序
        self.pairread(clonedfragmentList)
1 G$ ]2 ]8 g' ]0 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]
        SeqIO.write(readsList_1, sequencingResult_1, "fasta")
7 R9 M/ v, `& J5 E        SeqIO.write(readsList_2, sequencingResult_2, "fasta")

    def resultsummary(self):
        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
        print("N值:" + str(self.N))
/ _) g* o* c6 w! x/ a* y' Z. @        print("期望片段长度:" + str(self.averagefragmentlength))
        print("克隆保留率:" + str(self.cloneRetainprobability))
        print("片段数量:" + str(len(self.fragmentList)))
        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
8 q: z2 @8 R9 C. M        print("reads总数量:" + str(len(self.readsList)))
        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
+ W& X1 r# o$ ?3 Q7 R8 d        m = self.allreadslength / self.genomeLength
        print("覆盖度(m值):" + str(round(m, 5)))
        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
: {0 v& A- {# D/ u) U# -------------------------------------------主程序-------------------------------------------
+ k: j( ]' _) {0 }# 模拟单端测序
  s. |3 z/ }! ]5 x/ usequencingObj = Sequencing()
3 H5 e, E( q) _+ d4 XsequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
sequencingObj.resultsummary()

# 模拟双端测序
sequencingObj = Sequencing()
sequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
sequencingObj.resultsummary()
( i8 {* a: k4 U" |from Bio import SeqIO
from math import exp
, d$ ~* X  P+ I0 s8 timport random

8 R# H0 ?% d- ^1 Yclass Sequencing:
    # N代表拷贝份数
    def __init__(self):
. T- @3 w# ~+ {' G( J+ A        self.fragmentList = []
) g  C6 J9 c5 _2 u) k$ B7 t& h8 S        self.readsID = 1
        self.readsList = []
        self.averagefragmentlength = 650
        self.minfragmentlength = 500
3 j! }5 y* }. Y  e        self.maxfragmentlength = 800
7 }9 _6 ^5 v* g. s- p1 L        self.cloneRetainprobability = 1
  ]  w1 b4 F) o7 H        self.minreadslength = 50
        self.maxreadslength = 150
        self.N = 10
3 I& e; R% y% E- R7 K        self.genomeLength = 0
        self.allreadslength = 0

    # 生成断裂点
* f4 H, s# ?$ j- u5 c    def generatebreakpoint(self, seqlen, averageLength):
0 Y. b4 R1 r# O( z- I6 P9 h, Z        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
+ k# t7 z! @4 @6 L        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
        breakpoint.append(seqlen)
; }1 T) v9 }* I. g; ]0 r. m( a; g        breakpoint.append(0)
        # 把随机断裂点从小到大排序
3 q! Z4 d% Y- ~$ b4 m9 @        breakpoint.sort()
        return breakpoint
( Y$ n# J4 l2 ^8 M' e; C
' N6 J% i+ t# u    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
" t" Y" ~! S7 ]6 F0 u    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
0 v9 n( S8 C3 t2 ]% ^' J6 l; p        for i in range(len(breakpoint) - 1):
0 D1 Y$ ?$ s. @$ |1 t            fragment = seq[breakpoint:breakpoint[i + 1]]
! d' K6 S8 J/ d( `& A4 r9 d% \! C9 }/ c            if maxfragmentlength > len(fragment) > minfragmentlength:
. r) n* S+ z: t2 E' |* t2 h$ h, a                self.fragmentList.append(fragment)
        return self.fragmentList
3 l9 M8 C8 V6 Z9 h
    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
    def clonefragment(self, fragmentList, cloneRetainprobability):
        clonedfragmentList = []
6 n8 U1 o: S8 Z( p2 o; }6 W2 o        Lossprobability = [random.random() for _ in range(len(fragmentList))]
        for i in range(len(fragmentList)):
$ W4 @. [" S* H+ w% i7 l            if Lossprobability <= cloneRetainprobability:
+ n' [& \% \3 f0 }8 Y  X# v$ w2 z                clonedfragmentList.append(fragmentList)
        return clonedfragmentList
7 @9 t/ g( O5 E: _: c7 `* q2 r/ n5 h
; B$ G0 ]7 L0 ]4 L1 m# o    # 模拟单端测序,并修改reads的ID号
* {, o/ t4 n+ p8 o    def singleread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
( t0 W9 o9 L# q. S3 P. B            fragment.id = ""
, D, b7 \' G1 T) _* Q            fragment.name = ""
            fragment.description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + fragment.description
            self.readsID += 1
3 M1 ^! T1 p+ Y9 W1 Y            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])
2 \7 M/ p4 a7 f5 O6 r
6 B! B8 e' r# _0 o    def singlereadsequencing(self, genomedata, sequencingResult):
        for seq_record in SeqIO.parse(genomedata, "fasta"):
) }& c6 D2 P8 i2 @" ^            seqlen = len(seq_record)
& \* g, L4 s" }% s* L5 B. Y' m            self.genomeLength += seqlen
8 u) P% K0 k1 |. s/ _6 `            for i in range(self.N):
                # 生成断裂点
0 @& n2 c6 r- W$ B: w9 A7 U( M% l; \4 W                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
' h" H7 E/ w0 |                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟单端测序
% P/ c; U3 n: P8 T8 w        self.singleread(clonedfragmentList)
4 m" e) F8 Q" t! w9 k% O        SeqIO.write(self.readsList, sequencingResult, "fasta")
1 z* W' A, N! J) ^  |
2 c6 U1 w" @0 O    def pairread(self, clonedfragmentList):
& }( I1 z! @6 K        for fragment in clonedfragmentList:
            fragment.id = ""
            fragment.name = ""
0 e; k; h3 r* Z5 J/ x, [  M, }            description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + description
9 V" Q  s! Q) S: R            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])
+ a* @& _5 i$ S7 Q. e* H. o
2 E; k4 O  M+ I  k% J$ ]            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength

" H4 v% b0 S2 R0 t9 C2 Z            fragmentcomplement = fragment.reverse_complement()
            fragmentcomplement.id = ""
            fragmentcomplement.name = ""
9 i- C+ I" _( V: K, ]. a+ J  u" z            fragmentcomplement.description = str(self.readsID) + "." + description
            self.readsList.append(fragmentcomplement[:readslength])
- X" t. X7 N' N8 p/ ?+ G
            self.readsID += 1
/ c/ P) Z- P& ?/ D- f
    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
        for seq_record in SeqIO.parse(genomedata, "fasta"):
3 @) H+ m7 X5 M/ D            seqlen = len(seq_record)
            self.genomeLength += seqlen
            for i in range(self.N):
                # 生成断裂点
2 x+ v0 f2 E1 W& J& v$ ?                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
  J8 x" _/ @/ g5 S& Y3 w2 x                # 沿断裂点打断基因组
0 {6 Y" M! Y4 a  b$ q  @                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
, R. Z. Y( \; e1 ~1 i& e5 V$ ~; F+ {        # 模拟克隆时的随机丢失情况
$ E" A, f. T% X: Y        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
) D- e, B9 s! T6 {        # 模拟双端测序
        self.pairread(clonedfragmentList)
7 N( v' A2 }, M3 h6 s) m% P        readsList_1 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 0]
2 \: c% j8 j2 y, C5 V: {4 j! P        readsList_2 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 1]
        SeqIO.write(readsList_1, sequencingResult_1, "fasta")
9 _4 v4 j3 }4 v7 B+ ?( B        SeqIO.write(readsList_2, sequencingResult_2, "fasta")

    def resultsummary(self):
0 R+ Z7 T2 ~3 P3 r' o        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
7 {4 {/ ?6 {- ?/ g: s" u: {6 L        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
) D$ `: `8 C! m' }6 m% R        print("N值:" + str(self.N))
        print("期望片段长度:" + str(self.averagefragmentlength))
0 ~) s1 f0 W4 F. Y  L! O        print("克隆保留率:" + str(self.cloneRetainprobability))
& J3 R" j0 N7 P$ A% U7 a$ _! a        print("片段数量:" + str(len(self.fragmentList)))
        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
* b/ t& A, m( T) d        print("reads总数量:" + str(len(self.readsList)))
6 s% C; D0 f$ |& w! x. i% V& ]        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
        m = self.allreadslength / self.genomeLength
9 C* w4 g0 }( e5 C) o: X        print("覆盖度(m值):" + str(round(m, 5)))
2 K& |& c. E# _! J        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
, N! R) C6 _. F- i9 K# @        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
# -------------------------------------------主程序-------------------------------------------
1 h4 [; ?, U. W( {# 模拟单端测序
sequencingObj = Sequencing()
$ i. e& l2 k% OsequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
$ }' I/ z: [9 _1 h- N4 |sequencingObj.resultsummary()
& f% c- p8 H5 J. T
# 模拟双端测序
* Q3 U1 Y% j4 o! WsequencingObj = Sequencing()
- r4 G2 O% s" V8 IsequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
, `" J7 q# p& f" H- }* XsequencingObj.resultsummary()



$ @( t& U, @& d) s* R
% y3 @. N' U+ i" m+ s

3477959497

不错。。。。。。。。。。。。。。。。。
$ J( C- x; Q  H+ ]! N4 X