基因组测序模拟

杨利霞

基因组测序模拟
基因组测序模拟
" J: ?8 R$ c/ @* {
4 c3 a) Y8 s% c一、摘要
& H5 T% O/ @; P
% X6 R# Z) M/ e* J/ w3 _3 v通过熟悉已有的基因组测序模拟和评估程序，加深全基因组鸟枪法测序原理的理解，并且能够编写程序模拟全基因组鸟枪法测序，理解覆盖度、测序深度、拷贝数等概念，设置测序相关参数，生成单端/双端测序结果文件

, ?9 ?9 C# l2 ~; `9 `二、材料和方法

1、硬件平台

处理器:Intel（R） Core（TM）i7-4710MQ CPU @ 2.50GHz
7 e" k" w) V, k" J安装内存(RAM)：16.0GB

$ u- ~6 @: Q% @& U5 t" r0 a  A2、系统平台
Windows 8.1，Ubuntu

$ D2 O. i; i; ~1 E3、软件平台
! U% w. R; E4 ]  _7 b
3 g& p" t9 B4 a4 ?* C/ M' Sart_454
GenomeABC http://crdd.osdd.net/raghava/genomeabc/
5 r; `) X$ i3 ]9 g! APython3.5
Biopython
4、数据库资源

NCBI数据库：https://www.ncbi.nlm.nih.gov/
' F6 ]. X$ l. H& F* I7 _
9 E  K( {; Z; t. M3 v4 S3 c5、研究对象
+ [( S4 d% t3 ^8 T: c8 A6 Q
5 M" L! T* l; k" \% o  E酵母基因组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

6、方法

art_454的使用
& J  j, g8 F+ f6 _/ p1 @首先至art系列软件的官网，下载软件，在ubuntu系统安装，然后阅读相关参数设置的帮助文档，运行程序。
GenomeABC
进入GenomeABC(http://crdd.osdd.net/raghava/genomeabc/)，输入参数，获得模拟测序结果。
" ]2 l0 O. H( `编程模拟测序
0 L! E/ ^. N: m) d2 U! Z7 z下载安装python，并且安装biopython扩展模块，编写程序，模拟单端/双端测序。
+ t  C" F' K9 m# G! }1 ~三、结果
) _- x" P& ^9 t# K1 [
* S& T' k/ r8 H1、art_454的运行结果
  Q3 z2 W) J+ \7 D
无参数art_454运行，阅读帮助文档
% F: q; s6 B2 H6 j
5 B2 f" h4 m+ o! ]5 x& S$ g: Z2 a图表 1无参数art_454运行
对酵母基因组进行基因组单端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20.
8 N2 @/ v7 t* u4 x% Q下图为模拟单端测序，程序运行过程及结果

图表 2 art454单端测序
! I. ]+ G  M4 v, Z
3 {) f" y  V' ]0 }) k( m图表 3 art454单端模拟结果
双端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20；MEAN_FRAG_LEN设为1500，即平均片段长度为1500；STD_DEV设为20，即长度的标准差为20
5 G' t' k8 h1 n% H3 B! i1 s1 n# i; M" X下图为模拟双端测序，程序运行过程及结果

; H/ ^. H% T2 h5 z8 z图表 4 art454双端测序
) [! ^, ?: k4 R1 ^% w5 q# K5 s
图表 5 art454双端模拟结果
2、GenomeABC
下图为设置参数页面

下图为结果下载页面

5 @3 _" I1 |" \. [图表 6 结果下载页面
3、编程模拟测序结果
/ [4 ?2 d$ f/ B拷贝数是这里的N值；覆盖度是m，测序深度是宏观的量，在这里与覆盖度意思相同，就是测序仪10X，20X。
& e1 l, X' M) v! x: C单端测序
6 @4 k4 w$ j1 S
& O/ g4 F2 b! i: X2 c7 q0 J图表 7 程序模拟单端测序
双端测序

图表 8 程序模拟双端测序
测序结果
; u* `; ~, V! o0 N* d& g
! l! B+ ]5 |: v# r: U图表 9 结果文件

& B& K' D$ ?4 x' `% `. d1 A因为期望片段长度是600bp，在片段长度区间200-1000bp内，所以大部分的片段都没有删除。
测序结果统计表
5 o! f- G' N+ G: p, g3 t0 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
' O* P# K4 {$ T单端        12157kb        200-1000        20        600        0.95        213722        50-100        202996        15227.882kb        1.25259        0.28576        0.71424
5 ?$ L) f) a& _2 `$ e8 D双端        12157kb        200-1000        10        600        0.95        106704        50-100        202770        15212.662kb        1.25134        0.28612        0.71388
  |, A" g  a2 u/ `# Q# M双端        12157kb        200-1000        20        600        0.95        214212        50-100        407186        30534.265kb        2.51164        0.08114        0.91886
四、讨论和结论
9 g% V& N. ?+ G
1 L: W% h3 n6 ^- c6 J' u程序运行方法
! \, T8 C3 p, X
在类的构造方法init()中，调整参数。
Averagefragmentlength为片段平均的长度；
& p3 Y$ D( J& e' e" x1 zminfragmentlength和maxfragmentlength是保留片段的范围；
' c8 L/ X' |9 FcloneRetainprobability是克隆的保留率；
) T) {+ h& q+ Aminreadslength和maxreadslength是测序reads的长度范围
, R! _# [; q5 V9 V. O
  Z( ]. O6 H) }" u( R. u$ W1 j0 T模拟测序的诸多方法都封装成了Sequencing类，只需要创建类，并调用singlereadsequencing()和pairreadsequencing()方法，传入文件名的参数即可。

( N  ?. l! b- ?附录
+ ^- g4 J# [& J) O. ^) ]0 ?# l3 G% X
' F% a$ a' w. [$ P! H3 ]: Ufrom Bio import SeqIO
from math import exp
  e! K# A  d' J' ?7 \$ pimport random
. X0 z: P8 ~% V: j8 C+ g. R
class Sequencing:
    # N代表拷贝份数
0 I( u2 z& H" g0 h; H    def __init__(self)
0 [  S6 `4 L- @( ]. U        self.fragmentList = []
& z7 ]. e8 |0 f        self.readsID = 1
; J2 f0 T: k5 Z% J; Y/ y        self.readsList = []
        self.averagefragmentlength = 650
        self.minfragmentlength = 500
        self.maxfragmentlength = 800
* I) O, U1 {8 r        self.cloneRetainprobability = 1
        self.minreadslength = 50
- L5 g! Q+ l$ P: B9 b5 y        self.maxreadslength = 150
        self.N = 10
        self.genomeLength = 0
        self.allreadslength = 0
1 S2 n- k6 R5 Y- I% V) d
2 @: [- F( M- h/ p    # 生成断裂点
/ F3 H$ ^8 U% D& M    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)
        breakpoint.append(0)
) [5 D3 T9 E. B7 g4 p  |        # 把随机断裂点从小到大排序
6 V- W1 k' d; I8 l/ f" I6 D        breakpoint.sort()
        return breakpoint
$ R) \+ o5 o- z' Q) q* v
) v8 o% H6 R3 i    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
' E' x1 ~2 U% n& n* ^0 o+ M; @    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
. S; ^0 \& a0 Z" s4 I) i$ N& P        for i in range(len(breakpoint) - 1):
) P: m9 i# ?: }( l0 L. ?+ y4 L, N            fragment = seq[breakpoint:breakpoint[i + 1]]
            if maxfragmentlength > len(fragment) > minfragmentlength:
: \: U& {, m( z                self.fragmentList.append(fragment)
7 Z" r# C) F6 n: \        return self.fragmentList

    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
2 n& I$ `4 e# W0 P/ G: h6 W    def clonefragment(self, fragmentList, cloneRetainprobability):
        clonedfragmentList = []
        Lossprobability = [random.random() for _ in range(len(fragmentList))]
        for i in range(len(fragmentList)):
9 q+ R0 F7 I/ Y, \            if Lossprobability <= cloneRetainprobability:
+ ]' v6 Z4 f; b% r                clonedfragmentList.append(fragmentList)
8 Z  L, o3 F+ Y& M        return clonedfragmentList

/ N' }0 u& p* z% |  g5 W    # 模拟单端测序,并修改reads的ID号
4 c5 p! n/ G4 A1 Q+ [: ?    def singleread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
! a+ j% }! Q( M, g# X            fragment.id = ""
            fragment.name = ""
0 ?) B6 K$ U$ p2 c; G            fragment.description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + fragment.description
7 s  E5 t5 X0 A            self.readsID += 1
; \/ C1 D, P  a6 n( t            readslength = random.randint(self.minreadslength, self.maxreadslength)
% i( f6 h; q. L" \            self.allreadslength += readslength
/ n9 D4 U8 ^% j+ Z* M            self.readsList.append(fragment[:readslength])
4 @3 B( n- {1 z7 ]/ x6 u
    def singlereadsequencing(self, genomedata, sequencingResult):
* V& |+ W) r0 }. n8 i1 x        for seq_record in SeqIO.parse(genomedata, "fasta"):
+ h1 T. Q& w. ^' }1 `! l            seqlen = len(seq_record)
0 |1 K9 o0 Z4 W5 \% y5 ?            self.genomeLength += seqlen
; C2 i1 z! n8 b; [1 i            for i in range(self.N):
                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
: J+ N- R- a$ j& A3 q, B' D5 i                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
/ \- j4 `* }6 D' r8 E# ?$ F9 h/ D        # 模拟克隆时的随机丢失情况
2 G4 m2 N9 N8 E5 s0 B        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟单端测序
8 Z2 ~4 u: i! f        self.singleread(clonedfragmentList)
        SeqIO.write(self.readsList, sequencingResult, "fasta")

    def pairread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
9 R: W" _9 {& o$ q" T; I' Z+ Q            fragment.id = ""
) V4 @: I4 B" n9 t: I6 {: z% W            fragment.name = ""
9 L, m  |; c; N. C5 A            description = fragment.description[12:].split(",")[0]
. H1 p1 ^$ F0 U0 x( M# l& N            fragment.description = str(self.readsID) + "." + description
            readslength = random.randint(self.minreadslength, self.maxreadslength)
. M7 `' B& o& ?9 n. G            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])
, ^$ t" z0 N9 w, J2 O/ o8 N7 D
6 o  g, {9 E8 U! F8 K: U            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength

            fragmentcomplement = fragment.reverse_complement()
7 \) U6 N0 t8 t, X2 C$ u$ m            fragmentcomplement.id = ""
, `" d% N. n' S2 ~3 w5 i            fragmentcomplement.name = ""
4 D% [' j- Y# w& B+ n/ w            fragmentcomplement.description = str(self.readsID) + "." + description
0 r0 Z9 k: r6 k% }& j9 R            self.readsList.append(fragmentcomplement[:readslength])

            self.readsID += 1
# S% ^- v/ h. P- m9 n
/ f1 r/ w# T! r" e  D1 @2 q( V    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
        for seq_record in SeqIO.parse(genomedata, "fasta"):
  G( N  `9 `' v            seqlen = len(seq_record)
6 q5 b! V; c; G  {8 ]# p            self.genomeLength += seqlen
' h: W8 J' Q) u7 d/ J6 y# {            for i in range(self.N):
- o; O8 C; [1 V5 d7 f- y                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
+ n! e% o7 M: ]5 _- a4 s2 T* J2 k        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟双端测序
        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]
! J: t! ~/ I; P        SeqIO.write(readsList_1, sequencingResult_1, "fasta")
        SeqIO.write(readsList_2, sequencingResult_2, "fasta")

    def resultsummary(self):
  R, s: _9 L4 d( E3 \8 s, j        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))
* v2 E4 P9 \8 b! m* F        print("reads总数量:" + str(len(self.readsList)))
8 z7 I/ M) Q' V# q3 n        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
        m = self.allreadslength / self.genomeLength
        print("覆盖度(m值):" + str(round(m, 5)))
7 ], _$ s8 q7 b0 |3 D7 i" m- C        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
; y) @% O1 N! I: e        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
( d; N* D$ P# O/ R( k" Q/ x1 j  T# -------------------------------------------主程序-------------------------------------------
# 模拟单端测序
sequencingObj = Sequencing()
sequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
sequencingObj.resultsummary()

& ]+ S2 q0 w" D& }! |# 模拟双端测序
sequencingObj = Sequencing()
$ r, Q7 ]& \/ X5 psequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
sequencingObj.resultsummary()
from Bio import SeqIO
from math import exp
import random

; L- d& T2 z7 Y) g: _% D7 V  Q8 rclass Sequencing:
    # N代表拷贝份数
    def __init__(self):
$ w7 N( L$ x. p) z3 n2 e        self.fragmentList = []
        self.readsID = 1
        self.readsList = []
        self.averagefragmentlength = 650
! l' V6 K9 W: r0 _, }2 H/ P        self.minfragmentlength = 500
2 y6 A" D: |  L" Y/ F) p0 s        self.maxfragmentlength = 800
- ~: S, R* A% @        self.cloneRetainprobability = 1
0 d# Y; L! q9 [+ g2 R5 F4 b        self.minreadslength = 50
% _* g" U6 ~* N        self.maxreadslength = 150
        self.N = 10
        self.genomeLength = 0
        self.allreadslength = 0

    # 生成断裂点
) a% I6 z$ Y" k+ Y6 Q    def generatebreakpoint(self, seqlen, averageLength):
5 h5 M; E: v8 \9 f* x5 b        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
& q/ J3 G" `8 E( P        breakpoint.append(seqlen)
3 s- H$ q4 K* ?- y" ?% ~        breakpoint.append(0)
        # 把随机断裂点从小到大排序
1 B& J3 W# Z" R/ h        breakpoint.sort()
        return breakpoint
. }% @+ n! p0 J* X$ A/ ?6 x
) Z9 ~* c& G. a    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
0 W' I, Q" y. s+ ?3 {. K    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
4 L3 P9 e1 {7 Z        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
2 S6 A9 e% n( [% t0 C' Y
    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
    def clonefragment(self, fragmentList, cloneRetainprobability):
, u- U8 P/ g5 B. c: ]; f        clonedfragmentList = []
        Lossprobability = [random.random() for _ in range(len(fragmentList))]
        for i in range(len(fragmentList)):
            if Lossprobability <= cloneRetainprobability:
2 d4 H) j0 X7 e3 K3 I& S" \                clonedfragmentList.append(fragmentList)
        return clonedfragmentList
5 [3 S1 l. B  \, h8 N' }
0 t4 ~9 C- m  y+ g; \0 |/ K9 y. C# l    # 模拟单端测序,并修改reads的ID号
    def singleread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
            fragment.id = ""
3 c5 |6 S. Q9 ]9 O6 ~4 i            fragment.name = ""
            fragment.description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + fragment.description
            self.readsID += 1
            readslength = random.randint(self.minreadslength, self.maxreadslength)
- w. D# R* x7 w5 ~            self.allreadslength += readslength
' R# A5 P8 F9 R& ^            self.readsList.append(fragment[:readslength])
' p; [! b: h: N: k1 I$ O
    def singlereadsequencing(self, genomedata, sequencingResult):
/ ^! K$ l) k2 W0 U; Y, i        for seq_record in SeqIO.parse(genomedata, "fasta"):
  x2 u! H8 e4 O4 X3 d( M  `            seqlen = len(seq_record)
            self.genomeLength += seqlen
* d+ G9 L) }& ?. n            for i in range(self.N):
                # 生成断裂点
/ ^9 Q; h0 v0 N% q( i$ ?                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
        # 模拟克隆时的随机丢失情况
! i% D8 `! A' f! B- x0 g3 o2 \  m        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟单端测序
        self.singleread(clonedfragmentList)
+ x' c+ E% O' k" C; A        SeqIO.write(self.readsList, sequencingResult, "fasta")
* f# F" y* d! X" h2 b
    def pairread(self, clonedfragmentList):
5 ~: ^; f  w! Q% H' J        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)
* N( B, c4 D! R& x  O            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])

- ^" H3 x$ x3 x; o$ K$ \            readslength = random.randint(self.minreadslength, self.maxreadslength)
( M3 r7 z7 v8 Z" h            self.allreadslength += readslength
6 d: m( M: S4 N$ H8 l4 v7 {; v" w: X$ I
            fragmentcomplement = fragment.reverse_complement()
            fragmentcomplement.id = ""
            fragmentcomplement.name = ""
+ k. w# _8 f6 F            fragmentcomplement.description = str(self.readsID) + "." + description
            self.readsList.append(fragmentcomplement[:readslength])
) d' Q$ A- @( R8 J/ S
            self.readsID += 1

/ \; Y9 w5 P, Y2 R3 W4 a    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
( X1 ^( P  F  Y  q8 U  m* B, P        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
            self.genomeLength += seqlen
            for i in range(self.N):
                # 生成断裂点
3 ^1 z8 F; E& z$ j& O: c                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
6 H8 D2 @$ w; k9 f, Z9 M/ A+ P                # 沿断裂点打断基因组
5 c) x9 X, `8 B2 ?! X8 M                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
0 [# B& D8 X! s$ A$ H1 `- ]        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟双端测序
" X' }( @$ A* R0 g        self.pairread(clonedfragmentList)
        readsList_1 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 0]
& e+ z! S) W  f1 O* L; r        readsList_2 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 1]
$ u& _# @! B$ t5 v; _        SeqIO.write(readsList_1, sequencingResult_1, "fasta")
5 T6 h: F! P4 ^        SeqIO.write(readsList_2, sequencingResult_2, "fasta")

" E5 ?3 K# k% `3 ~" k/ u    def resultsummary(self):
2 v9 r  O$ W" }! e1 n7 w, U# J        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
! _# k! @) g7 ]# h        print("N值:" + str(self.N))
        print("期望片段长度:" + str(self.averagefragmentlength))
        print("克隆保留率:" + str(self.cloneRetainprobability))
, _" K- k, b' Y" g7 `        print("片段数量:" + str(len(self.fragmentList)))
* b! c5 C1 L4 F4 O% s; o# B, i        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
# ]' c# ^# c8 T. Y. f        print("reads总数量:" + str(len(self.readsList)))
+ m+ Z7 s% Y( |: w3 L3 K0 ^        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
        m = self.allreadslength / self.genomeLength
1 ]/ c. g2 h8 P! s$ E        print("覆盖度(m值):" + str(round(m, 5)))
        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
" V# S8 I8 d! ]& ^  A        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
# -------------------------------------------主程序-------------------------------------------
# 模拟单端测序
sequencingObj = Sequencing()
7 D% E) h& U* l0 F/ G3 w. TsequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
sequencingObj.resultsummary()
% o$ I  N# Q! I9 e6 k. W
; ~7 [8 ?4 ]1 z# 模拟双端测序
sequencingObj = Sequencing()
sequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
sequencingObj.resultsummary()



% q+ q2 ~# h1 t( b& ?

3477959497

不错。。。。。。。。。。。。。。。。。
6 T. S# c9 v1 S: Z5 h- H