基因组测序模拟

杨利霞

基因组测序模拟
7 ^6 q. a4 M8 s! Q1 I基因组测序模拟

9 @, k- q7 c. P5 Y1 R) D) P一、摘要

通过熟悉已有的基因组测序模拟和评估程序，加深全基因组鸟枪法测序原理的理解，并且能够编写程序模拟全基因组鸟枪法测序，理解覆盖度、测序深度、拷贝数等概念，设置测序相关参数，生成单端/双端测序结果文件

二、材料和方法

1、硬件平台

处理器:Intel（R） Core（TM）i7-4710MQ CPU @ 2.50GHz
安装内存(RAM)：16.0GB
* ~  @* ~8 l& Z
2、系统平台
0 _; p% w) l. ^# O# MWindows 8.1，Ubuntu
# S( }  J$ \% Y# m5 H
7 O% B' n# s% v3、软件平台

art_454
GenomeABC http://crdd.osdd.net/raghava/genomeabc/
Python3.5
( c' d. J. [0 A# x' S5 DBiopython
8 S) \  y! t, w4、数据库资源

NCBI数据库：https://www.ncbi.nlm.nih.gov/
, q/ c" L1 t( s8 S/ a" W# U
5、研究对象
* C0 n* I' E& Z6 N: V
酵母基因组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
, x+ s4 o6 ~, P9 ?7 U
6、方法

art_454的使用
- m  N2 |, @% Q: R  w% x首先至art系列软件的官网，下载软件，在ubuntu系统安装，然后阅读相关参数设置的帮助文档，运行程序。
  x- ?5 s* f% B, ~' }GenomeABC
进入GenomeABC(http://crdd.osdd.net/raghava/genomeabc/)，输入参数，获得模拟测序结果。
6 q( J7 Y( m8 m; @! P- g( V编程模拟测序
6 {9 F% ?8 F3 v. p# F下载安装python，并且安装biopython扩展模块，编写程序，模拟单端/双端测序。
三、结果

1、art_454的运行结果

无参数art_454运行，阅读帮助文档
5 `3 M, Q' ?2 c) Z6 ]  `" B
图表 1无参数art_454运行
0 [2 l6 m) {+ H7 \对酵母基因组进行基因组单端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20.
下图为模拟单端测序，程序运行过程及结果

! E7 o% j" L4 o8 h' r; ]图表 2 art454单端测序

图表 3 art454单端模拟结果
: w' i2 y% i# l; j- J0 x双端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20；MEAN_FRAG_LEN设为1500，即平均片段长度为1500；STD_DEV设为20，即长度的标准差为20
: H0 E; ^( p3 ?下图为模拟双端测序，程序运行过程及结果

图表 4 art454双端测序
. f1 ?1 }# I3 b/ B
0 S, z+ ~* t2 t  V3 F图表 5 art454双端模拟结果
2、GenomeABC
下图为设置参数页面

下图为结果下载页面
- l2 \! P+ U2 ?" B5 P
1 F8 H+ x2 U$ k1 \* k% R图表 6 结果下载页面
3、编程模拟测序结果
拷贝数是这里的N值；覆盖度是m，测序深度是宏观的量，在这里与覆盖度意思相同，就是测序仪10X，20X。
单端测序
1 v# x# _: c4 r. N1 |
图表 7 程序模拟单端测序
双端测序
6 [- a' k; Q5 V5 T) a8 j, R% m& Z
) p5 }- K# x9 y# o7 l图表 8 程序模拟双端测序
测序结果
% j2 D- F, Y4 o8 J1 _0 d
图表 9 结果文件
3 W3 m' j5 [$ Q$ y
因为期望片段长度是600bp，在片段长度区间200-1000bp内，所以大部分的片段都没有删除。
. I3 g  C( @/ l! f! d8 U测序结果统计表

测序方式        基因组大小(bp)        片段长度区间 (bp)        N值        期望片段长度        克隆保留率        片段数量        Reads长度范围（bp）        Reads总数量        Reads总长度        覆盖度(m值)        理论丢失率(e-m)        覆盖率(1-e-m)
( p# S7 y' P1 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
双端        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
四、讨论和结论

6 i  g1 j& b1 r0 ]1 J* G! J程序运行方法
; p" G! [, ?, j) d. n
在类的构造方法init()中，调整参数。
Averagefragmentlength为片段平均的长度；
minfragmentlength和maxfragmentlength是保留片段的范围；
. ?- N4 z* k% ^2 L5 f2 P% |5 x$ f) VcloneRetainprobability是克隆的保留率；
, [0 A, s7 W. `' C. D! K5 p) o& ^minreadslength和maxreadslength是测序reads的长度范围
; V+ |+ N6 Q" S" T3 {. O
模拟测序的诸多方法都封装成了Sequencing类，只需要创建类，并调用singlereadsequencing()和pairreadsequencing()方法，传入文件名的参数即可。

' L$ D& p9 i! @% q4 U附录
4 R3 R' F! t# n2 W  N+ x. ~
from Bio import SeqIO
/ v# b4 V( f2 s2 Gfrom math import exp
3 ?# w1 S# q  M( l/ nimport random

# _  r7 A8 g; F( s. lclass Sequencing:
    # N代表拷贝份数
( N- T; N+ q, ]3 `    def __init__(self)
        self.fragmentList = []
        self.readsID = 1
1 w+ S: [7 K% D( V6 P/ O; p0 y        self.readsList = []
# E1 ^( a2 p( ^) v% \8 j/ v+ t1 {, b* f0 {        self.averagefragmentlength = 650
! L& \, Y2 r3 A. ~2 l2 z        self.minfragmentlength = 500
        self.maxfragmentlength = 800
6 T4 Z8 _" x- L! s        self.cloneRetainprobability = 1
3 c$ v2 K$ U: ?  R) c) S+ f8 H        self.minreadslength = 50
        self.maxreadslength = 150
0 e9 }6 x# D0 T$ V! H        self.N = 10
; K/ w: u$ L3 W: a/ s; h        self.genomeLength = 0
        self.allreadslength = 0

    # 生成断裂点
- U  j' E! a$ ?$ ~: D! p5 f4 `6 x    def generatebreakpoint(self, seqlen, averageLength):
3 Y2 v, `6 w5 B+ {        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
        breakpoint.append(seqlen)
        breakpoint.append(0)
        # 把随机断裂点从小到大排序
. [" R3 c( l6 z7 L: A        breakpoint.sort()
2 Q9 ?" Y/ A* \' R% l/ m% O        return breakpoint
8 ~" W! d2 R1 ^8 W1 s3 @
    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
6 k$ J7 z& M( `/ }0 h6 a1 \        for i in range(len(breakpoint) - 1):
% Y1 Q. ^' ]& A  E, K            fragment = seq[breakpoint:breakpoint[i + 1]]
            if maxfragmentlength > len(fragment) > minfragmentlength:
                self.fragmentList.append(fragment)
! I6 H# a" n( ?        return self.fragmentList

    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
    def clonefragment(self, fragmentList, cloneRetainprobability):
        clonedfragmentList = []
        Lossprobability = [random.random() for _ in range(len(fragmentList))]
0 W4 ^' D: `% g/ r& O" {3 K        for i in range(len(fragmentList)):
            if Lossprobability <= cloneRetainprobability:
9 N# }) t1 a& j                clonedfragmentList.append(fragmentList)
        return clonedfragmentList
) g6 i& m9 q9 u; A
    # 模拟单端测序,并修改reads的ID号
    def singleread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
6 B9 R  x- W: r) y' m) @' [            fragment.id = ""
8 _$ |8 E" l- p3 a& o" u) T            fragment.name = ""
            fragment.description = fragment.description[12:].split(",")[0]
: P# X: }* a9 N. t) T5 T, i            fragment.description = str(self.readsID) + "." + fragment.description
            self.readsID += 1
2 K$ \* q5 D& [# {1 G5 |  B            readslength = random.randint(self.minreadslength, self.maxreadslength)
  d7 W+ T8 O8 F9 i+ q2 n0 C: ], C            self.allreadslength += readslength
; Y( a# S  O, @+ j+ z2 {# `: X' u$ S            self.readsList.append(fragment[:readslength])

5 d/ w7 s0 I% y, `5 ?    def singlereadsequencing(self, genomedata, sequencingResult):
$ T3 Y! b1 [. }( x5 b        for seq_record in SeqIO.parse(genomedata, "fasta"):
* G. j' Q! |, [& J5 M* M* ^( Z            seqlen = len(seq_record)
            self.genomeLength += seqlen
+ a0 m7 K% h- ~1 v  K1 E  I% [0 e            for i in range(self.N):
                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
" s- o( W& p0 Q  ^, e                # 沿断裂点打断基因组
( \9 W: L/ v  @' Z0 z: ^                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
2 k) k0 R7 ]" Y2 V" z* i" s( ]+ d        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟单端测序
        self.singleread(clonedfragmentList)
        SeqIO.write(self.readsList, sequencingResult, "fasta")

! w, V  f- Y: h$ B6 {) s, r    def pairread(self, clonedfragmentList):
% D6 T$ j9 P' X6 P) S" u; T' }2 Y        for fragment in clonedfragmentList:
4 q1 K9 L& P7 ~* K% {1 K! w6 C            fragment.id = ""
" h( B% a- }  z; G0 D            fragment.name = ""
            description = fragment.description[12:].split(",")[0]
9 R( O2 B9 O& T& I9 ?            fragment.description = str(self.readsID) + "." + description
  A* X  A* C. R2 w            readslength = random.randint(self.minreadslength, self.maxreadslength)
1 x, o3 ~( H! x- S            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])

3 k7 ]$ _6 {: J2 Z& [4 P            readslength = random.randint(self.minreadslength, self.maxreadslength)
+ C7 x/ i: V% n( j            self.allreadslength += readslength
3 l5 Y. d+ [+ J+ y- O/ j
+ ~: _: V( s$ C1 m+ c            fragmentcomplement = fragment.reverse_complement()
            fragmentcomplement.id = ""
9 W$ _0 y* I( `& t: W- o# F            fragmentcomplement.name = ""
            fragmentcomplement.description = str(self.readsID) + "." + description
% A) s# M4 ^5 _8 ?; d5 V$ ~            self.readsList.append(fragmentcomplement[:readslength])

, o, I# G4 Y  ~4 `: L& U! Z" o            self.readsID += 1
' j1 \) o# R1 R( U. l& l  f) g
! O( p# ], d2 b1 ~: D8 l    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
        for seq_record in SeqIO.parse(genomedata, "fasta"):
, M, \( k& z1 Y. a3 ~2 t5 ]2 ~9 k/ X            seqlen = len(seq_record)
            self.genomeLength += seqlen
( R* Y! O' C$ y4 O. x" \! M            for i in range(self.N):
; ^, ^5 M2 b' e( ~2 Y                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
        # 模拟克隆时的随机丢失情况
        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]
        SeqIO.write(readsList_1, sequencingResult_1, "fasta")
  t8 O+ L% A( F; G4 P( u        SeqIO.write(readsList_2, sequencingResult_2, "fasta")

: G- C% F: g3 r# G+ q6 H    def resultsummary(self):
/ X' A- [' n5 m- Q        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
7 A: |9 L( u) B" _! f: n        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
# _# K2 e$ X. o$ S        print("N值:" + str(self.N))
! M) K" K8 r9 m4 ?5 T8 h* |        print("期望片段长度:" + str(self.averagefragmentlength))
        print("克隆保留率:" + str(self.cloneRetainprobability))
        print("片段数量:" + str(len(self.fragmentList)))
  z& Y& E4 _! p7 J3 T; Y        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
        print("reads总数量:" + str(len(self.readsList)))
        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
7 o* t' h3 |$ M0 ^6 s        m = self.allreadslength / self.genomeLength
        print("覆盖度(m值):" + str(round(m, 5)))
, d" m+ p; @) X" I6 H" \1 `        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
3 ]9 Y% \1 J; _5 N1 ~# -------------------------------------------主程序-------------------------------------------
8 i# a7 X7 o" S2 ~) h% g# 模拟单端测序
sequencingObj = Sequencing()
2 o: |- l' m* y! i  A! e& ]sequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
# v/ O; K' d% ~3 Q+ DsequencingObj.resultsummary()

; d5 W5 Q; n0 {# 模拟双端测序
sequencingObj = Sequencing()
sequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
sequencingObj.resultsummary()
from Bio import SeqIO
; X6 l5 _/ Q* n9 n9 U" afrom math import exp
import random
: x4 o# k$ ^" e/ \: E& K
class Sequencing:
. u6 _. Z: D! g1 l" r9 H    # N代表拷贝份数
" K# x' ]6 K) [* C! W    def __init__(self):
6 F. W" q+ o; C/ |$ d        self.fragmentList = []
6 J+ U% Q# C! s; l$ A: p        self.readsID = 1
) z6 F$ u2 l& z4 O, Y2 u        self.readsList = []
        self.averagefragmentlength = 650
        self.minfragmentlength = 500
, z- @/ w/ w8 u1 }; |! r4 q2 g        self.maxfragmentlength = 800
9 Z8 e& h& F& f+ N8 e, R        self.cloneRetainprobability = 1
" O" j7 u1 E7 p2 |* C: H- {. _  ]# c        self.minreadslength = 50
        self.maxreadslength = 150
        self.N = 10
        self.genomeLength = 0
9 u1 k' [% ]# u/ s2 x        self.allreadslength = 0

' `0 u+ v4 P1 V; L    # 生成断裂点
. ]. H# U, u0 \    def generatebreakpoint(self, seqlen, averageLength):
; b8 y( ]! e  r7 Z        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
        breakpoint.append(seqlen)
  I; K' _. N# F2 s& A8 ?/ O* e& C; Q        breakpoint.append(0)
0 g3 s0 w9 e) q9 Y* l0 ?* M        # 把随机断裂点从小到大排序
        breakpoint.sort()
        return breakpoint
' m& o2 p  x7 T' \% |  O( Z) B
    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
1 Z: k5 @4 @. _5 V- G: u! y# }        for i in range(len(breakpoint) - 1):
4 K- `) t; _: q1 z4 i6 k6 f            fragment = seq[breakpoint:breakpoint[i + 1]]
            if maxfragmentlength > len(fragment) > minfragmentlength:
! J+ p2 v2 M! Y% p% p2 G                self.fragmentList.append(fragment)
8 W3 K* f( S2 }: D! X1 l+ {        return self.fragmentList
, Q/ e) f, S( {. x& E* X
6 N* `# B8 s& V    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
    def clonefragment(self, fragmentList, cloneRetainprobability):
$ B/ ?. F& A2 R- \        clonedfragmentList = []
        Lossprobability = [random.random() for _ in range(len(fragmentList))]
        for i in range(len(fragmentList)):
7 G9 }1 b! C+ s: U0 s! e            if Lossprobability <= cloneRetainprobability:
. H# v8 Z# Y6 K# b# J; P                clonedfragmentList.append(fragmentList)
        return clonedfragmentList

    # 模拟单端测序,并修改reads的ID号
  p" |4 S; M  F) C6 w* R    def singleread(self, clonedfragmentList):
6 p9 w( j% ?" [        for fragment in clonedfragmentList:
5 d  Z9 w# K% H& x& o/ g            fragment.id = ""
            fragment.name = ""
            fragment.description = fragment.description[12:].split(",")[0]
. _0 D- M. C7 f/ i2 B            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])
. ~- V, v/ n) x( J/ @- x
    def singlereadsequencing(self, genomedata, sequencingResult):
        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
            self.genomeLength += seqlen
" O6 [- o, a: o& I+ `            for i in range(self.N):
                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
' u3 F: z8 s; F5 Q4 l. T$ l$ Y2 m% d                # 沿断裂点打断基因组
7 J. o& y; t4 K5 `; Z                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
& A2 x* ]9 {1 D  m/ S7 q$ L" P        # 模拟克隆时的随机丢失情况
; x4 U% N9 B2 H- C# W        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟单端测序
  {2 `* o7 k. ?        self.singleread(clonedfragmentList)
        SeqIO.write(self.readsList, sequencingResult, "fasta")
( w: ^3 A5 f( p# b
: y) f& I" K- K) u    def pairread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
            fragment.id = ""
5 c# m/ k9 c7 A% _; e            fragment.name = ""
            description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + description
$ b/ X. S8 T( z7 F' _            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])

            readslength = random.randint(self.minreadslength, self.maxreadslength)
6 a5 @- }3 B! t- X( a: a2 ~: \            self.allreadslength += readslength
4 J- B3 @$ o. x; p9 B# h
+ E. K+ [" U1 i( y            fragmentcomplement = fragment.reverse_complement()
            fragmentcomplement.id = ""
            fragmentcomplement.name = ""
            fragmentcomplement.description = str(self.readsID) + "." + description
            self.readsList.append(fragmentcomplement[:readslength])
0 y* [1 S! V6 Q- d8 X7 n# K2 y
. n- q/ m. h/ g3 t            self.readsID += 1

    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
8 |6 \* {; y% [. [! x        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
% V0 w- D  |& @6 I            self.genomeLength += seqlen
            for i in range(self.N):
( j9 c# C3 [+ b9 B, g                # 生成断裂点
. o3 V3 j1 i4 K2 _                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
- c, D9 T2 u' ^        # 模拟克隆时的随机丢失情况
$ @- h0 H& g( N& S0 ^* W        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟双端测序
3 |2 x' ~. \2 y% _  `        self.pairread(clonedfragmentList)
( n* K8 X7 d1 K5 J        readsList_1 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 0]
$ T* b. n* j3 \, s4 Z0 t' |        readsList_2 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 1]
        SeqIO.write(readsList_1, sequencingResult_1, "fasta")
/ G* N. @6 o1 f( l        SeqIO.write(readsList_2, sequencingResult_2, "fasta")
! X2 L8 \) G# P% D; N. c4 e
    def resultsummary(self):
4 T7 N% ?* r- R( p8 m* r        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
% Q' E6 Z6 _0 ^' M  i        print("N值:" + str(self.N))
  N% S# X5 w' l) d        print("期望片段长度:" + str(self.averagefragmentlength))
" [4 ?3 G5 v% t/ X. z, n        print("克隆保留率:" + str(self.cloneRetainprobability))
* B% j( X9 Q$ |  K        print("片段数量:" + str(len(self.fragmentList)))
5 L* L* T, C7 q* J1 H2 J7 P        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
" `- G. I- h) ~9 i# v        print("reads总数量:" + str(len(self.readsList)))
8 j9 y- }, [1 s+ Q& j1 P        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
6 I8 Q4 a. l9 \% s% L        m = self.allreadslength / self.genomeLength
! r, j6 [  c4 y4 h& \3 Y( T        print("覆盖度(m值):" + str(round(m, 5)))
        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
" i$ z. u$ i& H$ S! i5 H+ J1 s        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
# -------------------------------------------主程序-------------------------------------------
# 模拟单端测序
sequencingObj = Sequencing()
3 v& ?7 R8 E2 @" gsequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
! `8 s8 M- K; B* O0 rsequencingObj.resultsummary()
+ f, L+ F+ y% Q+ S1 M% w5 u
2 e' a5 e+ |( l2 B! k/ M: L) B# 模拟双端测序
sequencingObj = Sequencing()
sequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
) n8 v  g2 H# }3 O% P% B6 ssequencingObj.resultsummary()
0 j5 ]0 g, m0 e% T' n
# w- P% Q' n  X) ]( S/ S0 j

% P. i/ [9 L" f9 c0 i3 N" t

3477959497

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