基因组测序模拟

杨利霞

基因组测序模拟
9 b2 {  B7 v' J) s4 s基因组测序模拟
% ~  z4 A9 B0 M" W; K+ D3 A( _# p
. l, x: ]0 }5 J+ E一、摘要

3 F$ b( P! q1 v. U9 R4 C( b# N通过熟悉已有的基因组测序模拟和评估程序，加深全基因组鸟枪法测序原理的理解，并且能够编写程序模拟全基因组鸟枪法测序，理解覆盖度、测序深度、拷贝数等概念，设置测序相关参数，生成单端/双端测序结果文件
0 [0 Y8 U" ~2 Y2 R# {
  ]# N+ q* M. A) |4 P二、材料和方法
$ x; P: O1 A$ V& E& ?$ t
, T% |/ n2 c. p# B. a: {/ B1、硬件平台

- n/ E5 `. ~2 F8 j处理器:Intel（R） Core（TM）i7-4710MQ CPU @ 2.50GHz
. Q1 H6 w  B% H, C! ~安装内存(RAM)：16.0GB

2、系统平台
! Q' ?. ?# A# z, OWindows 8.1，Ubuntu
5 h" w* Z( q/ Y6 Z. S3 U
3、软件平台

art_454
2 M# q, {- g6 e8 FGenomeABC http://crdd.osdd.net/raghava/genomeabc/
: ?2 ]9 g; J% |( Y+ h7 t9 L2 KPython3.5
$ K' f; d  t! i4 {Biopython
) z' q" D1 n& E4 `4、数据库资源

- I$ g  [8 S4 y  `& F1 PNCBI数据库：https://www.ncbi.nlm.nih.gov/
5 d; S+ `1 @* P( J
$ ]: f) l+ m+ [7 h  u8 x5、研究对象

酵母基因组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、方法
) H6 p) [. R9 m1 f1 ~- s
  o! t* A7 h3 ~: Cart_454的使用
首先至art系列软件的官网，下载软件，在ubuntu系统安装，然后阅读相关参数设置的帮助文档，运行程序。
GenomeABC
进入GenomeABC(http://crdd.osdd.net/raghava/genomeabc/)，输入参数，获得模拟测序结果。
编程模拟测序
下载安装python，并且安装biopython扩展模块，编写程序，模拟单端/双端测序。
. m( o; n( U( A三、结果
) Q3 X8 r& X( a
1、art_454的运行结果

) }- P  |  x# P* j7 d; \无参数art_454运行，阅读帮助文档
" ^  D8 ]& t6 H
图表 1无参数art_454运行
0 O( C& l" y9 a对酵母基因组进行基因组单端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20.
4 a* V1 z* G+ w- h- W下图为模拟单端测序，程序运行过程及结果

# i  k4 H8 v/ P6 W# J( Z: s图表 2 art454单端测序
" Y$ ?2 v& w& u% ?" j  b2 O- A
& f/ \- b) E* }8 W! Y/ s* |5 V图表 3 art454单端模拟结果
双端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20；MEAN_FRAG_LEN设为1500，即平均片段长度为1500；STD_DEV设为20，即长度的标准差为20
. P& v& }' O# N+ w8 o0 Q2 [( `下图为模拟双端测序，程序运行过程及结果
  m: j, C7 C. [* [0 Q/ Q
图表 4 art454双端测序

$ y( S; [7 w; Q4 t& u图表 5 art454双端模拟结果
2、GenomeABC
6 ^( A) v0 ]" L# E6 y下图为设置参数页面
( p  F4 i# k/ U# C! M! r
下图为结果下载页面

. d6 t# q' }+ u! Q图表 6 结果下载页面
* b1 S2 H) u$ s) ^/ y3、编程模拟测序结果
( w! ?3 c* H9 V/ q/ e$ r' F拷贝数是这里的N值；覆盖度是m，测序深度是宏观的量，在这里与覆盖度意思相同，就是测序仪10X，20X。
& t# L' f$ r9 R8 v单端测序
3 x  P, I7 {0 Y! J. R. h# G( U
图表 7 程序模拟单端测序
6 `6 u# l7 f% x0 G双端测序
* F- q. l5 l- W- Z* S
: P& ^" L. q1 R$ k图表 8 程序模拟双端测序
  R3 e% T) w9 X5 g& u) H# D' d' L测序结果
2 F. V1 W$ l; l- U0 {
; o0 `+ D( u; u/ Q  I# Z! `$ u; p图表 9 结果文件

因为期望片段长度是600bp，在片段长度区间200-1000bp内，所以大部分的片段都没有删除。
% {2 X" d* X% q( L测序结果统计表

测序方式        基因组大小(bp)        片段长度区间 (bp)        N值        期望片段长度        克隆保留率        片段数量        Reads长度范围（bp）        Reads总数量        Reads总长度        覆盖度(m值)        理论丢失率(e-m)        覆盖率(1-e-m)
! p- I; ^: x) B9 r1 Q) l' v0 e单端        12157kb        200-1000        10        600        0.95        107378        50-100        101968        7645.541kb        0.62889        0.53318        0.46682
" ~7 c- T% c( `, f- y9 @单端        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
四、讨论和结论

程序运行方法
; _, x( o; a9 g
在类的构造方法init()中，调整参数。
Averagefragmentlength为片段平均的长度；
0 }8 D9 F. i* k* q2 _3 gminfragmentlength和maxfragmentlength是保留片段的范围；
cloneRetainprobability是克隆的保留率；
& c9 q! H+ B, f- p6 E6 aminreadslength和maxreadslength是测序reads的长度范围

模拟测序的诸多方法都封装成了Sequencing类，只需要创建类，并调用singlereadsequencing()和pairreadsequencing()方法，传入文件名的参数即可。

附录

from Bio import SeqIO
) x$ l7 b- K% c% bfrom math import exp
import random
6 {/ j8 x9 N9 x: |+ p( {
class Sequencing:
+ T5 C! m6 m9 a    # N代表拷贝份数
    def __init__(self)
; K% ~& a1 p% G% F" I; k        self.fragmentList = []
; s, N9 X" o9 \& f. H. B8 @        self.readsID = 1
( H+ N0 B+ f1 n+ l% r/ ]" U        self.readsList = []
        self.averagefragmentlength = 650
0 m1 |! l  u/ w& ~" B& q6 Z% ^        self.minfragmentlength = 500
1 m+ I) Z4 M/ `* k        self.maxfragmentlength = 800
( i9 x) W5 j: ?9 L. e0 h        self.cloneRetainprobability = 1
        self.minreadslength = 50
6 @, v; {% e5 ~0 T- p" U# ~        self.maxreadslength = 150
6 R7 P9 l8 z$ A6 B6 z* J6 b  R! V        self.N = 10
        self.genomeLength = 0
        self.allreadslength = 0
, Q; M! A0 c0 L/ e5 W
    # 生成断裂点
    def generatebreakpoint(self, seqlen, averageLength):
- C; d) ^6 x4 R. I1 j$ k, P        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
8 A4 X& k2 U- a4 M        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
        breakpoint.append(seqlen)
        breakpoint.append(0)
* ]" ]* D4 ^; g5 B% o        # 把随机断裂点从小到大排序
: J5 u" @* S9 g1 ^& i* s) d# r        breakpoint.sort()
1 t& v' Q. B, [7 |; V. f+ K# z- B        return breakpoint

    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
2 q# O4 k; O. E        for i in range(len(breakpoint) - 1):
            fragment = seq[breakpoint:breakpoint[i + 1]]
            if maxfragmentlength > len(fragment) > minfragmentlength:
                self.fragmentList.append(fragment)
4 y" [- Q( u6 e        return self.fragmentList

/ v2 |6 h$ T" J    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
* z, t5 i& [9 p    def clonefragment(self, fragmentList, cloneRetainprobability):
        clonedfragmentList = []
2 C# k: g; v4 w4 ^# p        Lossprobability = [random.random() for _ in range(len(fragmentList))]
: F) d' j6 k4 r4 h0 H- b        for i in range(len(fragmentList)):
8 `- h' Y! X) y* Y0 ~8 T% D7 u% w7 v            if Lossprobability <= cloneRetainprobability:
+ D4 U7 q" \$ [# }4 P/ O7 A                clonedfragmentList.append(fragmentList)
        return clonedfragmentList
5 I# r. W- Q8 F6 u: w
    # 模拟单端测序,并修改reads的ID号
" n# `  ^. t' F: q2 o' ?    def singleread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
1 F/ @% u0 F# O3 T0 x            fragment.id = ""
, L2 b+ I2 h* L9 Z            fragment.name = ""
            fragment.description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + fragment.description
            self.readsID += 1
* r, E9 Q6 \! Y  L1 o            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
. L$ o. g) l; x. ?5 R# k            self.readsList.append(fragment[:readslength])
% k, l( _& K1 W! C1 p/ P  m$ a7 _% u
    def singlereadsequencing(self, genomedata, sequencingResult):
7 R" X, v4 ^; v3 u$ r6 d9 m, W! r1 z        for seq_record in SeqIO.parse(genomedata, "fasta"):
- r- r" q( A. _* a( y            seqlen = len(seq_record)
# W0 `3 S1 l+ [3 K# T  Z; O            self.genomeLength += seqlen
  q- H# \0 i: T            for i in range(self.N):
; p3 ]! p# ]/ a! m+ n# s2 C                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
" M4 L, q4 q4 R* K, V; I* }                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
        # 模拟克隆时的随机丢失情况
+ \$ q' M: `/ d$ Y( X! |9 Q" V        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
3 b' `( p4 ]& o( u) w% q        # 模拟单端测序
        self.singleread(clonedfragmentList)
- S9 C! c  T2 _8 Y9 |( B/ p/ A        SeqIO.write(self.readsList, sequencingResult, "fasta")

. ^0 {4 l, v% M* X    def pairread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
% h% Y( o: f( g9 ~3 q            fragment.id = ""
            fragment.name = ""
4 B4 M# g$ [$ H' H9 ?2 a            description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + description
            readslength = random.randint(self.minreadslength, self.maxreadslength)
8 f& B, E" V3 f" L* E            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])
1 ?% y4 {8 q( H
            readslength = random.randint(self.minreadslength, self.maxreadslength)
: [: N: s! m) z. B" k, ^: r            self.allreadslength += readslength

            fragmentcomplement = fragment.reverse_complement()
            fragmentcomplement.id = ""
            fragmentcomplement.name = ""
            fragmentcomplement.description = str(self.readsID) + "." + description
8 C! c8 R' [) ]            self.readsList.append(fragmentcomplement[:readslength])
2 Q: {% X2 i3 `( ^0 B/ p
            self.readsID += 1

# X0 O% k7 G: E$ @, }    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
. g1 [7 H- d! _9 ]  m; I4 G& i9 l0 F7 T        for seq_record in SeqIO.parse(genomedata, "fasta"):
$ S- S5 @& r- P( v8 [            seqlen = len(seq_record)
7 p" D0 d" d9 \5 n- c7 b) B            self.genomeLength += seqlen
            for i in range(self.N):
; G9 t3 H5 o8 K' O9 s1 P6 N' p: [' L, ?                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
; p! f* f8 N" \  I  @                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
4 T& f# f5 J! H; M7 b0 ]. c        # 模拟双端测序
2 H' ?% K! v' j3 E0 }1 C5 A; ]! ]        self.pairread(clonedfragmentList)
5 o$ O# V" n% t        readsList_1 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 0]
) n! q/ h2 H. c9 _( c/ ?, o1 y6 s& I" t        readsList_2 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 1]
        SeqIO.write(readsList_1, sequencingResult_1, "fasta")
0 J3 [7 y) H  b. K1 x0 v2 r1 M        SeqIO.write(readsList_2, sequencingResult_2, "fasta")

    def resultsummary(self):
; _; ~& p7 `3 }7 {3 X5 F: ]        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
# L& O" K+ X& {; k% j        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
        print("N值:" + str(self.N))
        print("期望片段长度:" + str(self.averagefragmentlength))
4 H* V! P+ W6 X  k: t  a        print("克隆保留率:" + str(self.cloneRetainprobability))
        print("片段数量:" + str(len(self.fragmentList)))
        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
/ K  N/ A. l7 Y2 z        print("reads总数量:" + str(len(self.readsList)))
        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
        m = self.allreadslength / self.genomeLength
        print("覆盖度(m值):" + str(round(m, 5)))
        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
: w( Y" a$ O* |& B5 f        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
# -------------------------------------------主程序-------------------------------------------
# 模拟单端测序
sequencingObj = Sequencing()
sequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
- m' u1 C' i1 \9 ^/ \1 LsequencingObj.resultsummary()
- o/ p$ {, [8 Z# q# J  S
' G4 a+ X* M3 [8 k+ [  d- [0 P: K9 h# 模拟双端测序
/ {2 F" y. c# S* Z) SsequencingObj = Sequencing()
3 l5 F# w! U% L1 u* ^3 g$ S# [sequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
sequencingObj.resultsummary()
from Bio import SeqIO
6 f$ h" K2 N- I1 r) Hfrom math import exp
0 }% _" \) I* ^import random

0 b" E8 ]' |7 }. mclass Sequencing:
    # N代表拷贝份数
) i' @0 c; ?9 v- y    def __init__(self):
        self.fragmentList = []
        self.readsID = 1
        self.readsList = []
        self.averagefragmentlength = 650
        self.minfragmentlength = 500
        self.maxfragmentlength = 800
        self.cloneRetainprobability = 1
/ D  w" W. R8 b+ P8 z        self.minreadslength = 50
        self.maxreadslength = 150
8 l6 Q3 r+ \& [        self.N = 10
6 G* U; e& M& _% ^        self.genomeLength = 0
        self.allreadslength = 0
# i2 }9 d' p+ `
. I9 t/ {% u+ D) B# l  e  `    # 生成断裂点
    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)
        # 把随机断裂点从小到大排序
        breakpoint.sort()
        return breakpoint

    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
1 g% J8 R/ p$ Z; \/ q    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
( H) D: j# q/ R. \; Q9 N* Y        for i in range(len(breakpoint) - 1):
            fragment = seq[breakpoint:breakpoint[i + 1]]
            if maxfragmentlength > len(fragment) > minfragmentlength:
2 Z+ m  x3 l' ^7 [- [                self.fragmentList.append(fragment)
        return self.fragmentList
7 z8 M+ n5 D# U
    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
' |- H+ \  ^( G+ C# V    def clonefragment(self, fragmentList, cloneRetainprobability):
        clonedfragmentList = []
$ w# T: M& v) V        Lossprobability = [random.random() for _ in range(len(fragmentList))]
. B( q" b. s) W, F8 a( l+ H/ K1 b9 N5 ^        for i in range(len(fragmentList)):
            if Lossprobability <= cloneRetainprobability:
                clonedfragmentList.append(fragmentList)
        return clonedfragmentList

    # 模拟单端测序,并修改reads的ID号
# t) W: u! P* G/ z7 o4 v    def singleread(self, clonedfragmentList):
8 A9 L/ u" B* Y; x# s5 t        for fragment in clonedfragmentList:
& Q3 x: B5 V1 \7 d1 S7 o5 @$ }9 }            fragment.id = ""
" e0 {" d6 _5 k) V            fragment.name = ""
            fragment.description = fragment.description[12:].split(",")[0]
1 {4 W# B1 }7 T% `0 S8 Z            fragment.description = str(self.readsID) + "." + fragment.description
            self.readsID += 1
2 q8 L) ?5 G- ~0 k% ?3 r* w            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
4 q; ]5 k* }& e) p/ Z8 {9 ]. X: v            self.readsList.append(fragment[:readslength])

( M/ R% M, w) r  n" @0 G    def singlereadsequencing(self, genomedata, sequencingResult):
& B3 B5 X' O5 ~6 x, A2 W. @8 A% w        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
            self.genomeLength += seqlen
            for i in range(self.N):
                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
5 T1 m/ o8 n% E+ I5 y; d                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
        # 模拟克隆时的随机丢失情况
7 Q2 g: e5 L; K4 V1 e        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟单端测序
        self.singleread(clonedfragmentList)
        SeqIO.write(self.readsList, sequencingResult, "fasta")

    def pairread(self, clonedfragmentList):
. h. x- S% V/ w7 N1 w4 M1 O. [* @1 ~        for fragment in clonedfragmentList:
            fragment.id = ""
% [/ U+ V, }2 O7 ^% p            fragment.name = ""
            description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + description
# K8 G" I+ ?# H" R  A            readslength = random.randint(self.minreadslength, self.maxreadslength)
. p% H- r1 X- _6 Z% E            self.allreadslength += readslength
0 A, A: X* Z# }8 S            self.readsList.append(fragment[:readslength])
; `7 o8 X9 T7 M/ D
            readslength = random.randint(self.minreadslength, self.maxreadslength)
* f, t9 ]1 ]% {3 ~            self.allreadslength += readslength

            fragmentcomplement = fragment.reverse_complement()
8 `& U' u  A1 v+ z* U! m            fragmentcomplement.id = ""
            fragmentcomplement.name = ""
            fragmentcomplement.description = str(self.readsID) + "." + description
) q' U1 H6 c" f6 Y            self.readsList.append(fragmentcomplement[:readslength])
; k5 A5 f8 w0 F/ o1 {
            self.readsID += 1
. ?7 T+ j( c/ Z; X3 `
    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
7 A3 z( [3 G7 W        for seq_record in SeqIO.parse(genomedata, "fasta"):
# ~! q$ \9 t# V9 s. z            seqlen = len(seq_record)
            self.genomeLength += seqlen
5 ^3 u& C; k5 D  B: |( I3 R2 X( i            for i in range(self.N):
! w. C/ Y( V! V3 M4 L                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
  ^- q( n2 u8 Z        # 模拟克隆时的随机丢失情况
+ Z/ r% e+ x! U! I        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
3 B- A- T3 M# Y" O8 O3 k, u& _        # 模拟双端测序
0 |# K- \$ O: t( m9 f        self.pairread(clonedfragmentList)
% l* @( D, B  ?- Y        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]
% J1 m2 k- Q1 d        SeqIO.write(readsList_1, sequencingResult_1, "fasta")
        SeqIO.write(readsList_2, sequencingResult_2, "fasta")
) \, ~- R6 |* @* p
4 ]  n: t: \! Q# u6 i5 {: g- J' h: n    def resultsummary(self):
        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
% L4 t5 N% v3 h& R" _$ B        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
        print("N值:" + str(self.N))
        print("期望片段长度:" + str(self.averagefragmentlength))
8 W7 W* M5 K9 G3 C* {4 d& ~        print("克隆保留率:" + str(self.cloneRetainprobability))
1 {: ]  }" z* v1 ]8 B        print("片段数量:" + str(len(self.fragmentList)))
* l! W( v! h5 O7 h+ D" |        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
1 X( ^: V- Z2 s: y5 a! Q/ c        print("reads总数量:" + str(len(self.readsList)))
        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
        m = self.allreadslength / self.genomeLength
6 ~3 Z8 ~0 n- D% a- M        print("覆盖度(m值):" + str(round(m, 5)))
        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
6 z! K1 L! P& H: [2 N+ k        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
( I' D0 u8 t' c. \! ^% c9 V+ g9 W4 w# -------------------------------------------主程序-------------------------------------------
& A/ T7 @, S5 Z# 模拟单端测序
0 M3 h0 x- _! U. MsequencingObj = Sequencing()
sequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
sequencingObj.resultsummary()
  Y9 A; H6 C3 J, _6 g3 s! Y& L2 r6 S
# 模拟双端测序
: Q4 ~* t, A) c0 b/ d) qsequencingObj = Sequencing()
sequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
# O! h6 a7 d1 G1 OsequencingObj.resultsummary()



; ^  b) A' i8 _; l+ ^6 B

3477959497

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