基因组测序模拟

杨利霞

基因组测序模拟
基因组测序模拟
7 E; j) D- Z! h$ G6 e
, d* a4 i% Z. d1 t一、摘要

通过熟悉已有的基因组测序模拟和评估程序，加深全基因组鸟枪法测序原理的理解，并且能够编写程序模拟全基因组鸟枪法测序，理解覆盖度、测序深度、拷贝数等概念，设置测序相关参数，生成单端/双端测序结果文件
/ t# n% s( _9 J1 m: T+ K
0 V, E/ E2 h) G( b二、材料和方法

1、硬件平台

处理器:Intel（R） Core（TM）i7-4710MQ CPU @ 2.50GHz
3 `( h4 d: w! N. }! T$ x安装内存(RAM)：16.0GB
0 B- E7 u* ^' M/ T- K
2、系统平台
9 {& E; t$ n' Z  h- ?8 _Windows 8.1，Ubuntu
: @# ~4 _2 Q2 @! A2 b1 c5 P' ^
5 H5 g/ \+ x, o! [3、软件平台

art_454
( x1 Z; K+ A1 xGenomeABC http://crdd.osdd.net/raghava/genomeabc/
" Z; `$ W9 C# CPython3.5
Biopython
4、数据库资源
! f' S( y% _1 z  }& m2 W
# Y' Q& Q! R5 ^9 D- V# C& v+ cNCBI数据库：https://www.ncbi.nlm.nih.gov/

/ [! X/ r5 i) Q$ r/ ]5 [5、研究对象

4 b# P; G! `8 M6 g0 v: e* X( y2 P酵母基因组Saccharomyces cerevisiae S288c (assembly R64)
" U, |: C  W3 w! W( a# v" q: W9 X' yftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_genomic.fna.gz
6 x. ?2 }+ K( p1 z* S% x" z/ X
6、方法
5 o$ j+ h8 g0 g& q% }4 A2 t
art_454的使用
首先至art系列软件的官网，下载软件，在ubuntu系统安装，然后阅读相关参数设置的帮助文档，运行程序。
' I6 {4 D$ I  m3 Q) X! u0 r; K5 t. HGenomeABC
' o7 L4 L4 u$ Q% G# [* T4 c进入GenomeABC(http://crdd.osdd.net/raghava/genomeabc/)，输入参数，获得模拟测序结果。
编程模拟测序
: W/ k& s9 z% i9 L下载安装python，并且安装biopython扩展模块，编写程序，模拟单端/双端测序。
三、结果

! }) x( i: H* g; \3 Q) V" {1、art_454的运行结果
8 C" t4 U9 q/ C: d
( V8 H3 h1 D+ ~3 E4 a1 D无参数art_454运行，阅读帮助文档

图表 1无参数art_454运行
对酵母基因组进行基因组单端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20.
下图为模拟单端测序，程序运行过程及结果

4 a' o% N0 `; i+ ]( G图表 2 art454单端测序
1 K8 N" k9 a6 A; I5 I; w
2 y* M* J* d- L9 S2 ?  y5 w3 t图表 3 art454单端模拟结果
双端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20；MEAN_FRAG_LEN设为1500，即平均片段长度为1500；STD_DEV设为20，即长度的标准差为20
下图为模拟双端测序，程序运行过程及结果

! t# W, w* A8 v# C/ d图表 4 art454双端测序

图表 5 art454双端模拟结果
- D4 T/ U  _  Y1 @; e2、GenomeABC
下图为设置参数页面
4 D8 M+ e; y0 s% A& S+ w
. b; \( D2 R' i0 }2 _2 L6 }下图为结果下载页面

. ?3 c& Z) i. \! o1 n图表 6 结果下载页面
7 B3 G: q* Y; T+ T. f! K; G1 K3、编程模拟测序结果
2 B2 y, h2 y$ W* U% A; k/ u/ V拷贝数是这里的N值；覆盖度是m，测序深度是宏观的量，在这里与覆盖度意思相同，就是测序仪10X，20X。
单端测序
- [' S5 L9 k% Y0 k( ~! ~/ j9 n" x! f# u
  r" A5 b, d! V0 Z* x7 C; x图表 7 程序模拟单端测序
双端测序
# p3 [9 K: r: g  d0 j
图表 8 程序模拟双端测序
测序结果
1 ?+ q* d! V( Z; n3 S, m+ h0 s
图表 9 结果文件

! g- V# _3 w  V8 Q6 e0 e) {- {# o$ g因为期望片段长度是600bp，在片段长度区间200-1000bp内，所以大部分的片段都没有删除。
' B% T& }& h' C" C测序结果统计表

测序方式        基因组大小(bp)        片段长度区间 (bp)        N值        期望片段长度        克隆保留率        片段数量        Reads长度范围（bp）        Reads总数量        Reads总长度        覆盖度(m值)        理论丢失率(e-m)        覆盖率(1-e-m)
1 w& X# Z5 b! s3 a单端        12157kb        200-1000        10        600        0.95        107378        50-100        101968        7645.541kb        0.62889        0.53318        0.46682
- ?/ L  e* l6 ?9 r: Q单端        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
$ {: }8 W6 p6 d2 ^双端        12157kb        200-1000        20        600        0.95        214212        50-100        407186        30534.265kb        2.51164        0.08114        0.91886
四、讨论和结论

- E( r' v0 x1 ]* v2 q; x程序运行方法
# ]6 h  ]1 L. Q  Y6 c; l1 X1 E2 c
: }, v: f& N, W3 i在类的构造方法init()中，调整参数。
- f; J& m0 m7 e, L2 O- i$ k0 }- QAveragefragmentlength为片段平均的长度；
minfragmentlength和maxfragmentlength是保留片段的范围；
5 {2 {! I% M4 o. w- B1 w+ |cloneRetainprobability是克隆的保留率；
- ]# q" Z3 F/ g1 Y8 \% `minreadslength和maxreadslength是测序reads的长度范围
, O+ [' v$ w" {: V" b0 x8 A
8 G; r  n% L, O9 P9 V3 w模拟测序的诸多方法都封装成了Sequencing类，只需要创建类，并调用singlereadsequencing()和pairreadsequencing()方法，传入文件名的参数即可。

附录
$ O4 t1 P$ O5 s" O  x5 b  n9 z0 B; Y
6 ]+ L6 M' {7 b  S0 qfrom Bio import SeqIO
6 }2 ~2 g( p" X+ ~" H6 z! X' v6 N" sfrom math import exp
import random
* O0 m5 W3 A% y; i/ u1 B: m
class Sequencing:
    # N代表拷贝份数
    def __init__(self)
# G5 n* n# H  M' o6 C1 m# d        self.fragmentList = []
  z( T, R% [, B9 S        self.readsID = 1
        self.readsList = []
        self.averagefragmentlength = 650
        self.minfragmentlength = 500
4 Z+ }! e, j" {; E2 P4 h; m        self.maxfragmentlength = 800
        self.cloneRetainprobability = 1
        self.minreadslength = 50
        self.maxreadslength = 150
$ g7 v, F9 J% C        self.N = 10
        self.genomeLength = 0
        self.allreadslength = 0
4 E1 ^/ D" U8 S- t! [+ u
    # 生成断裂点
    def generatebreakpoint(self, seqlen, averageLength):
        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
0 R/ {: B+ I; j+ D8 L        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
# W' L6 Z6 M: K1 W& T* G$ N        breakpoint.append(seqlen)
        breakpoint.append(0)
7 h6 D/ s" _$ M3 F5 Z        # 把随机断裂点从小到大排序
# n! r( q6 s7 w        breakpoint.sort()
        return breakpoint
2 }" p, ]! z2 J  e( j) ~4 H
, \; g, h2 r8 M6 `6 e    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围: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:
. B) M- Y! ~9 T; U! f" t" Y2 ?                self.fragmentList.append(fragment)
        return self.fragmentList

    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
: c* v7 A) W% n. F( \$ \, a# Z    def clonefragment(self, fragmentList, cloneRetainprobability):
        clonedfragmentList = []
        Lossprobability = [random.random() for _ in range(len(fragmentList))]
/ I" i4 W. Z. @8 K1 R. s        for i in range(len(fragmentList)):
            if Lossprobability <= cloneRetainprobability:
                clonedfragmentList.append(fragmentList)
" p% ^3 B  h" _, D; f        return clonedfragmentList

; r7 i0 ~6 t# E3 z! `  X  w: t3 a    # 模拟单端测序,并修改reads的ID号
    def singleread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
            fragment.id = ""
            fragment.name = ""
* q' N3 v1 I, t+ i; q) K            fragment.description = fragment.description[12:].split(",")[0]
4 |6 u  m3 |9 B6 c6 @) U+ k/ l            fragment.description = str(self.readsID) + "." + fragment.description
1 C, W& q/ O" z; X5 U            self.readsID += 1
            readslength = random.randint(self.minreadslength, self.maxreadslength)
; Y/ N. {) z/ @! _            self.allreadslength += readslength
1 {1 v8 l% U! U0 o8 F# \            self.readsList.append(fragment[:readslength])

1 M2 @1 m2 g) y5 F& y    def singlereadsequencing(self, genomedata, sequencingResult):
        for seq_record in SeqIO.parse(genomedata, "fasta"):
7 \* F1 P3 p9 B" K" f; Z            seqlen = len(seq_record)
0 E9 \, P" \1 z4 T            self.genomeLength += seqlen
; ^3 H$ B% x3 A, E$ r) t            for i in range(self.N):
                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
8 H! Z) s+ O( L# W8 M- N        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟单端测序
        self.singleread(clonedfragmentList)
% z/ B' f% p* ]. P, B        SeqIO.write(self.readsList, sequencingResult, "fasta")
2 c; w; ?6 G8 k# M
/ Q2 _0 I; U6 @- I8 t5 ]% p9 ^    def pairread(self, clonedfragmentList):
* Q. W, @3 N! o) j# Q  D        for fragment in clonedfragmentList:
: S/ V- \6 A! n, q            fragment.id = ""
( s- F4 D/ X& h1 F            fragment.name = ""
            description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + description
            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
1 d+ t7 W3 t: C6 n; w: a            self.readsList.append(fragment[:readslength])
9 Z/ c9 `2 y0 d; b! p. E" T
            readslength = random.randint(self.minreadslength, self.maxreadslength)
" K) a/ _$ u1 p9 r0 I  P) G9 j            self.allreadslength += readslength

            fragmentcomplement = fragment.reverse_complement()
' n1 B/ v6 \, @6 A- p( u            fragmentcomplement.id = ""
            fragmentcomplement.name = ""
" ~5 a+ b& c6 U& n            fragmentcomplement.description = str(self.readsID) + "." + description
/ r# }( b2 f# z3 k; T/ e% D* n            self.readsList.append(fragmentcomplement[:readslength])

# V- Z, c* O- [( M; ?* m9 A1 {            self.readsID += 1

    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
6 Q  _' S/ u4 E5 O        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
. p( \8 n3 _$ o/ t            self.genomeLength += seqlen
6 @! m$ v3 L( ^            for i in range(self.N):
                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟双端测序
" ?/ e9 i' G3 }! |. y        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")
        SeqIO.write(readsList_2, sequencingResult_2, "fasta")
5 x1 u$ b+ ?. R& n: q% M- L
    def resultsummary(self):
        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
# ?: E" W: o2 i! ^) e* c- E* ~. W        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
. K3 g8 }5 Q* B* \3 n; U        print("N值:" + str(self.N))
0 z& S- x" H5 P' k" ^* `: X8 y: u2 j        print("期望片段长度:" + str(self.averagefragmentlength))
* v( T% K/ ]3 H$ l/ R6 Q        print("克隆保留率:" + str(self.cloneRetainprobability))
        print("片段数量:" + str(len(self.fragmentList)))
        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
  J! [# F4 \. H, v        print("reads总数量:" + str(len(self.readsList)))
7 z" d. p- {2 a        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
        m = self.allreadslength / self.genomeLength
. |* j2 l) y4 L) ~8 B! s        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 Y- I# G) z- V  ^sequencingObj = Sequencing()
sequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
: Q; e# e% l% A6 n1 zsequencingObj.resultsummary()
8 k5 g2 V" j1 i/ Q' k5 q
1 h- K, }, W9 @( v4 D# 模拟双端测序
( [; s% D% `9 a4 |+ jsequencingObj = Sequencing()
2 m, @: p, r0 J4 L* g: ~% k' ssequencingObj.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
: P/ Z8 h* c/ m4 Z1 t" P- {import random
1 f* H) b/ v, h* @' Z
9 q' ^* d; ^; V1 m% q, Kclass Sequencing:
- ^. `; W3 ~# [& V# q" j    # N代表拷贝份数
6 H1 R, d# I7 w0 h    def __init__(self):
        self.fragmentList = []
        self.readsID = 1
$ t2 r  Z: |; c- z# |        self.readsList = []
( m6 H" w1 m" D        self.averagefragmentlength = 650
! |& S& `4 h# f        self.minfragmentlength = 500
        self.maxfragmentlength = 800
        self.cloneRetainprobability = 1
. }3 R4 |; P: r* b+ ?        self.minreadslength = 50
) s8 {$ i* a4 ~- _        self.maxreadslength = 150
        self.N = 10
, H, X9 m6 {5 [% H8 W        self.genomeLength = 0
        self.allreadslength = 0
/ x; X" C0 S+ ]3 d3 a8 |) K
    # 生成断裂点
, o$ R! \) J% {# E+ o    def generatebreakpoint(self, seqlen, averageLength):
        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
' ?* S3 F* M& y1 v$ Y- k4 o6 d        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
) o0 {8 J0 ?3 z. P        breakpoint.append(seqlen)
        breakpoint.append(0)
        # 把随机断裂点从小到大排序
        breakpoint.sort()
        return breakpoint
4 d3 M! B- n5 R; s) G& A, _
    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
; @( p6 \, q9 }2 D  p5 W    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
* z3 O' E6 Q3 \- {7 r0 ~4 L        for i in range(len(breakpoint) - 1):
            fragment = seq[breakpoint:breakpoint[i + 1]]
            if maxfragmentlength > len(fragment) > minfragmentlength:
  v4 _% `; w8 K% w+ g5 q# F                self.fragmentList.append(fragment)
. D6 Z& e  Q6 P0 {        return self.fragmentList

    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
6 S: F3 j; O% B! }% y$ J* @    def clonefragment(self, fragmentList, cloneRetainprobability):
4 z8 T% r" P$ D$ K8 s        clonedfragmentList = []
5 p- K7 Y) N# A; J2 U        Lossprobability = [random.random() for _ in range(len(fragmentList))]
4 Z& s) }$ O  X  w( Q; ~! `        for i in range(len(fragmentList)):
, y3 M5 |1 f8 b; w* f' F- x' s: n            if Lossprobability <= cloneRetainprobability:
                clonedfragmentList.append(fragmentList)
9 V/ D* s( f  @0 b: |' e, [( n% A& U        return clonedfragmentList
4 V1 c; E( y- u8 w. I
    # 模拟单端测序,并修改reads的ID号
    def singleread(self, clonedfragmentList):
' o( V1 E# Z( H: d5 {        for fragment in clonedfragmentList:
            fragment.id = ""
0 H1 w" D& z2 n7 p/ R! ~            fragment.name = ""
0 o7 f. C4 O+ i6 `* Q: v" J3 o- B            fragment.description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + fragment.description
            self.readsID += 1
" u! H9 o  C3 L            readslength = random.randint(self.minreadslength, self.maxreadslength)
6 W7 l: z1 l2 v& H            self.allreadslength += readslength
& ?6 q$ |. f) p/ I            self.readsList.append(fragment[:readslength])

    def singlereadsequencing(self, genomedata, sequencingResult):
        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
) {7 G( Y: F6 c* Y. m: q$ U            self.genomeLength += seqlen
" ~8 Q2 _9 v9 W            for i in range(self.N):
2 v6 S1 R" ~+ R) l& h# Q                # 生成断裂点
3 E9 k* w6 b& o                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
9 d4 m4 P0 o6 b. M; \                # 沿断裂点打断基因组
' M% [, y+ l5 ]- Z" g8 B. y- l                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
" @8 u- ^, {) E- N. w( K1 Y. E2 `        # 模拟单端测序
( B2 w) @( L$ T; w% P        self.singleread(clonedfragmentList)
' ^" x* z* G$ v' ?: Y        SeqIO.write(self.readsList, sequencingResult, "fasta")

$ ^. f) @* n9 z    def pairread(self, clonedfragmentList):
        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)
            self.allreadslength += readslength
4 D8 D3 D8 f, ^! c4 }            self.readsList.append(fragment[:readslength])
1 ?, }1 C) n6 U; l1 t1 N2 ^! E
" `0 `/ Z( t' h& S3 a7 n5 K* [            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
4 z5 J4 N, ]3 Y) m
            fragmentcomplement = fragment.reverse_complement()
4 r$ {/ _4 j  w: S! l: o  i            fragmentcomplement.id = ""
5 M+ [5 K  J2 ]' W4 w            fragmentcomplement.name = ""
. y6 A  _0 m0 f0 ]6 ^8 N% T6 Y            fragmentcomplement.description = str(self.readsID) + "." + description
            self.readsList.append(fragmentcomplement[:readslength])

            self.readsID += 1

    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
0 w) W% i- ^6 a1 |& h6 C$ C: ^" V        for seq_record in SeqIO.parse(genomedata, "fasta"):
4 w5 F( b8 E- j  K& r8 O, G! U            seqlen = len(seq_record)
            self.genomeLength += seqlen
7 Q! y5 i& j$ a6 h            for i in range(self.N):
: T) |: d$ R7 \                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
        # 模拟克隆时的随机丢失情况
  {8 ?2 e; a& h        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
, n4 q; o& ^5 k        # 模拟双端测序
9 d4 a9 b0 C3 n, w        self.pairread(clonedfragmentList)
: l+ K) I5 X! k$ C, g5 \2 x        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")
        SeqIO.write(readsList_2, sequencingResult_2, "fasta")

# z0 g9 T# j7 ^( a* N    def resultsummary(self):
        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
' m, E1 s, h7 I7 M9 w        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)))
$ `# N$ q% z8 [0 {0 @3 ^        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
( {# d7 I$ @8 p# l        print("reads总数量:" + str(len(self.readsList)))
        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
; ~' A2 E4 P* J2 b! G        m = self.allreadslength / self.genomeLength
        print("覆盖度(m值):" + str(round(m, 5)))
        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
# -------------------------------------------主程序-------------------------------------------
# 模拟单端测序
# m- r1 y" P" S3 c' x9 O( c5 p5 B) w0 asequencingObj = Sequencing()
sequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
' |" t0 \9 A2 b; QsequencingObj.resultsummary()
' R6 g0 H) V/ z6 L7 t- o( |: _
3 A  B# x5 n; n9 R: n) e# 模拟双端测序
sequencingObj = Sequencing()
) \" J$ g2 R; L# f. R8 |sequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
sequencingObj.resultsummary()

3477959497

不错。。。。。。。。。。。。。。。。。
' h# j3 \# P2 c9 a& C3 [' ]