基因组测序模拟

杨利霞

基因组测序模拟
: z9 M- R+ s8 D. r6 ^- I0 N基因组测序模拟

* [0 \- L) F, V( n( g, I5 |一、摘要
  ~; x0 }" }6 F1 u  d
& y- h  p* \0 O/ K& x; m- V6 c通过熟悉已有的基因组测序模拟和评估程序，加深全基因组鸟枪法测序原理的理解，并且能够编写程序模拟全基因组鸟枪法测序，理解覆盖度、测序深度、拷贝数等概念，设置测序相关参数，生成单端/双端测序结果文件

: {! P* Y% h/ F! g$ F9 k二、材料和方法

1、硬件平台

处理器:Intel（R） Core（TM）i7-4710MQ CPU @ 2.50GHz
安装内存(RAM)：16.0GB
, D8 e3 \' y2 l, ~
* i& L$ M0 w+ M4 Y3 B% K. Q2、系统平台
Windows 8.1，Ubuntu
/ c7 g+ r, k% U
3、软件平台

art_454
GenomeABC http://crdd.osdd.net/raghava/genomeabc/
0 s6 I5 @9 [; Y$ Q' |1 dPython3.5
Biopython
4、数据库资源
, F# @$ C% m8 i* d$ B3 R, K
NCBI数据库：https://www.ncbi.nlm.nih.gov/
1 j- m# u: q9 o/ K/ b& A
5、研究对象
# x# @/ ~/ a$ x6 E  j& |& i$ @
酵母基因组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
4 j$ S% h" u: D0 o: S, ?
6、方法
% l: f" p) u. E4 ~6 k. k
) V6 {  n$ V4 y, ~art_454的使用
首先至art系列软件的官网，下载软件，在ubuntu系统安装，然后阅读相关参数设置的帮助文档，运行程序。
GenomeABC
# G, H/ E# E# E0 n. s进入GenomeABC(http://crdd.osdd.net/raghava/genomeabc/)，输入参数，获得模拟测序结果。
3 Z# F' O; v2 t6 R, b- d4 ~编程模拟测序
下载安装python，并且安装biopython扩展模块，编写程序，模拟单端/双端测序。
三、结果

0 t) }0 d' t0 j3 \) W' g& E1、art_454的运行结果

无参数art_454运行，阅读帮助文档
8 t7 S% k0 `6 k) U+ P" O
图表 1无参数art_454运行
" p/ h2 l) {0 i3 `$ k对酵母基因组进行基因组单端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20.
下图为模拟单端测序，程序运行过程及结果

- g9 T% ^. M; s5 V图表 2 art454单端测序
- B$ e* _; h& P
图表 3 art454单端模拟结果
/ T/ }4 {4 c  y& q双端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20；MEAN_FRAG_LEN设为1500，即平均片段长度为1500；STD_DEV设为20，即长度的标准差为20
下图为模拟双端测序，程序运行过程及结果
9 @6 N: o, P, V( m1 a  c/ R2 [
& {3 Q0 m4 q+ u+ i& X+ P, w% D$ E5 p图表 4 art454双端测序

图表 5 art454双端模拟结果
2、GenomeABC
下图为设置参数页面
: J, }5 }6 J: a5 t& M
6 J; O# P/ s& M( D; G下图为结果下载页面
4 e* J1 x1 Q% K$ L) r! W3 L" k7 r
图表 6 结果下载页面
3、编程模拟测序结果
% a$ ^: O% j; X拷贝数是这里的N值；覆盖度是m，测序深度是宏观的量，在这里与覆盖度意思相同，就是测序仪10X，20X。
单端测序

图表 7 程序模拟单端测序
1 z. A' F, A2 O/ m( v双端测序

图表 8 程序模拟双端测序
测序结果

图表 9 结果文件
% J7 P8 U; R# \* h2 J
4 {, g, X* k. j$ S3 P( J因为期望片段长度是600bp，在片段长度区间200-1000bp内，所以大部分的片段都没有删除。
测序结果统计表
+ E- `' n9 O- G/ k8 W6 n' _& }
测序方式        基因组大小(bp)        片段长度区间 (bp)        N值        期望片段长度        克隆保留率        片段数量        Reads长度范围（bp）        Reads总数量        Reads总长度        覆盖度(m值)        理论丢失率(e-m)        覆盖率(1-e-m)
7 g% X8 t9 v, q# i单端        12157kb        200-1000        10        600        0.95        107378        50-100        101968        7645.541kb        0.62889        0.53318        0.46682
3 T/ @: B& }, C8 C. E单端        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
四、讨论和结论
- c' R$ i# s8 @8 W
+ {1 N) ^9 C& ~$ }2 I程序运行方法
6 y8 Y- w! Z7 n1 A4 R, L
. x( A; j( {0 k% M  r5 v  b在类的构造方法init()中，调整参数。
# h3 y$ k2 E4 t* ]Averagefragmentlength为片段平均的长度；
minfragmentlength和maxfragmentlength是保留片段的范围；
/ ?: ^7 |9 m6 O3 r7 QcloneRetainprobability是克隆的保留率；
minreadslength和maxreadslength是测序reads的长度范围
: ]! \9 A. w+ i9 Q" l! r$ u1 `- g
+ ~0 L2 Y' D. \& z, ^模拟测序的诸多方法都封装成了Sequencing类，只需要创建类，并调用singlereadsequencing()和pairreadsequencing()方法，传入文件名的参数即可。

! T  ?" K: z1 q2 I* @  N. R附录
) ^+ {# K( |, O* y
from Bio import SeqIO
from math import exp
import random

class Sequencing:
    # N代表拷贝份数
) x, I: d! B# w6 F    def __init__(self)
        self.fragmentList = []
) h* @5 l) h- y7 m7 E9 t. U8 ?        self.readsID = 1
! A) V3 V# _2 X9 U, c+ a( V        self.readsList = []
' P* p7 ], Q+ u. r# F5 y+ d/ i        self.averagefragmentlength = 650
4 `) C1 X! X5 n/ K& i9 r3 h! F        self.minfragmentlength = 500
        self.maxfragmentlength = 800
/ |0 y6 |1 q; Z8 g0 p' w        self.cloneRetainprobability = 1
        self.minreadslength = 50
        self.maxreadslength = 150
8 X% m, Q) J9 b2 J        self.N = 10
        self.genomeLength = 0
. T; C5 @6 G$ O        self.allreadslength = 0

7 @( p% @# u8 V2 W! I! S    # 生成断裂点
    def generatebreakpoint(self, seqlen, averageLength):
        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
9 B  {; C2 t5 e# K+ g        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
        breakpoint.append(seqlen)
        breakpoint.append(0)
8 F2 U9 N/ I1 d        # 把随机断裂点从小到大排序
' b9 T, n) h- s* O7 z# |        breakpoint.sort()
        return breakpoint

. }0 ^0 K  B" p0 ^  a1 t+ P    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
( o" ^9 Q$ Y* \7 u3 f, h    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
! k8 N8 j; ?3 {+ ]' {7 r( R        for i in range(len(breakpoint) - 1):
+ H; I$ G( t3 y! d7 p  d6 |            fragment = seq[breakpoint:breakpoint[i + 1]]
# Q' z- Q- `' N' L            if maxfragmentlength > len(fragment) > minfragmentlength:
( u$ [  X- |7 Q) t" I4 l1 a& _3 l8 p                self.fragmentList.append(fragment)
        return self.fragmentList

    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
    def clonefragment(self, fragmentList, cloneRetainprobability):
" h: w. D0 x6 P, G6 U        clonedfragmentList = []
        Lossprobability = [random.random() for _ in range(len(fragmentList))]
        for i in range(len(fragmentList)):
8 j3 i/ `3 i, v/ M5 `) }: D            if Lossprobability <= cloneRetainprobability:
9 I# T5 [% F1 t( V& T1 H  k                clonedfragmentList.append(fragmentList)
* R0 [" M# L$ @1 O% C. [        return clonedfragmentList
. v" C5 G4 f9 o/ k* R
+ }" T3 g* ]. d6 m' \    # 模拟单端测序,并修改reads的ID号
. `+ m# D5 k% S' y: D5 b    def singleread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
, b' t4 `+ ~" A1 R# R            fragment.id = ""
# \4 P! J  D& e6 G; r            fragment.name = ""
            fragment.description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + fragment.description
% `/ D' i7 q" O9 L2 Q! O$ H4 e' g            self.readsID += 1
2 ]. s: J; F0 l4 ~            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
/ B- d' j* W7 i1 G4 u            self.readsList.append(fragment[:readslength])
- D8 s5 L( T1 p
- l- U/ {' x0 `* h' m1 x    def singlereadsequencing(self, genomedata, sequencingResult):
        for seq_record in SeqIO.parse(genomedata, "fasta"):
$ V2 _! V4 I8 A            seqlen = len(seq_record)
            self.genomeLength += seqlen
            for i in range(self.N):
                # 生成断裂点
: l0 d5 {8 f3 c- M2 P: j( ]8 {0 X                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
- Y7 _2 |7 j0 I% }                # 沿断裂点打断基因组
# W8 {; P5 {* ^/ q" o                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
1 b( ?, q3 N% G' v: Y5 g" s  Q; k        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
! a% ?9 v3 M; l5 `) l4 l        # 模拟单端测序
1 g# a( S  a- K. s' Y6 P. N        self.singleread(clonedfragmentList)
; ~3 `3 Q+ m/ ]1 A. C7 n        SeqIO.write(self.readsList, sequencingResult, "fasta")

    def pairread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
) _3 X+ a9 j; _# s% S2 ~- ]+ W) C7 f- K            fragment.id = ""
            fragment.name = ""
            description = fragment.description[12:].split(",")[0]
5 _4 K: F+ O* j8 `. M            fragment.description = str(self.readsID) + "." + description
            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])

0 z  t6 U% q- z" J% {" ~            readslength = random.randint(self.minreadslength, self.maxreadslength)
! u1 l0 T' `' m/ s, r1 B            self.allreadslength += readslength

            fragmentcomplement = fragment.reverse_complement()
            fragmentcomplement.id = ""
            fragmentcomplement.name = ""
            fragmentcomplement.description = str(self.readsID) + "." + description
2 W4 D) O. U" u            self.readsList.append(fragmentcomplement[:readslength])

, X8 N! L: B+ O! }2 g            self.readsID += 1
* G8 [: ^+ A. ?- E
    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
7 ?4 C& m8 l2 [8 a" z            self.genomeLength += seqlen
            for i in range(self.N):
                # 生成断裂点
9 m, v* f( `3 ~$ j! |1 S; X4 N& [                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
6 f; u1 N5 \: d) D                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
" ~8 S/ T6 k' n2 o. o( \/ v% s        # 模拟克隆时的随机丢失情况
/ C+ ^. F2 E( }        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
* |# Y* {4 K' r' m        # 模拟双端测序
& Y# p$ g1 J8 d) O% L- Y        self.pairread(clonedfragmentList)
/ K+ x9 N, }; K) }+ ?* U! R" a, s        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")
+ c9 z5 M: a  W& u3 Z1 z, l* ?: ~        SeqIO.write(readsList_2, sequencingResult_2, "fasta")

    def resultsummary(self):
        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
7 o8 Q, ~  n$ T* z        print("N值:" + str(self.N))
  J6 b+ r! }$ o        print("期望片段长度:" + str(self.averagefragmentlength))
        print("克隆保留率:" + str(self.cloneRetainprobability))
        print("片段数量:" + str(len(self.fragmentList)))
" `+ [, e; ?6 e, l        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
        print("reads总数量:" + str(len(self.readsList)))
& |) {; m/ @0 O. m& ?+ C2 |        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
$ G; k5 |, O! o; I        m = self.allreadslength / self.genomeLength
% j8 B9 \8 w/ m( c7 M        print("覆盖度(m值):" + str(round(m, 5)))
' g/ n( R- o5 D  t) L        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
# -------------------------------------------主程序-------------------------------------------
" q) k% {: r' J. f  t8 V" y# 模拟单端测序
9 O% g$ t; b; h! }sequencingObj = Sequencing()
+ E7 p- c; ]4 ^sequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
8 B- J* l& F7 x; k6 p& b$ F5 ZsequencingObj.resultsummary()

  f3 ]7 N# D; I) R4 h0 k4 `+ M# 模拟双端测序
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
0 ~* O) N$ K" U( r; }from math import exp
/ q  a) [1 b3 ~5 x, o& cimport random

; p- N5 S3 D4 C  ~1 z) a# T7 Aclass Sequencing:
( G8 O3 X9 V& j; w* i; [    # N代表拷贝份数
    def __init__(self):
! r  y, d. O% Z  J5 D" z$ x6 d; g$ z$ c        self.fragmentList = []
5 J  U3 T! z7 k8 ?4 @( @        self.readsID = 1
        self.readsList = []
$ ~8 Q- N: ?  x9 F: T        self.averagefragmentlength = 650
        self.minfragmentlength = 500
1 i6 F1 D) Y4 g        self.maxfragmentlength = 800
1 C( ]. `0 E3 e8 Q( V' p        self.cloneRetainprobability = 1
        self.minreadslength = 50
5 V0 \  ^/ c6 z" ?0 c        self.maxreadslength = 150
        self.N = 10
        self.genomeLength = 0
+ E. i8 s, @% s$ y( P; r  p5 Q        self.allreadslength = 0

    # 生成断裂点
    def generatebreakpoint(self, seqlen, averageLength):
9 X0 U: A7 w# Z. C' _& H, [        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
( r3 `* _' N* D) _. d& c        breakpoint.append(seqlen)
        breakpoint.append(0)
        # 把随机断裂点从小到大排序
        breakpoint.sort()
        return breakpoint
0 J0 c, f% m! d, Y; G6 k
/ s0 C; h2 I3 G) ]3 B    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
        for i in range(len(breakpoint) - 1):
            fragment = seq[breakpoint:breakpoint[i + 1]]
            if maxfragmentlength > len(fragment) > minfragmentlength:
4 ]( f' c4 H! p" S, e                self.fragmentList.append(fragment)
( E2 @6 @8 H8 v0 s3 F  A        return self.fragmentList

    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
7 J, v# F% Z, j8 r& `* x/ A  [    def clonefragment(self, fragmentList, cloneRetainprobability):
        clonedfragmentList = []
; B1 H* i1 I% r6 ]6 p+ R        Lossprobability = [random.random() for _ in range(len(fragmentList))]
        for i in range(len(fragmentList)):
            if Lossprobability <= cloneRetainprobability:
                clonedfragmentList.append(fragmentList)
        return clonedfragmentList

    # 模拟单端测序,并修改reads的ID号
) R0 K" Y! t' L% O. h9 ^    def singleread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
            fragment.id = ""
6 j: V: \; J, l$ t% i, ]* J            fragment.name = ""
' `% ]: f& R$ h" `1 E            fragment.description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + fragment.description
3 z5 k# J  h& G9 `5 m9 n3 `2 U8 @            self.readsID += 1
5 t% V0 Z% v! m: B) ]            readslength = random.randint(self.minreadslength, self.maxreadslength)
: Q% f# G: N4 C0 Z6 G            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])
5 n6 a: O$ y# q- ^# r# e
    def singlereadsequencing(self, genomedata, sequencingResult):
$ K' P4 q1 @/ L+ n: H        for seq_record in SeqIO.parse(genomedata, "fasta"):
4 ]7 k/ `0 Z. K! B( s) n% E            seqlen = len(seq_record)
1 y* n8 S, B5 P" i- M. k            self.genomeLength += seqlen
            for i in range(self.N):
                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
# v  `1 `2 F- K3 _                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
6 q" q) g: _! h  e0 b1 t        # 模拟克隆时的随机丢失情况
, t" _: C2 n4 c. M8 @        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
7 ~5 ?0 y# @2 ^' p        # 模拟单端测序
        self.singleread(clonedfragmentList)
        SeqIO.write(self.readsList, sequencingResult, "fasta")
" D5 P5 b4 Z) D7 P. R
    def pairread(self, clonedfragmentList):
& Q1 Y* Z$ a  j; V: f6 l" S        for fragment in clonedfragmentList:
            fragment.id = ""
- h$ j7 a+ x3 A# P/ R; S" e            fragment.name = ""
            description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + description
            readslength = random.randint(self.minreadslength, self.maxreadslength)
# |6 A- [2 F9 {! d2 l5 ]" H$ w            self.allreadslength += readslength
/ {2 g3 L4 }7 r0 T& {& z, E' W! }            self.readsList.append(fragment[:readslength])
  n( {0 u$ W7 A# x; X8 D2 j7 G
            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
8 r7 I5 X! U8 {% V* j
% \6 I; ^- v( Q) n            fragmentcomplement = fragment.reverse_complement()
9 I+ e9 v5 t! C' i. t+ A3 {            fragmentcomplement.id = ""
. _7 n9 h# t9 g1 e            fragmentcomplement.name = ""
" I/ A  V, Q! r            fragmentcomplement.description = str(self.readsID) + "." + description
            self.readsList.append(fragmentcomplement[:readslength])

            self.readsID += 1
& K0 J+ l8 @( t# l, N  z
/ o4 E6 `/ [$ F$ ]( F& N    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
9 X& Z, g, z9 L        for seq_record in SeqIO.parse(genomedata, "fasta"):
% ^( \; _4 u, J& F: j4 o" x            seqlen = len(seq_record)
            self.genomeLength += seqlen
# l$ i( \) t6 |/ i) q  S  S            for i in range(self.N):
/ A' b1 c, T, u( f% O                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
        # 模拟克隆时的随机丢失情况
8 w1 b3 {4 t3 S- B+ v        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]
5 E/ C, s1 D, S8 X        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")
, M9 q( o5 L& g  v( l, Y
" o, X0 y  f7 m& l& h* @4 B6 j; Y    def resultsummary(self):
        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
2 d' L5 B. w. E; M+ m  W        print("N值:" + str(self.N))
        print("期望片段长度:" + str(self.averagefragmentlength))
        print("克隆保留率:" + str(self.cloneRetainprobability))
        print("片段数量:" + str(len(self.fragmentList)))
  u. m  Z+ T) _9 d/ ~2 x( H        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
        print("reads总数量:" + str(len(self.readsList)))
        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
  ]' F, k. x) i  y+ t' _        m = self.allreadslength / self.genomeLength
        print("覆盖度(m值):" + str(round(m, 5)))
7 O/ q! O# e4 d; H# ?3 B$ E        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
, I- v" _7 g$ E% C        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
% i# P& G3 ^( T' c# -------------------------------------------主程序-------------------------------------------
# 模拟单端测序
1 h" V  K% B' \, o! S% _sequencingObj = Sequencing()
# V2 r3 P5 H" csequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
0 `6 o0 R1 Q& L( H. \! |sequencingObj.resultsummary()
8 ?8 G: l3 s6 _1 t8 E: Y
; Y$ m# W& S+ y8 d  j# 模拟双端测序
sequencingObj = Sequencing()
( `' w$ {) f, l0 I% xsequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
1 j: O3 [8 O$ S. |( s6 o+ v; ssequencingObj.resultsummary()

2 s* I' }+ r4 T2 ^$ }0 u

; U+ g, m& N& O# t2 g3 W9 I
9 I5 n! H  m, Q5 }$ j2 O" i

3477959497

不错。。。。。。。。。。。。。。。。。
$ d- V' o1 t: Y1 L. m5 K$ W