基因组测序模拟

杨利霞

基因组测序模拟
. E" n% |8 s9 j# P; l基因组测序模拟

$ V6 [: v! }: {. D2 o一、摘要

通过熟悉已有的基因组测序模拟和评估程序，加深全基因组鸟枪法测序原理的理解，并且能够编写程序模拟全基因组鸟枪法测序，理解覆盖度、测序深度、拷贝数等概念，设置测序相关参数，生成单端/双端测序结果文件
" R% p+ H' C6 k' ?
) F$ l  r" \# F* s二、材料和方法
  x2 [1 ]9 ?. n9 k
1、硬件平台

9 ^6 \; s# r9 {& k0 Q处理器:Intel（R） Core（TM）i7-4710MQ CPU @ 2.50GHz
安装内存(RAM)：16.0GB

( |) H5 h2 d3 x. J9 r( W2、系统平台
Windows 8.1，Ubuntu
- ^6 }& J6 r0 K  [
3、软件平台

art_454
# _# I# n' U# u8 d" @GenomeABC http://crdd.osdd.net/raghava/genomeabc/
Python3.5
Biopython
4、数据库资源

3 u* s0 k; f- C% D- kNCBI数据库：https://www.ncbi.nlm.nih.gov/
- U# \, y" E! D0 Z. R+ w
5、研究对象

( i5 ?  J. S8 I3 ~' ?- Y7 \4 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
  t1 U. e' l- ?. H# F9 t1 s
0 h7 q0 S- H2 H+ ?' T( L6、方法

, T& `* e9 R9 A0 xart_454的使用
首先至art系列软件的官网，下载软件，在ubuntu系统安装，然后阅读相关参数设置的帮助文档，运行程序。
+ g- {1 A8 n, y: T( RGenomeABC
9 e/ }) k5 U6 e. K% E) a进入GenomeABC(http://crdd.osdd.net/raghava/genomeabc/)，输入参数，获得模拟测序结果。
编程模拟测序
. y9 _3 B8 W1 v' t. R下载安装python，并且安装biopython扩展模块，编写程序，模拟单端/双端测序。
( `$ A. g& F: K& A' {三、结果
1 ?" k: F* M9 {4 n8 U
7 X) h: G! N) b6 i1 `/ {4 W1、art_454的运行结果
* Q7 t) }/ Q# L! s8 x
9 k- q0 P; c$ e无参数art_454运行，阅读帮助文档

图表 1无参数art_454运行
对酵母基因组进行基因组单端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20.
下图为模拟单端测序，程序运行过程及结果
' j8 e% K0 Z& \; Y' d3 F
图表 2 art454单端测序

2 O) v2 E& x, ~' S9 C% w3 D图表 3 art454单端模拟结果
双端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20；MEAN_FRAG_LEN设为1500，即平均片段长度为1500；STD_DEV设为20，即长度的标准差为20
1 i. A3 V7 D. v9 k! B% ]0 @/ l下图为模拟双端测序，程序运行过程及结果
* V0 s6 R3 |8 K
! W& t# \/ o$ y3 s5 Y: t图表 4 art454双端测序
( a& `, r% u7 ~# V" W" z7 e& j
图表 5 art454双端模拟结果
2、GenomeABC
下图为设置参数页面
; D0 S4 Y7 ~; m' U
4 Y5 m  r+ q* Y1 ~4 v5 I下图为结果下载页面

# M7 H, F' d' g6 H" V. p图表 6 结果下载页面
3、编程模拟测序结果
拷贝数是这里的N值；覆盖度是m，测序深度是宏观的量，在这里与覆盖度意思相同，就是测序仪10X，20X。
1 @" Y+ b( n( x0 Y; G( f单端测序
. |4 v, R% z: ?2 q
* k7 Q- ~7 V+ ?# o% ~* d图表 7 程序模拟单端测序
5 N7 u) ]* @; n: l6 U( `8 N双端测序
/ u! e' V% Z4 N
图表 8 程序模拟双端测序
! k' n. B9 e  ~测序结果

图表 9 结果文件

" I7 |- O, i! ]  r( g因为期望片段长度是600bp，在片段长度区间200-1000bp内，所以大部分的片段都没有删除。
测序结果统计表

测序方式        基因组大小(bp)        片段长度区间 (bp)        N值        期望片段长度        克隆保留率        片段数量        Reads长度范围（bp）        Reads总数量        Reads总长度        覆盖度(m值)        理论丢失率(e-m)        覆盖率(1-e-m)
4 H. T' D4 A3 q9 p" r6 ^( A# {单端        12157kb        200-1000        10        600        0.95        107378        50-100        101968        7645.541kb        0.62889        0.53318        0.46682
% _$ y, |! Y! `% R5 P8 K单端        12157kb        200-1000        20        600        0.95        213722        50-100        202996        15227.882kb        1.25259        0.28576        0.71424
8 w7 V3 e3 ^" w4 r, v双端        12157kb        200-1000        10        600        0.95        106704        50-100        202770        15212.662kb        1.25134        0.28612        0.71388
2 g* a5 T0 h7 l双端        12157kb        200-1000        20        600        0.95        214212        50-100        407186        30534.265kb        2.51164        0.08114        0.91886
四、讨论和结论
7 G; g  f5 C  Q" b! E$ F, w+ p
程序运行方法
: z7 Q" h2 g& Y; i% A
4 K8 ^# g0 o( a1 Z5 |8 t在类的构造方法init()中，调整参数。
Averagefragmentlength为片段平均的长度；
minfragmentlength和maxfragmentlength是保留片段的范围；
cloneRetainprobability是克隆的保留率；
minreadslength和maxreadslength是测序reads的长度范围
; m- E9 D4 z% ?7 E  Q
模拟测序的诸多方法都封装成了Sequencing类，只需要创建类，并调用singlereadsequencing()和pairreadsequencing()方法，传入文件名的参数即可。

附录

from Bio import SeqIO
- a' M9 g0 y, ~9 dfrom math import exp
import random
: W6 R- i. n7 k7 o( H# P5 J5 s+ y
class Sequencing:
    # N代表拷贝份数
7 [( W7 i+ T- W8 Z    def __init__(self)
        self.fragmentList = []
        self.readsID = 1
        self.readsList = []
4 j5 W, J2 E+ b! u" K; X) G5 U        self.averagefragmentlength = 650
        self.minfragmentlength = 500
        self.maxfragmentlength = 800
        self.cloneRetainprobability = 1
  r" y3 K: }0 @        self.minreadslength = 50
        self.maxreadslength = 150
# Q6 V- b; I" Q  }( o- k        self.N = 10
( K7 I$ F2 C$ d) B  s# M! v        self.genomeLength = 0
2 U5 _7 `4 ?* L8 X( R# x        self.allreadslength = 0

    # 生成断裂点
    def generatebreakpoint(self, seqlen, averageLength):
$ g, m; z* @- N6 |        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
        breakpoint.append(seqlen)
3 b1 _8 Z" g9 }( o/ A. O2 i- E        breakpoint.append(0)
        # 把随机断裂点从小到大排序
        breakpoint.sort()
' B* R& q+ Q: V  a0 i$ T        return breakpoint

' d( w( O2 N6 Y4 N/ p) O    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
3 w& C( w5 O2 ~! D    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
        for i in range(len(breakpoint) - 1):
. d! u  L0 ~# }7 i! b7 ]4 N            fragment = seq[breakpoint:breakpoint[i + 1]]
            if maxfragmentlength > len(fragment) > minfragmentlength:
2 v: b4 N3 j1 Y* |% {% P) M/ W. F5 l                self.fragmentList.append(fragment)
: y4 ?8 x: K3 E7 r" q/ A        return self.fragmentList

8 N: ]6 {3 d- _    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
    def clonefragment(self, fragmentList, cloneRetainprobability):
        clonedfragmentList = []
        Lossprobability = [random.random() for _ in range(len(fragmentList))]
' I# @6 w4 n) |6 S        for i in range(len(fragmentList)):
            if Lossprobability <= cloneRetainprobability:
/ t4 z$ X7 {2 B; U                clonedfragmentList.append(fragmentList)
3 w  c, R: b; p        return clonedfragmentList

    # 模拟单端测序,并修改reads的ID号
    def singleread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
1 v8 b; Y5 Q: H8 T1 g! B- E            fragment.id = ""
+ J4 T. W8 V# H: x            fragment.name = ""
9 E  n  R* B; r9 d+ ^! @! d1 k5 ^            fragment.description = fragment.description[12:].split(",")[0]
1 n' `+ e6 }5 Y" l  u            fragment.description = str(self.readsID) + "." + fragment.description
            self.readsID += 1
            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])
/ }/ N0 L! d% s/ f) K, K  E
    def singlereadsequencing(self, genomedata, sequencingResult):
        for seq_record in SeqIO.parse(genomedata, "fasta"):
9 V# g5 a. t9 f; g. Z' M            seqlen = len(seq_record)
1 J0 Q5 c: z7 Y( E* j3 U6 q            self.genomeLength += seqlen
$ {( G; D8 x  z5 ?, l            for i in range(self.N):
                # 生成断裂点
$ L; d9 z/ e8 l: S5 c                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
; |" f- \7 K: l3 Y' i5 @. B                # 沿断裂点打断基因组
- i) M5 A3 p9 t2 y, w  p/ T: |                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟单端测序
- f" P  X4 ]1 D. a) Q; @        self.singleread(clonedfragmentList)
$ E& F; B! H0 S2 K- L        SeqIO.write(self.readsList, sequencingResult, "fasta")

    def pairread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
8 B1 m+ s! ], [, B) o            fragment.id = ""
            fragment.name = ""
            description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + description
9 r; n6 I1 |8 A, H) F            readslength = random.randint(self.minreadslength, self.maxreadslength)
1 ?2 v7 f1 h9 _, y2 h            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])

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

" R7 i: c* x0 e$ C7 u; C+ Q            fragmentcomplement = fragment.reverse_complement()
            fragmentcomplement.id = ""
  i$ U4 Y( w$ R  c' Q" g            fragmentcomplement.name = ""
            fragmentcomplement.description = str(self.readsID) + "." + description
: G) D& }! h* p8 a5 h            self.readsList.append(fragmentcomplement[:readslength])
9 O. M4 a1 E$ J2 l- v
7 [; I( A  a( q' e3 D! N1 _+ z            self.readsID += 1
( G: J) C6 O- v) ]
( R8 P. ^2 s7 O8 i5 u  ^/ h    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
+ R1 p+ a: D* _' T5 r            self.genomeLength += seqlen
            for i in range(self.N):
                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
2 }2 g5 L* z( d( k# i, k" J        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟双端测序
" A. c6 H. H% \# f        self.pairread(clonedfragmentList)
4 A1 N) @2 L7 Y& p- O/ ]) X        readsList_1 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 0]
3 ~8 y/ g0 \0 }# i# K        readsList_2 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 1]
        SeqIO.write(readsList_1, sequencingResult_1, "fasta")
. T, i3 S( ?3 W& V        SeqIO.write(readsList_2, sequencingResult_2, "fasta")
/ C2 G! Z  O: `
8 A1 ]. w7 ~# S9 w5 ?- _: x% b    def resultsummary(self):
& V/ ?( n: Z; @- w: E        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
7 e5 k$ X, q1 L& t5 ^6 i        print("N值:" + str(self.N))
, E3 T# c" d1 Y5 `' c( _        print("期望片段长度:" + str(self.averagefragmentlength))
1 S8 j. w( ?, f+ e. f3 Y        print("克隆保留率:" + str(self.cloneRetainprobability))
# w+ q3 g5 N* o7 Q; w: Z/ b* l        print("片段数量:" + str(len(self.fragmentList)))
" X+ O) y2 I! Q% [& L        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
& i% P, ?$ v1 h6 H        print("reads总数量:" + str(len(self.readsList)))
! C* }; j2 H; g        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
8 t6 ]2 x3 D( n1 @0 x2 g) B% t        m = self.allreadslength / self.genomeLength
6 t0 J* q( n% U) X        print("覆盖度(m值):" + str(round(m, 5)))
( X$ E! `+ |# I3 }% G, L' J        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
) F. ^$ J0 N, e% G        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
- ^$ i) V! h$ f+ C  N: I5 ]3 o# -------------------------------------------主程序-------------------------------------------
* a- A7 }, R$ t( @# 模拟单端测序
% _5 i4 y4 L% i' ZsequencingObj = Sequencing()
sequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
6 E2 ?0 O: l* A& I3 K3 G9 R: Y2 vsequencingObj.resultsummary()

# 模拟双端测序
sequencingObj = Sequencing()
" w% Y' m9 M/ m% WsequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
4 i% @  S9 O- QsequencingObj.resultsummary()
from Bio import SeqIO
  T! I1 U3 q3 A* Qfrom math import exp
import random

class Sequencing:
    # N代表拷贝份数
3 z1 Q4 }  g1 X* X    def __init__(self):
" ]* c' |  h1 y# e' X, ~/ n" c        self.fragmentList = []
        self.readsID = 1
: J4 Y' I9 }8 D! `8 a: W        self.readsList = []
        self.averagefragmentlength = 650
        self.minfragmentlength = 500
7 B% ^  s7 r! O8 ]. s5 U        self.maxfragmentlength = 800
8 B: D  Q- N; l9 v2 q        self.cloneRetainprobability = 1
        self.minreadslength = 50
% P/ M! ?0 s+ T* a4 h4 b6 m        self.maxreadslength = 150
' a1 S7 N& q8 p& u, F4 p7 N9 x        self.N = 10
/ T" n! u9 j+ A) J9 q: c        self.genomeLength = 0
9 @: d6 @2 Z9 ^0 J0 U        self.allreadslength = 0
% k) h) W+ _2 w& t" t+ [8 Q
    # 生成断裂点
    def generatebreakpoint(self, seqlen, averageLength):
        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
+ o$ x9 {7 B5 c7 ]        breakpoint.append(seqlen)
        breakpoint.append(0)
        # 把随机断裂点从小到大排序
        breakpoint.sort()
2 S: X3 [0 ~. i! T- o) e+ _        return breakpoint
% m) J/ X, g: |& U4 R
( I0 v' M5 d/ D' }# X" N- S    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
0 k* Z1 @" A9 n) }8 ]    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
4 w3 \& {4 x# F/ u# e' n        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
+ O# `* l3 r0 z7 @. A
    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
+ ^" L8 f& W6 f; P& R/ e    def clonefragment(self, fragmentList, cloneRetainprobability):
. p+ d' q( j- s( J4 B) V- S        clonedfragmentList = []
1 }9 o: x5 R3 l( a9 ]2 D        Lossprobability = [random.random() for _ in range(len(fragmentList))]
        for i in range(len(fragmentList)):
0 Z+ ]# R; X4 J1 ?  J. `8 l5 t" Z            if Lossprobability <= cloneRetainprobability:
                clonedfragmentList.append(fragmentList)
        return clonedfragmentList
' w  V2 Q  t  x; g
    # 模拟单端测序,并修改reads的ID号
! ~0 n, M$ `1 a7 Y# ]    def singleread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
            fragment.id = ""
            fragment.name = ""
            fragment.description = fragment.description[12:].split(",")[0]
. e4 Q: m  l/ g7 w            fragment.description = str(self.readsID) + "." + fragment.description
            self.readsID += 1
+ N1 A: m7 v* [; ~# D/ `6 T            readslength = random.randint(self.minreadslength, self.maxreadslength)
% N% \/ D+ C$ v            self.allreadslength += readslength
5 n1 a1 O  m, V( X            self.readsList.append(fragment[:readslength])

0 @1 o" c1 L& k( {' `    def singlereadsequencing(self, genomedata, sequencingResult):
        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
$ @% i# D! a: G7 X+ Y            self.genomeLength += seqlen
            for i in range(self.N):
: b2 v' _5 k  c                # 生成断裂点
# S% [9 ~. C7 r2 p                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
1 H" g. q' G) I, k' Y! r                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
/ t0 c6 }  l9 I- [8 m. a        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟单端测序
        self.singleread(clonedfragmentList)
        SeqIO.write(self.readsList, sequencingResult, "fasta")

8 _" v- K# `" z& ?7 x    def pairread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
            fragment.id = ""
+ ?1 J$ C0 m3 m2 [            fragment.name = ""
6 N, S* v1 ]2 Y9 e            description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + description
3 M& {/ D) C5 ]9 w            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
- S! s7 W% }8 j! D            self.readsList.append(fragment[:readslength])
* G% ?6 n: {6 k; w1 X& ?* p6 [
            readslength = random.randint(self.minreadslength, self.maxreadslength)
) ?: e" W4 H0 p0 v* l5 ^" ]% V            self.allreadslength += readslength

; I( n2 ]7 K: s4 {8 I7 G3 D- E            fragmentcomplement = fragment.reverse_complement()
* w, F8 k( l9 T8 m5 z            fragmentcomplement.id = ""
& D% Q* o' U# ~: R& j" j$ [! c            fragmentcomplement.name = ""
            fragmentcomplement.description = str(self.readsID) + "." + description
            self.readsList.append(fragmentcomplement[:readslength])
7 q) L5 n1 J5 m5 O" n
            self.readsID += 1
& J$ l8 u/ K5 F0 @% h
. X8 B* z( v. Q3 ?* W) V0 o    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):
                # 生成断裂点
4 _# R& l3 V, f2 h6 D. ]; A1 O& j                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
% p8 [9 r3 |- Y7 g7 N        # 模拟克隆时的随机丢失情况
( A5 K: G6 X( S+ N1 ?  v        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟双端测序
        self.pairread(clonedfragmentList)
( I# t7 y& _3 f; n* D1 R        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]
+ K: ]# |9 e3 ~& u" u        SeqIO.write(readsList_1, sequencingResult_1, "fasta")
        SeqIO.write(readsList_2, sequencingResult_2, "fasta")
/ ~, q% K" A7 ]) o3 m/ j
    def resultsummary(self):
        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
' x' D; y2 Y/ [/ e        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
# W/ {7 Q% f- K' K( O. x! E        print("N值:" + str(self.N))
        print("期望片段长度:" + str(self.averagefragmentlength))
        print("克隆保留率:" + str(self.cloneRetainprobability))
: Y- [( Y: |& a5 y# d- Q# M3 E; p        print("片段数量:" + str(len(self.fragmentList)))
' p- J: ?7 o  q4 [5 s        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
+ v" \7 W; X! O/ p4 @0 Z/ Y9 B        print("reads总数量:" + str(len(self.readsList)))
: G. i1 E* |/ x: o* j9 y+ H7 ~3 {        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
. h0 ^' g) h1 n$ N( e; b4 f        m = self.allreadslength / self.genomeLength
        print("覆盖度(m值):" + str(round(m, 5)))
7 D! m6 h& _  ^4 j        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
0 e: D6 H9 x7 A+ z& e# -------------------------------------------主程序-------------------------------------------
! |" T1 B! S" Z& |0 L# 模拟单端测序
sequencingObj = Sequencing()
0 R# P4 h' K( S! e8 d' FsequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
1 w% a1 d( m& M" `3 k+ _! ?+ A6 zsequencingObj.resultsummary()
' O3 {* E$ i& e6 ?
8 e, y# [. Y0 T/ ]# 模拟双端测序
, L: k. ?8 `- D' Y" J0 w; fsequencingObj = Sequencing()
sequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
' p6 n, ]! l0 C; P$ {  BsequencingObj.resultsummary()

3477959497

不错。。。。。。。。。。。。。。。。。
: ?3 E  e* E% F; T/ A  h# D