基因组测序模拟

杨利霞

基因组测序模拟
基因组测序模拟
) s$ R- S/ ~5 w  O
一、摘要

+ e, C5 {; r9 O; M% C通过熟悉已有的基因组测序模拟和评估程序，加深全基因组鸟枪法测序原理的理解，并且能够编写程序模拟全基因组鸟枪法测序，理解覆盖度、测序深度、拷贝数等概念，设置测序相关参数，生成单端/双端测序结果文件
9 C- d2 G7 K4 N/ a) U
* f; c1 ^: Y! Z* \; s3 O, ~二、材料和方法

/ d: A$ Y4 F1 t6 b* z1、硬件平台

处理器:Intel（R） Core（TM）i7-4710MQ CPU @ 2.50GHz
安装内存(RAM)：16.0GB

2、系统平台
Windows 8.1，Ubuntu
- t7 n% d; M0 D2 w! m
3 ~' I- ^0 I* f$ g3、软件平台
- m# z  E4 F: M! v) Y+ `3 H! {8 f
art_454
  Z$ {- c2 h4 X8 I/ c& P# aGenomeABC http://crdd.osdd.net/raghava/genomeabc/
: [, C; x0 n7 P9 H, ^Python3.5
Biopython
) d, p1 w7 }6 B5 D4、数据库资源

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

5、研究对象
' u' @# p6 ]( z; v% I' S4 B4 Y+ z
; S6 Z# v8 E# B! x( |1 y  Y酵母基因组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

) K2 M# e9 \# O$ Z6、方法
9 Z' K% {# {( o' r8 ]
5 o# k4 M" Y9 D5 y1 gart_454的使用
首先至art系列软件的官网，下载软件，在ubuntu系统安装，然后阅读相关参数设置的帮助文档，运行程序。
GenomeABC
进入GenomeABC(http://crdd.osdd.net/raghava/genomeabc/)，输入参数，获得模拟测序结果。
! m; R' M! x/ y! t  L* D编程模拟测序
# x" f4 D  [: q; ~4 Q" V下载安装python，并且安装biopython扩展模块，编写程序，模拟单端/双端测序。
. h7 d- ]9 x; ?. [( k" p+ {0 ~" @三、结果

) s. k& H( a' ~4 i$ k' T" G1、art_454的运行结果

2 ~$ r4 T) A5 S, s6 `$ Q4 B; M, \5 q6 g无参数art_454运行，阅读帮助文档

& F6 J3 b3 ]% h- O) d/ B图表 1无参数art_454运行
对酵母基因组进行基因组单端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20.
下图为模拟单端测序，程序运行过程及结果
4 k: A3 m( \  {! m: j, I: p: Q
图表 2 art454单端测序

图表 3 art454单端模拟结果
双端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20；MEAN_FRAG_LEN设为1500，即平均片段长度为1500；STD_DEV设为20，即长度的标准差为20
! p2 o9 Y% O$ d3 C! k: x' W2 e下图为模拟双端测序，程序运行过程及结果

, e" j) i7 i) X6 O3 j7 v图表 4 art454双端测序

图表 5 art454双端模拟结果
2、GenomeABC
, K& X' C. U- G$ n* z) H* z下图为设置参数页面

下图为结果下载页面

- l4 v+ V4 g+ W: W6 o图表 6 结果下载页面
* t. u( s2 T' o0 X4 n9 F" e3、编程模拟测序结果
拷贝数是这里的N值；覆盖度是m，测序深度是宏观的量，在这里与覆盖度意思相同，就是测序仪10X，20X。
; Y2 S& v3 c0 b3 y" d1 x5 [* f单端测序
& Q" `- n' o3 y; V4 ?8 ?+ A. \
图表 7 程序模拟单端测序
* I) ^: o+ Y/ a' H  |8 n6 U双端测序

图表 8 程序模拟双端测序
/ {8 U0 S- p& L6 u: A测序结果

5 d$ k) M! `; `4 B& P5 X: `图表 9 结果文件

因为期望片段长度是600bp，在片段长度区间200-1000bp内，所以大部分的片段都没有删除。
1 j# l# J, d5 c. Y测序结果统计表

4 P. o6 {6 I7 ^$ r( t测序方式        基因组大小(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
! e5 a: H: U5 i; E7 l' n四、讨论和结论
$ g. w. k' q0 ?) v5 z! R% p8 L
程序运行方法
; I) S5 Y  A+ x- y8 b" z) W
在类的构造方法init()中，调整参数。
! G$ P9 V: p3 ?2 L, X. cAveragefragmentlength为片段平均的长度；
$ b5 ~, |7 K2 ], u) k3 u8 Nminfragmentlength和maxfragmentlength是保留片段的范围；
3 b/ z+ S& Z% a2 ncloneRetainprobability是克隆的保留率；
minreadslength和maxreadslength是测序reads的长度范围
9 I+ B3 T0 N/ u+ l; u
& w' \# S$ |8 \模拟测序的诸多方法都封装成了Sequencing类，只需要创建类，并调用singlereadsequencing()和pairreadsequencing()方法，传入文件名的参数即可。

$ x1 `) A$ a6 `7 z/ r* R( x附录
% q6 m# I& p# e) m4 _4 S
from Bio import SeqIO
from math import exp
7 k8 D% r& v2 q0 E, C" }$ pimport random

class Sequencing:
( o( d0 P+ ~/ t5 _) P    # N代表拷贝份数
    def __init__(self)
        self.fragmentList = []
: i8 Y9 p7 J4 D9 \        self.readsID = 1
        self.readsList = []
- G; T6 o; I# d; Y- Z, O. V/ R        self.averagefragmentlength = 650
        self.minfragmentlength = 500
2 @; x+ ~% K9 g7 c, J        self.maxfragmentlength = 800
) O2 y, z3 e1 e& L9 S        self.cloneRetainprobability = 1
0 e4 w  S0 r9 C! A. [! n5 Y. q        self.minreadslength = 50
        self.maxreadslength = 150
        self.N = 10
" c2 q4 K7 X* r  T, R' f        self.genomeLength = 0
% W4 Q$ w6 c$ p9 c: s) V  K        self.allreadslength = 0
6 S: v0 s& |( Q: V
    # 生成断裂点
    def generatebreakpoint(self, seqlen, averageLength):
+ P, P1 S+ C7 l3 C1 ?- Q1 J7 C7 [2 n3 Z        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
( h, ^% \( o8 g+ @. A        breakpoint.append(seqlen)
5 W5 W/ D/ P6 e        breakpoint.append(0)
. C8 A+ l  ^  q! ?        # 把随机断裂点从小到大排序
        breakpoint.sort()
        return breakpoint
& B& q8 Y, a- T
    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
        for i in range(len(breakpoint) - 1):
8 N4 a& F  T" P            fragment = seq[breakpoint:breakpoint[i + 1]]
0 w' {( V1 v. D, K' h/ `; u: [" P            if maxfragmentlength > len(fragment) > minfragmentlength:
6 }2 A: M2 C0 t; ~" A3 F* X5 K                self.fragmentList.append(fragment)
/ c6 `4 J+ `) ], B/ I( i5 i- l        return self.fragmentList
  h) e) S# ]" i: u, S8 F4 l9 W
    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
; o* H* ]; w$ v+ w    def clonefragment(self, fragmentList, cloneRetainprobability):
) K7 K/ e) o$ s) O, ^0 n+ i        clonedfragmentList = []
/ d: V" e1 p* S& o* d        Lossprobability = [random.random() for _ in range(len(fragmentList))]
7 o5 F/ G6 i1 t. _0 M1 Y        for i in range(len(fragmentList)):
" b) J! C  R+ l  n            if Lossprobability <= cloneRetainprobability:
: H0 n1 a; e: t1 t4 Q) b                clonedfragmentList.append(fragmentList)
        return clonedfragmentList
% Z& Z9 Z) g4 [, _* G( y( |5 w( V
    # 模拟单端测序,并修改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
  ?& p7 i. d$ ]. {1 Q+ |  V            self.readsID += 1
            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])

* w9 X, }  s6 G* n1 Q    def singlereadsequencing(self, genomedata, sequencingResult):
. b( z2 N+ M& T4 e9 v, l! r        for seq_record in SeqIO.parse(genomedata, "fasta"):
& v! M* |' i& B            seqlen = len(seq_record)
            self.genomeLength += seqlen
            for i in range(self.N):
                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
% k7 c. o: \* D# j        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟单端测序
2 j& E1 [$ r/ b: f5 m: f        self.singleread(clonedfragmentList)
4 Z4 E0 E, w- Y( {        SeqIO.write(self.readsList, sequencingResult, "fasta")

2 |- i% C( n9 g$ D/ u    def pairread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
            fragment.id = ""
            fragment.name = ""
" w% B- l3 N; Y6 Q1 h% v            description = fragment.description[12:].split(",")[0]
/ m" Z3 Y1 w1 G; L' o            fragment.description = str(self.readsID) + "." + description
            readslength = random.randint(self.minreadslength, self.maxreadslength)
6 [0 G# b9 J! u            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])
6 T: g7 t5 w7 s7 M- O9 K1 `# ~) c( u
            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
. D  X3 k  c4 u
            fragmentcomplement = fragment.reverse_complement()
( l$ |; s7 V# E+ P6 h  m4 O& L            fragmentcomplement.id = ""
! W* }0 D$ j+ U* @8 y5 A8 G5 V            fragmentcomplement.name = ""
+ F5 S+ z- K- _8 v" f            fragmentcomplement.description = str(self.readsID) + "." + description
            self.readsList.append(fragmentcomplement[:readslength])
9 k1 H; w6 L! D. W
$ I% E/ N1 S- v8 ?            self.readsID += 1

& x1 Z/ Q3 J1 P) O0 f    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
( j# K8 V5 V, n2 A, X, h0 G        for seq_record in SeqIO.parse(genomedata, "fasta"):
. q4 e8 Y7 Y7 U% i9 }) c            seqlen = len(seq_record)
/ G: p% A9 @* l/ K$ b! @            self.genomeLength += seqlen
4 q+ l5 ?8 m1 x" q' s3 E0 ?            for i in range(self.N):
                # 生成断裂点
9 p7 o) m1 H& k8 A& D! `% A8 W$ K                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
) ~" f$ `+ f, g. i& V1 i  u" n                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
3 ?! M! F" X1 }' `% u- _        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟双端测序
        self.pairread(clonedfragmentList)
3 v0 G: i" A$ \, v. q: j, O& F        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]
- E- C# W, ]3 P% ^0 F" |        SeqIO.write(readsList_1, sequencingResult_1, "fasta")
! j& V" Q- v. G9 B% {        SeqIO.write(readsList_2, sequencingResult_2, "fasta")
+ R$ l4 @# E! ~/ I6 h$ A* X; ^, N
    def resultsummary(self):
        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
        print("N值:" + str(self.N))
        print("期望片段长度:" + str(self.averagefragmentlength))
: l5 n, T7 c: G- g        print("克隆保留率:" + str(self.cloneRetainprobability))
        print("片段数量:" + str(len(self.fragmentList)))
: f4 _6 P; w0 a2 H5 @        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
        print("reads总数量:" + str(len(self.readsList)))
        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
9 a5 h5 B, p8 R. N: \        m = self.allreadslength / self.genomeLength
& |5 a! P" E6 L' Q% t* |; M        print("覆盖度(m值):" + str(round(m, 5)))
        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
- Y6 e3 n* [/ _5 l4 F7 b" h4 \        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
/ V3 `+ ~; X( t1 K4 `- H# -------------------------------------------主程序-------------------------------------------
) R1 O% _+ g9 c. {# 模拟单端测序
sequencingObj = Sequencing()
7 V  H- }$ K2 i3 F* ysequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
2 D; F) m0 d' x; ~, S6 YsequencingObj.resultsummary()

# 模拟双端测序
) a  S6 Z1 w6 Z* q' y0 TsequencingObj = Sequencing()
sequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
sequencingObj.resultsummary()
from Bio import SeqIO
# x$ E" D8 [& ~' A3 b$ Z, P7 kfrom math import exp
import random
) V1 ~! I. n4 K# Q. i
class Sequencing:
6 s' O5 i; s' y6 l$ v    # N代表拷贝份数
' t8 N0 e9 \+ l. k6 W8 u& ^& r    def __init__(self):
7 Y2 _" N" P7 L. {        self.fragmentList = []
        self.readsID = 1
0 P4 J2 M8 Q, t! Z+ F2 L        self.readsList = []
        self.averagefragmentlength = 650
        self.minfragmentlength = 500
        self.maxfragmentlength = 800
        self.cloneRetainprobability = 1
        self.minreadslength = 50
        self.maxreadslength = 150
        self.N = 10
        self.genomeLength = 0
/ @, H  i. b& F6 P1 k3 k! a: j        self.allreadslength = 0
4 [2 W7 N- [" @/ A
8 p3 J3 L6 V6 v+ }    # 生成断裂点
    def generatebreakpoint(self, seqlen, averageLength):
        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
+ y: y% u) ?" r3 C8 x) v" W        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
        breakpoint.append(seqlen)
        breakpoint.append(0)
) O7 V) j/ o2 N, v4 G        # 把随机断裂点从小到大排序
        breakpoint.sort()
        return breakpoint
0 p5 ~! e# U, R
    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
- o- \' I- Z" q; r. W    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
        for i in range(len(breakpoint) - 1):
            fragment = seq[breakpoint:breakpoint[i + 1]]
+ B) }6 w* [0 I! f            if maxfragmentlength > len(fragment) > minfragmentlength:
' I* z7 r1 ?7 q1 I9 I. t6 p3 t  z                self.fragmentList.append(fragment)
        return self.fragmentList

    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
    def clonefragment(self, fragmentList, cloneRetainprobability):
        clonedfragmentList = []
        Lossprobability = [random.random() for _ in range(len(fragmentList))]
4 w+ N) F! x6 W& q        for i in range(len(fragmentList)):
2 g4 ]) ?4 O. S& d# x) \  z            if Lossprobability <= cloneRetainprobability:
                clonedfragmentList.append(fragmentList)
+ c$ Y9 k$ |+ k2 Y! g        return clonedfragmentList

    # 模拟单端测序,并修改reads的ID号
    def singleread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
            fragment.id = ""
            fragment.name = ""
, z: u- D9 B- U# V6 X            fragment.description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + fragment.description
            self.readsID += 1
+ S+ ~3 h4 P6 `1 }            readslength = random.randint(self.minreadslength, self.maxreadslength)
; [; }* f7 Y! y$ t+ e8 y4 Z            self.allreadslength += readslength
; h& ~% M/ H4 q! u7 q8 x. m2 W            self.readsList.append(fragment[:readslength])
7 R" B# T  P$ N, C0 |% Y
$ ]6 e7 M* Y- ?. r; f' [# p    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):
' a9 o. g) E3 r) [, _& O                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
/ h" |& d9 q" r1 x) h9 K: n                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
5 G6 c/ G4 q$ T7 O. E9 @  ]) a        # 模拟克隆时的随机丢失情况
4 W' ]; I! C# o' @1 @1 ~& h% X        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟单端测序
& t, Y+ K) h# l% p7 I2 K; N7 q9 d        self.singleread(clonedfragmentList)
. F% ~0 T; P" a. ^        SeqIO.write(self.readsList, sequencingResult, "fasta")

: t9 g* l# T( {: Z$ y' {/ C0 |    def pairread(self, clonedfragmentList):
5 d7 T5 H  \! N3 w' N        for fragment in clonedfragmentList:
; J4 [# A7 A* e+ p) y            fragment.id = ""
            fragment.name = ""
            description = fragment.description[12:].split(",")[0]
  [$ X- `" g4 z( t  l  e            fragment.description = str(self.readsID) + "." + description
; \4 {5 K" t; r6 u* P2 z" E            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])

) C! W, i% D7 }+ N5 H/ U1 j+ y            readslength = random.randint(self.minreadslength, self.maxreadslength)
2 L+ g4 V% y# |( I            self.allreadslength += readslength

4 w* C" m" v. O; Z' l/ C  @            fragmentcomplement = fragment.reverse_complement()
            fragmentcomplement.id = ""
            fragmentcomplement.name = ""
            fragmentcomplement.description = str(self.readsID) + "." + description
            self.readsList.append(fragmentcomplement[:readslength])
3 r& b: `$ V$ X# B& y# o
; @( z' D4 L' C7 }# e            self.readsID += 1

/ z3 {# s- Y; I( H  r7 Z. u$ `    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
        for seq_record in SeqIO.parse(genomedata, "fasta"):
  H7 g5 G7 ~, B3 u# L2 k            seqlen = len(seq_record)
            self.genomeLength += seqlen
2 J+ x2 R  c  G) ]/ v) ]2 v            for i in range(self.N):
                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
- ^4 h; r, U% c* X                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
6 A. t7 ]% x0 j" M- ~1 F8 }        # 模拟双端测序
        self.pairread(clonedfragmentList)
        readsList_1 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 0]
/ @+ {$ f8 }2 ^- j! H: V' a7 ~        readsList_2 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 1]
        SeqIO.write(readsList_1, sequencingResult_1, "fasta")
        SeqIO.write(readsList_2, sequencingResult_2, "fasta")
- Y, T- m6 H% N3 |) x  t' O4 d
    def resultsummary(self):
# q5 n& J( M$ }# Z3 t        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
        print("N值:" + str(self.N))
        print("期望片段长度:" + str(self.averagefragmentlength))
        print("克隆保留率:" + str(self.cloneRetainprobability))
        print("片段数量:" + str(len(self.fragmentList)))
        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
        print("reads总数量:" + str(len(self.readsList)))
6 V: p9 w1 }: b% C$ ^5 a& {. G        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
& U  \! ^# |# Y0 F        m = self.allreadslength / self.genomeLength
( m* x$ L# ]' z% |- R& v2 U6 g        print("覆盖度(m值):" + str(round(m, 5)))
5 v4 P. ]# @! l9 b4 k5 s1 N* e2 N        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
( c+ [. h, O2 \0 [- p# -------------------------------------------主程序-------------------------------------------
( e# R5 M5 w" x5 P9 z# 模拟单端测序
* I# _2 m' x  w3 N' s' ^sequencingObj = Sequencing()
sequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
5 v; X, z) M! S4 z1 j6 W! D1 g) N- wsequencingObj.resultsummary()

9 I2 s5 O' T! P6 ]6 ^5 o0 }" n# 模拟双端测序
- U) w2 M  i3 w2 L$ `  a: vsequencingObj = Sequencing()
sequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
sequencingObj.resultsummary()
3 H: o- O  \3 g
( m, j$ Q7 p" u0 @0 f2 W
( X' L: j6 d8 W- j+ I( R2 Q# a
  ~- x8 _, I2 V

3477959497

不错。。。。。。。。。。。。。。。。。
: q5 B2 h) s" g3 f1 e( O) i& l