基因组测序模拟

杨利霞

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

一、摘要

通过熟悉已有的基因组测序模拟和评估程序，加深全基因组鸟枪法测序原理的理解，并且能够编写程序模拟全基因组鸟枪法测序，理解覆盖度、测序深度、拷贝数等概念，设置测序相关参数，生成单端/双端测序结果文件
/ B6 }/ H$ `3 o' f$ s$ C9 K  `3 Y
* y: U* w) B" ?$ \9 B二、材料和方法

+ ~  Z3 c' t& P3 x5 ?  O% Q1、硬件平台
- }  n% H! u* L1 L/ B+ O% K
处理器:Intel（R） Core（TM）i7-4710MQ CPU @ 2.50GHz
安装内存(RAM)：16.0GB

; L6 S" Z5 Y+ `7 j9 H% [6 U2、系统平台
Windows 8.1，Ubuntu

* G! H7 {7 `0 |; R  ?' T. N0 T3、软件平台
+ O" I) K; P9 G' y+ o$ f8 @( H
art_454
1 U- j) B& ~& E/ P: w0 P: EGenomeABC http://crdd.osdd.net/raghava/genomeabc/
Python3.5
Biopython
4 u% f2 `5 k; J' w4、数据库资源

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

5、研究对象

3 W% \: w- j; Y5 t/ t2 I7 U酵母基因组Saccharomyces cerevisiae S288c (assembly R64)
& H* x  ]. Q2 _; n! M( J5 Uftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_genomic.fna.gz

2 ^9 n. \* M6 l' A: G% V6、方法

" m$ w2 w. r& p8 N. kart_454的使用
0 t' k$ S* {5 S: x: Z3 U  ~. s, [首先至art系列软件的官网，下载软件，在ubuntu系统安装，然后阅读相关参数设置的帮助文档，运行程序。
' ?$ B& @, w/ K; rGenomeABC
; X: l' ^  B, h" Q6 u- a) Z0 {进入GenomeABC(http://crdd.osdd.net/raghava/genomeabc/)，输入参数，获得模拟测序结果。
. c4 L" a9 x! l# R' ^编程模拟测序
下载安装python，并且安装biopython扩展模块，编写程序，模拟单端/双端测序。
- ]5 E5 }9 i. \! J7 R1 ]- H三、结果
9 N" V* j0 x$ Q  P9 C( ^3 ^* V
% X2 c; B- k- s* V) v5 l* w1、art_454的运行结果
5 f6 a% k8 u/ C- _- i
无参数art_454运行，阅读帮助文档

图表 1无参数art_454运行
1 I% I% u, I* L+ D, e对酵母基因组进行基因组单端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20.
下图为模拟单端测序，程序运行过程及结果
) P( j4 ?) T; F# Q0 v/ _9 J; i
0 l% x! Y3 ^0 c/ n+ q$ M/ s. e5 P图表 2 art454单端测序
/ r% j8 M& k; h3 v4 }! x
3 U, h9 ]2 _; b图表 3 art454单端模拟结果
/ P! v0 u; w0 X+ J" z& K2 A+ S6 |/ E1 [% v双端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20；MEAN_FRAG_LEN设为1500，即平均片段长度为1500；STD_DEV设为20，即长度的标准差为20
下图为模拟双端测序，程序运行过程及结果

图表 4 art454双端测序
* K2 y- q" _! D. }5 B+ o
图表 5 art454双端模拟结果
2、GenomeABC
下图为设置参数页面

下图为结果下载页面

图表 6 结果下载页面
% m7 J3 W2 A2 V& o2 t  q9 G$ T6 F5 o3、编程模拟测序结果
; H8 s( t* E- Z! {* C拷贝数是这里的N值；覆盖度是m，测序深度是宏观的量，在这里与覆盖度意思相同，就是测序仪10X，20X。
! N: F: k7 i' D5 u. p单端测序

图表 7 程序模拟单端测序
双端测序

- z5 J6 Q1 l0 c& n, D% s图表 8 程序模拟双端测序
. ~( ^- {$ q4 I! `* `' Y) F测序结果

1 g9 ~( _( W  h% @( G图表 9 结果文件

& N! L! }0 |/ y" g# X因为期望片段长度是600bp，在片段长度区间200-1000bp内，所以大部分的片段都没有删除。
测序结果统计表
% A4 n' i9 \) h$ J- M! E9 [, j
! s2 f  R4 _+ Y. Y测序方式        基因组大小(bp)        片段长度区间 (bp)        N值        期望片段长度        克隆保留率        片段数量        Reads长度范围（bp）        Reads总数量        Reads总长度        覆盖度(m值)        理论丢失率(e-m)        覆盖率(1-e-m)
0 {, |* Q3 S+ h5 K2 E& z) U& c9 Y% W单端        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
% M; A8 d8 R9 O6 s双端        12157kb        200-1000        10        600        0.95        106704        50-100        202770        15212.662kb        1.25134        0.28612        0.71388
$ @! |& H8 u5 r+ G0 T双端        12157kb        200-1000        20        600        0.95        214212        50-100        407186        30534.265kb        2.51164        0.08114        0.91886
7 K9 m# e: ~# `: [) L- _8 X! ]四、讨论和结论

程序运行方法

在类的构造方法init()中，调整参数。
Averagefragmentlength为片段平均的长度；
( e& h& k$ P3 xminfragmentlength和maxfragmentlength是保留片段的范围；
cloneRetainprobability是克隆的保留率；
) F' E9 ]# O8 Ominreadslength和maxreadslength是测序reads的长度范围
& T8 S7 {# Q( Z! s! N; [
模拟测序的诸多方法都封装成了Sequencing类，只需要创建类，并调用singlereadsequencing()和pairreadsequencing()方法，传入文件名的参数即可。

附录
/ @, B, z+ Z' Z3 C5 B; I: g7 H
from Bio import SeqIO
; f  j4 q9 ?) K, K( q7 {( C; R% nfrom math import exp
import random

, I9 Z" o" M9 @# r5 G4 S- Hclass Sequencing:
    # N代表拷贝份数
    def __init__(self)
        self.fragmentList = []
        self.readsID = 1
/ \% T4 K: w9 g) R        self.readsList = []
        self.averagefragmentlength = 650
* Q- m" M- l( [  f% w! f        self.minfragmentlength = 500
        self.maxfragmentlength = 800
        self.cloneRetainprobability = 1
5 G" T$ A' Y5 N1 o+ G        self.minreadslength = 50
        self.maxreadslength = 150
        self.N = 10
        self.genomeLength = 0
        self.allreadslength = 0
: @/ t' O% v7 w$ f9 Y* z
    # 生成断裂点
. c* n7 p! H0 K$ B    def generatebreakpoint(self, seqlen, averageLength):
+ u) m* `. s: e& Q2 p# [; H" \0 j        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
        breakpoint.append(seqlen)
1 d' B4 j! F7 {' @) L5 D7 `4 R; Q        breakpoint.append(0)
        # 把随机断裂点从小到大排序
  w5 A- E# |0 g$ K        breakpoint.sort()
        return breakpoint
$ W' x0 H7 \' X) X) I1 ]$ ~9 a
    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
        for i in range(len(breakpoint) - 1):
4 \+ r8 Q# R- S: r            fragment = seq[breakpoint:breakpoint[i + 1]]
            if maxfragmentlength > len(fragment) > minfragmentlength:
8 C8 Y0 g3 A2 p& `& O                self.fragmentList.append(fragment)
        return self.fragmentList

    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
    def clonefragment(self, fragmentList, cloneRetainprobability):
        clonedfragmentList = []
, \8 R' g! E' L5 f& p' e5 B% t1 \        Lossprobability = [random.random() for _ in range(len(fragmentList))]
        for i in range(len(fragmentList)):
6 D' z$ f  b( \5 D9 L8 ~. Q            if Lossprobability <= cloneRetainprobability:
3 |9 |* Y& M9 m. F* ]                clonedfragmentList.append(fragmentList)
        return clonedfragmentList

& O. Q! P8 Z8 M6 m0 }. L    # 模拟单端测序,并修改reads的ID号
3 Z5 j0 m# }  z2 E- F. _    def singleread(self, clonedfragmentList):
  ^0 ~2 P% S: A" V1 f5 V! N) y$ p        for fragment in clonedfragmentList:
            fragment.id = ""
            fragment.name = ""
            fragment.description = fragment.description[12:].split(",")[0]
' f  Y  E0 v2 g3 R8 O4 v            fragment.description = str(self.readsID) + "." + fragment.description
            self.readsID += 1
            readslength = random.randint(self.minreadslength, self.maxreadslength)
& m6 D. Z6 b" X            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])

) l7 i- a3 u) s. P# L* V/ F    def singlereadsequencing(self, genomedata, sequencingResult):
        for seq_record in SeqIO.parse(genomedata, "fasta"):
/ |4 A4 \8 \; S8 k* R2 Y; l9 D            seqlen = len(seq_record)
            self.genomeLength += seqlen
            for i in range(self.N):
/ I  f8 y% ~- d. y0 U) q! i, E                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
& [, Q3 I) `* s                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
! I$ K8 y! E" @: K5 M6 F        # 模拟单端测序
& i' s* h: H  r# s- T' g% `        self.singleread(clonedfragmentList)
        SeqIO.write(self.readsList, sequencingResult, "fasta")

    def pairread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
( _" m, P# A3 s3 I% N) Q            fragment.id = ""
            fragment.name = ""
            description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + description
. W$ \8 ~. l2 D5 q7 K            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])
9 R6 {5 ]2 V/ L+ i: g/ `$ J) ~' q
            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
" j7 I) i0 V& F
            fragmentcomplement = fragment.reverse_complement()
* Y* Y+ |8 H- s9 l            fragmentcomplement.id = ""
            fragmentcomplement.name = ""
, u. J6 _& j0 a8 @2 ]  n            fragmentcomplement.description = str(self.readsID) + "." + description
8 _& F# ^4 X. y) m: T            self.readsList.append(fragmentcomplement[:readslength])
* [+ E9 L' e: Y, \3 c& O1 I/ @
5 u2 j# O. K) N7 T+ X8 X2 v            self.readsID += 1
- R  l# q" V; o3 U8 e
% U; G. t: h6 C. A7 B' B7 B9 ]    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
            self.genomeLength += seqlen
            for i in range(self.N):
                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
$ T% m) J+ }% {                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
5 ?3 L7 ~1 @% A. z        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
9 A- @; o3 e2 A/ Z; w( A) @1 A        # 模拟双端测序
4 O  p' [2 l, [/ \        self.pairread(clonedfragmentList)
        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")
% \5 {5 U# s" f* _7 j  Q% Z        SeqIO.write(readsList_2, sequencingResult_2, "fasta")
: k6 {- }  y( w  n. Y" r) M3 `
* A2 B8 k7 f: W6 H# d    def resultsummary(self):
% y$ }( p* ?1 r        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
        print("N值:" + str(self.N))
+ {8 U9 R! n8 f, p8 y) `8 l; z% y        print("期望片段长度:" + str(self.averagefragmentlength))
        print("克隆保留率:" + str(self.cloneRetainprobability))
        print("片段数量:" + str(len(self.fragmentList)))
6 ~3 E# Q: c1 h        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
        print("reads总数量:" + str(len(self.readsList)))
; Y; K' a# H# q9 G        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
        m = self.allreadslength / self.genomeLength
9 ?3 t( e& g* z9 c# I+ T7 G        print("覆盖度(m值):" + str(round(m, 5)))
        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
4 _! A& c! P8 G+ l6 ^, K. d/ D# -------------------------------------------主程序-------------------------------------------
/ @. a3 Y) M. ?6 e. d# B7 a; @# 模拟单端测序
sequencingObj = Sequencing()
* M. l- Q( u; y  i4 }9 KsequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
sequencingObj.resultsummary()

9 `  q2 Y( |" u6 d5 D- T# 模拟双端测序
$ O; U9 e* c, F! GsequencingObj = Sequencing()
sequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
sequencingObj.resultsummary()
* F+ F% F: i% r2 y; I( ^from Bio import SeqIO
from math import exp
: l, |/ C% x4 l0 I5 `import random
0 ]) ~; s! G+ h8 J
class Sequencing:
    # N代表拷贝份数
    def __init__(self):
- X+ R! Q; f( p8 l) R" Y5 Z1 Y( |        self.fragmentList = []
        self.readsID = 1
        self.readsList = []
1 J4 N  g' J- r$ T, i' X) ~        self.averagefragmentlength = 650
        self.minfragmentlength = 500
* ~3 S6 i$ A. R" Q. Z! ]1 x        self.maxfragmentlength = 800
9 q9 x# S: F" v# o& M# y, T1 b        self.cloneRetainprobability = 1
        self.minreadslength = 50
        self.maxreadslength = 150
. x  O& O  V" e        self.N = 10
        self.genomeLength = 0
9 `$ J6 U8 X3 J  ]! E! ^; U# e        self.allreadslength = 0

    # 生成断裂点
6 y- f) a. K; G8 w: `0 H9 {3 y    def generatebreakpoint(self, seqlen, averageLength):
2 J% H# q3 x) `9 M* A        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
) k9 N- y1 q9 x- w        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
5 a+ x6 K* V) g        breakpoint.append(seqlen)
        breakpoint.append(0)
! v1 _3 ]9 {" w6 @" u' U0 x0 Q        # 把随机断裂点从小到大排序
6 i' J6 N6 d7 k$ ~        breakpoint.sort()
        return breakpoint
, P+ k1 M1 b+ j+ b
    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
# J4 u' h, L+ g# L$ I+ ~        for i in range(len(breakpoint) - 1):
            fragment = seq[breakpoint:breakpoint[i + 1]]
1 D3 _7 Q- O( Y! D7 }2 z1 P2 s            if maxfragmentlength > len(fragment) > minfragmentlength:
* v. j/ F2 e+ a( V: p                self.fragmentList.append(fragment)
( ~7 D+ Q7 w$ T        return self.fragmentList

    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
    def clonefragment(self, fragmentList, cloneRetainprobability):
0 W! Y' Z9 K8 D  F: k        clonedfragmentList = []
        Lossprobability = [random.random() for _ in range(len(fragmentList))]
) E; G" Z: R( t; p        for i in range(len(fragmentList)):
: R  t! G4 L" |0 a! g' {9 x8 K* d            if Lossprobability <= cloneRetainprobability:
                clonedfragmentList.append(fragmentList)
: ^3 A/ z# F# {        return clonedfragmentList

    # 模拟单端测序,并修改reads的ID号
    def singleread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
            fragment.id = ""
            fragment.name = ""
+ W( q& i1 v0 T) Z" B            fragment.description = fragment.description[12:].split(",")[0]
3 z6 O9 ?# H8 X3 _7 H: P1 x            fragment.description = str(self.readsID) + "." + fragment.description
            self.readsID += 1
            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
6 A9 F- H! s" Q& `+ }: g9 N, X# r            self.readsList.append(fragment[:readslength])
/ m* P$ d  X; j6 a: \% l& q
    def singlereadsequencing(self, genomedata, sequencingResult):
* k3 N7 H/ E% z3 n# g8 |        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
2 _: O8 N) @; p8 }1 W+ V: Q+ o2 {7 q0 {            self.genomeLength += seqlen
5 f  t6 b4 u1 D* F            for i in range(self.N):
( e7 o! `( g) h5 [5 ~                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
8 A5 {' U( [: b) Z! C7 W$ e4 p: c                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
9 d% x; Q+ _) l0 R        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
$ r2 [1 R& P! R; y" q3 b2 |" y        # 模拟单端测序
        self.singleread(clonedfragmentList)
7 ?, b$ Y1 A% T) r  h. o( E' S* a        SeqIO.write(self.readsList, sequencingResult, "fasta")
: W" _. s$ g; N9 F7 U6 }6 a) ~
6 q) `! x! u: Z; B$ r    def pairread(self, clonedfragmentList):
- L) u) o1 U. }& m        for fragment in clonedfragmentList:
            fragment.id = ""
0 m" u. l/ b. n0 j  H" v            fragment.name = ""
            description = fragment.description[12:].split(",")[0]
% h& b) f9 }* K$ n' }) f& q, a  g# r9 Q: ]            fragment.description = str(self.readsID) + "." + description
            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
! ]: Z' s) A6 i            self.readsList.append(fragment[:readslength])

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

) I; U7 O. Y* Z% j' j+ r" i            fragmentcomplement = fragment.reverse_complement()
* ^% b8 D7 c- e7 m8 M0 F            fragmentcomplement.id = ""
            fragmentcomplement.name = ""
            fragmentcomplement.description = str(self.readsID) + "." + description
            self.readsList.append(fragmentcomplement[:readslength])

            self.readsID += 1
0 n3 U( W- O3 D7 N- k
    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
% @0 \$ l: Q& n. S4 v        for seq_record in SeqIO.parse(genomedata, "fasta"):
, r/ T) l( v- y. i8 r            seqlen = len(seq_record)
" Q1 {6 e' C; z1 R            self.genomeLength += seqlen
            for i in range(self.N):
/ q* U2 s  K. i$ C- j* u7 B                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
! h. t5 ?5 W2 p  S: O9 Q% k  W        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
& R% d" b/ T/ |6 t( m# T8 E        # 模拟双端测序
        self.pairread(clonedfragmentList)
        readsList_1 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 0]
/ V! L. [9 Y! K9 n' I+ m        readsList_2 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 1]
' o$ z: M/ H# q$ l# ~) J; `; A        SeqIO.write(readsList_1, sequencingResult_1, "fasta")
        SeqIO.write(readsList_2, sequencingResult_2, "fasta")

    def resultsummary(self):
2 C9 M) Z5 G3 A; [( R/ f% G        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
; y2 p6 r  W. g. ]        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
) U$ I$ f8 }1 F" p& [        print("N值:" + str(self.N))
, H0 q, A7 g) s        print("期望片段长度:" + str(self.averagefragmentlength))
* Z4 x  I  W; L/ A8 F        print("克隆保留率:" + str(self.cloneRetainprobability))
; K$ w- m% l+ F5 e* v! y& _        print("片段数量:" + str(len(self.fragmentList)))
9 f' F; F  Y3 a0 c        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
. @# C; `" C! I6 K1 z        print("reads总数量:" + str(len(self.readsList)))
        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
        m = self.allreadslength / self.genomeLength
$ v7 s3 [( B0 Z6 ]4 i3 G5 T( N        print("覆盖度(m值):" + str(round(m, 5)))
        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
9 b+ ?: `) k& C4 Z        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
# -------------------------------------------主程序-------------------------------------------
; W, t% |% ^5 T1 z  z# 模拟单端测序
) o) x/ q% k  w+ @sequencingObj = Sequencing()
) S# P! _! y- v+ F* {+ |sequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
1 Q+ q% ?  O/ x$ U( EsequencingObj.resultsummary()
# l, Q3 @! ^7 g. ?2 P
# 模拟双端测序
1 K7 ~2 J- j: usequencingObj = Sequencing()
/ J: ]* S0 M0 g4 l1 O  s% ~sequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
$ k# `% [3 D7 Q9 jsequencingObj.resultsummary()
( K7 p5 H7 |! b


& F. K* K$ g! k+ f6 ^  w5 u$ K
" ^6 ?; m2 `1 R: T3 H/ e  g( n

3477959497

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