基因组测序模拟

杨利霞

基因组测序模拟
0 f9 W5 y7 E' y基因组测序模拟

一、摘要
2 u3 J; v4 @( u% F, b; J- ^9 u2 |
  C7 B+ S7 J8 v通过熟悉已有的基因组测序模拟和评估程序，加深全基因组鸟枪法测序原理的理解，并且能够编写程序模拟全基因组鸟枪法测序，理解覆盖度、测序深度、拷贝数等概念，设置测序相关参数，生成单端/双端测序结果文件
  f9 f/ k# I3 C- L; \1 u
二、材料和方法
4 U9 l5 y9 M" b( [7 X
1、硬件平台
: h8 ^$ _* j; w6 g
) H1 ~4 `, n3 P1 z/ o处理器:Intel（R） Core（TM）i7-4710MQ CPU @ 2.50GHz
- X7 {( z4 S5 |3 Y5 p+ q; H: r" I安装内存(RAM)：16.0GB
  Z& Y) ~, _0 T* b
0 d9 s8 k  T* {: w5 L$ f6 I" ^% Y2、系统平台
Windows 8.1，Ubuntu

3、软件平台
, R/ W$ Y, W& r
art_454
GenomeABC http://crdd.osdd.net/raghava/genomeabc/
Python3.5
Biopython
4、数据库资源

NCBI数据库：https://www.ncbi.nlm.nih.gov/
6 S* I+ ?- ?' \
5、研究对象

酵母基因组Saccharomyces cerevisiae S288c (assembly R64)
, |5 `6 d' F" A% aftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_genomic.fna.gz

6、方法

art_454的使用
首先至art系列软件的官网，下载软件，在ubuntu系统安装，然后阅读相关参数设置的帮助文档，运行程序。
* `$ v" r* r8 C, ]' ^# fGenomeABC
进入GenomeABC(http://crdd.osdd.net/raghava/genomeabc/)，输入参数，获得模拟测序结果。
3 W9 @' S* l# q- t! e/ R. H# j" `编程模拟测序
下载安装python，并且安装biopython扩展模块，编写程序，模拟单端/双端测序。
三、结果

1、art_454的运行结果

无参数art_454运行，阅读帮助文档
+ i1 W/ X9 N5 T5 _/ w
图表 1无参数art_454运行
对酵母基因组进行基因组单端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20.
2 f$ |" v* |4 V* B, W/ O  e' @下图为模拟单端测序，程序运行过程及结果

图表 2 art454单端测序

图表 3 art454单端模拟结果
双端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20；MEAN_FRAG_LEN设为1500，即平均片段长度为1500；STD_DEV设为20，即长度的标准差为20
. O7 W8 T8 s5 W6 m下图为模拟双端测序，程序运行过程及结果

图表 4 art454双端测序
; w. L+ F! S7 E9 I4 ^
8 Y, f; n3 v1 r' m% |/ }图表 5 art454双端模拟结果
2、GenomeABC
; X) K; h4 U# R下图为设置参数页面
" ?. Q" v- b3 _5 D" R  C( c* S$ _2 Q
下图为结果下载页面
1 {9 {! s7 k: N- I* }
- @7 Z/ c# w! g+ H3 o. A# i0 ~图表 6 结果下载页面
7 c4 \/ g( o; a# C: u5 i3、编程模拟测序结果
拷贝数是这里的N值；覆盖度是m，测序深度是宏观的量，在这里与覆盖度意思相同，就是测序仪10X，20X。
单端测序
9 r9 [1 @( g' R: ], z
图表 7 程序模拟单端测序
双端测序
2 R, v% T' P  J/ d5 d6 K3 U+ {
, s6 a1 [/ M3 b  P) S0 S2 H. T图表 8 程序模拟双端测序
测序结果

图表 9 结果文件
1 v; F5 \" h" z/ D* ]- h; H
2 Q' p# W+ \- u5 S  z9 K5 C/ M因为期望片段长度是600bp，在片段长度区间200-1000bp内，所以大部分的片段都没有删除。
测序结果统计表

6 ^* i4 F; C# H7 V. Y2 e测序方式        基因组大小(bp)        片段长度区间 (bp)        N值        期望片段长度        克隆保留率        片段数量        Reads长度范围（bp）        Reads总数量        Reads总长度        覆盖度(m值)        理论丢失率(e-m)        覆盖率(1-e-m)
4 f' e" j# g# ]" P, {' [单端        12157kb        200-1000        10        600        0.95        107378        50-100        101968        7645.541kb        0.62889        0.53318        0.46682
( w6 p) G6 t+ a% B单端        12157kb        200-1000        20        600        0.95        213722        50-100        202996        15227.882kb        1.25259        0.28576        0.71424
) x3 _1 B1 c7 X! O+ ~双端        12157kb        200-1000        10        600        0.95        106704        50-100        202770        15212.662kb        1.25134        0.28612        0.71388
3 }. S9 Y: m1 i$ T: J6 \- |双端        12157kb        200-1000        20        600        0.95        214212        50-100        407186        30534.265kb        2.51164        0.08114        0.91886
% i6 @( C1 K) ]5 ?4 S8 }$ d四、讨论和结论

/ ?4 t+ {- z, U) y5 v6 e$ D程序运行方法
; b- g1 L; `) ^3 g3 U2 U
在类的构造方法init()中，调整参数。
Averagefragmentlength为片段平均的长度；
minfragmentlength和maxfragmentlength是保留片段的范围；
cloneRetainprobability是克隆的保留率；
; l% p8 ~3 ^" l# d( }minreadslength和maxreadslength是测序reads的长度范围
7 _- }' J9 M6 n. Z: ], Z. `- [
模拟测序的诸多方法都封装成了Sequencing类，只需要创建类，并调用singlereadsequencing()和pairreadsequencing()方法，传入文件名的参数即可。

. U  Y5 V: F! k附录

  b5 Q6 s  c, {" s7 I, O% ^; Gfrom Bio import SeqIO
$ i& x- c9 F6 @' w2 l- j8 Ffrom math import exp
import random
" b; U$ s( t: N: e% ]3 H
/ M/ {6 z- a1 Q  O7 pclass Sequencing:
    # N代表拷贝份数
/ `  z5 X" r- l4 `( H    def __init__(self)
6 q- Z5 D3 I# [! Z2 w8 ?        self.fragmentList = []
. g' d4 D" U9 P, U2 |: y, `$ T        self.readsID = 1
- |2 A" o& f1 j* g, C- {        self.readsList = []
        self.averagefragmentlength = 650
3 r% G( y" E6 `+ R        self.minfragmentlength = 500
8 y: M; K' W) {2 E6 j- \        self.maxfragmentlength = 800
) ?; F6 I( g  l. d3 [% l$ P        self.cloneRetainprobability = 1
        self.minreadslength = 50
        self.maxreadslength = 150
        self.N = 10
        self.genomeLength = 0
5 K% _2 ?: {" L! `        self.allreadslength = 0

    # 生成断裂点
    def generatebreakpoint(self, seqlen, averageLength):
# j0 |4 U( k- c$ N9 P9 V9 T7 z        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
3 A0 h" M8 I5 I" c  L, w: e        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
! H7 O# n$ G5 h# \3 v- L& T1 l# W        breakpoint.append(seqlen)
, ^: ^, Q6 P+ q! P+ n1 |        breakpoint.append(0)
" h1 d6 T- c1 y" {( o+ ]        # 把随机断裂点从小到大排序
        breakpoint.sort()
: f; y( F% o3 Z6 l! O! x: k1 L        return breakpoint

    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
  @2 ~5 P6 z* E! J  [% w& X    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
        for i in range(len(breakpoint) - 1):
8 Q8 C% F* Q$ Q+ {4 U2 ]/ m5 O( R            fragment = seq[breakpoint:breakpoint[i + 1]]
* [+ m7 L* _. g  C$ G$ ^( d            if maxfragmentlength > len(fragment) > minfragmentlength:
* ~8 {, s& S$ c: g" e3 i/ \                self.fragmentList.append(fragment)
1 B  c$ @+ }0 B0 m8 ~        return self.fragmentList

1 w; d+ y. }/ [: f* L    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
    def clonefragment(self, fragmentList, cloneRetainprobability):
# ^; G: @9 V1 ], v& F" [        clonedfragmentList = []
/ g9 F( v- y1 }* @) U8 K. M9 M        Lossprobability = [random.random() for _ in range(len(fragmentList))]
        for i in range(len(fragmentList)):
& @. M5 H' w/ V! o4 @            if Lossprobability <= cloneRetainprobability:
                clonedfragmentList.append(fragmentList)
# L! R5 ~) ~2 H7 A2 {- I        return clonedfragmentList

0 W* e, t5 m2 c! D    # 模拟单端测序,并修改reads的ID号
    def singleread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
9 i, p$ L% e: t1 I) J. {            fragment.id = ""
            fragment.name = ""
            fragment.description = fragment.description[12:].split(",")[0]
: d9 C4 E1 [7 R# _2 I            fragment.description = str(self.readsID) + "." + fragment.description
            self.readsID += 1
            readslength = random.randint(self.minreadslength, self.maxreadslength)
( d) U. K! ^8 Z! g/ y            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])

, r# q9 t! S" q- F: T    def singlereadsequencing(self, genomedata, sequencingResult):
        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
: S6 N2 S) |! s! `            self.genomeLength += seqlen
            for i in range(self.N):
                # 生成断裂点
& X  [* w, x! l4 Y) k3 p                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟单端测序
        self.singleread(clonedfragmentList)
) B( o) K0 k6 r6 `1 R0 q        SeqIO.write(self.readsList, sequencingResult, "fasta")

    def pairread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
            fragment.id = ""
            fragment.name = ""
' K3 k- X$ h8 m5 v: b; `            description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + description
$ q  x; W+ U% m, g$ f! b, [) J            readslength = random.randint(self.minreadslength, self.maxreadslength)
' a# `* F- J4 Q, S$ K7 z            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])

            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
  ?) T1 S& l  n7 F
2 f7 D5 L6 X' [0 g4 d            fragmentcomplement = fragment.reverse_complement()
            fragmentcomplement.id = ""
            fragmentcomplement.name = ""
* E: E" j8 U7 |) N/ T1 g' w            fragmentcomplement.description = str(self.readsID) + "." + description
" K& N( D* l: f5 z9 J0 W            self.readsList.append(fragmentcomplement[:readslength])

            self.readsID += 1

8 F$ c  A& V/ U+ I    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
- G- L1 D% C" H0 ]2 x$ \5 E        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
            self.genomeLength += seqlen
, S9 ^# N6 c6 i# e. [& L0 ~            for i in range(self.N):
* p3 ?: z8 l, ?% V% b: c0 F7 m  F                # 生成断裂点
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
0 q. ^$ K$ _; t! A2 ]% S                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟双端测序
: o+ F, H3 j9 Z& U* ]5 S        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]
5 f" f( @5 v6 j$ @" ~        SeqIO.write(readsList_1, sequencingResult_1, "fasta")
        SeqIO.write(readsList_2, sequencingResult_2, "fasta")

, v: x1 H' g: v. S+ k6 {    def resultsummary(self):
$ A2 q  x- W+ N" D% j        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
7 e& E0 H  r) C* K        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
9 N9 @: J  b$ g  Z        print("N值:" + str(self.N))
/ T* `% S& t3 g5 p& |, ~        print("期望片段长度:" + str(self.averagefragmentlength))
        print("克隆保留率:" + str(self.cloneRetainprobability))
0 H% |. K; s. X" `5 L9 o        print("片段数量:" + str(len(self.fragmentList)))
        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
7 l2 p1 U9 G" f$ J. D        print("reads总数量:" + str(len(self.readsList)))
  @/ w, |+ f! {4 Z' p        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
% u; W- `( N) [, c+ J4 G- I        m = self.allreadslength / self.genomeLength
        print("覆盖度(m值):" + str(round(m, 5)))
        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
# E5 W5 d9 }, M" S/ f* ]7 j. h1 @        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
# -------------------------------------------主程序-------------------------------------------
* _4 z7 E. n2 ~3 e/ t- Z, K# 模拟单端测序
5 \  @: c; ^) }9 \  [sequencingObj = Sequencing()
sequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
sequencingObj.resultsummary()
9 `3 D- _* g1 J$ B! Z4 N
# 模拟双端测序
sequencingObj = Sequencing()
. C7 T( ~( N3 `" e; z  r; vsequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
7 g1 _! L& [. Y! D9 QsequencingObj.resultsummary()
from Bio import SeqIO
0 Y* |( {+ I, T+ E% ^0 `3 M$ ^. Kfrom math import exp
import random

# S3 d; M2 [' @3 wclass Sequencing:
    # N代表拷贝份数
5 q# u7 z4 k" q* X7 {1 S5 J+ [: M! a    def __init__(self):
& `; u/ y+ J6 @( l+ \8 r        self.fragmentList = []
        self.readsID = 1
        self.readsList = []
        self.averagefragmentlength = 650
' _/ e" w1 n+ l7 O7 n# J6 E2 g8 |: V. ~        self.minfragmentlength = 500
* x( t7 Q5 r: [& B. X. N5 i6 d0 A        self.maxfragmentlength = 800
        self.cloneRetainprobability = 1
        self.minreadslength = 50
" Y5 X# _( P# t% V        self.maxreadslength = 150
; `8 n1 m3 b: z& `; f        self.N = 10
        self.genomeLength = 0
        self.allreadslength = 0
. c7 w/ y+ _+ z8 s7 l/ z4 x
    # 生成断裂点
    def generatebreakpoint(self, seqlen, averageLength):
        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
1 |4 `1 E; I3 r        breakpoint.append(seqlen)
        breakpoint.append(0)
- j/ Z) V3 T* ~0 K- A- {7 O        # 把随机断裂点从小到大排序
        breakpoint.sort()
        return breakpoint
2 `( ?8 ^: A6 n& u0 q3 }: j  k0 m
  c/ x4 g' ^1 B' d" u    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
1 q. Y9 G! R+ @! \- l! z6 [    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
        for i in range(len(breakpoint) - 1):
. K4 v; E2 E# ]9 i, Y# l  l            fragment = seq[breakpoint:breakpoint[i + 1]]
1 q5 {7 a: J( `3 U# r            if maxfragmentlength > len(fragment) > minfragmentlength:
) F; P7 v/ h3 J) ]! ^) M8 q6 ?  f                self.fragmentList.append(fragment)
        return self.fragmentList

* x0 E2 E5 O  U; E6 r    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
' x  d5 t5 {5 H4 U    def clonefragment(self, fragmentList, cloneRetainprobability):
. G* g6 [' N/ z5 m0 E9 h6 K; ]        clonedfragmentList = []
/ V8 v; m" L+ O3 U        Lossprobability = [random.random() for _ in range(len(fragmentList))]
        for i in range(len(fragmentList)):
            if Lossprobability <= cloneRetainprobability:
                clonedfragmentList.append(fragmentList)
        return clonedfragmentList

    # 模拟单端测序,并修改reads的ID号
  r) |  G0 t* ~0 ]# s& B    def singleread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
; h: `2 o; S, ]& q2 O4 c            fragment.id = ""
            fragment.name = ""
9 l$ j+ @# q$ K6 U            fragment.description = fragment.description[12:].split(",")[0]
            fragment.description = str(self.readsID) + "." + fragment.description
            self.readsID += 1
- |+ a5 g6 U8 V% \+ d            readslength = random.randint(self.minreadslength, self.maxreadslength)
: B8 C) T" d# l3 \6 _( O            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])

$ t& V4 i2 \* R1 k0 A# R6 s" V    def singlereadsequencing(self, genomedata, sequencingResult):
/ `8 e8 |% X) i* O4 V' R  ]2 N9 \7 z        for seq_record in SeqIO.parse(genomedata, "fasta"):
            seqlen = len(seq_record)
* u9 q- U; U  J. _% r( p: n. l            self.genomeLength += seqlen
            for i in range(self.N):
                # 生成断裂点
& ~+ o! ~- _: ^, l" ]$ a/ `. A                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
9 J# K+ d. \: p1 ^+ i# ~7 k# Z                # 沿断裂点打断基因组
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
7 @8 g# y0 ~. _; W# z, u  X        # 模拟克隆时的随机丢失情况
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
" e: Y$ Z4 }; B# V" S9 H7 M        # 模拟单端测序
        self.singleread(clonedfragmentList)
        SeqIO.write(self.readsList, sequencingResult, "fasta")

    def pairread(self, clonedfragmentList):
        for fragment in clonedfragmentList:
            fragment.id = ""
            fragment.name = ""
/ a5 e. T/ m" b2 ~. _$ a            description = fragment.description[12:].split(",")[0]
7 h# |% e6 i) U            fragment.description = str(self.readsID) + "." + description
            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
            self.readsList.append(fragment[:readslength])
6 [7 g* C; A4 _+ o1 M  s( `! R
( p! m% N; M3 n& R- p% S            readslength = random.randint(self.minreadslength, self.maxreadslength)
            self.allreadslength += readslength
1 I6 Z- Y( j3 e, y9 l6 D
            fragmentcomplement = fragment.reverse_complement()
  [9 T) K+ ]; V/ \) B            fragmentcomplement.id = ""
2 M5 `4 O0 F: {9 W; X0 t7 |            fragmentcomplement.name = ""
            fragmentcomplement.description = str(self.readsID) + "." + description
" V$ Z& s7 ?) L; F* o            self.readsList.append(fragmentcomplement[:readslength])
4 H7 N& z% I& A
% i% o5 O1 ?! [  K  R            self.readsID += 1
& `" p0 Q: B7 u4 J
5 Q- Y$ Z# w% a! ]4 a    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
6 U5 q! W3 e) S8 S( Q5 ^        for seq_record in SeqIO.parse(genomedata, "fasta"):
1 m4 @  Y+ T' J# f) u) K- j            seqlen = len(seq_record)
            self.genomeLength += seqlen
7 Z6 U- p6 w3 [            for i in range(self.N):
                # 生成断裂点
6 ?) D) F3 t, K$ W/ d' Y. G                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
. d/ {- z; E# J! x: l6 o                # 沿断裂点打断基因组
7 N2 k" K# I; w, @+ Y  E                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
2 C3 T* s# i+ V' p$ {$ T$ F  [6 n        # 模拟克隆时的随机丢失情况
* _, X. l1 {& r4 s5 u7 I        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
        # 模拟双端测序
6 d1 C$ v. J" ]6 b/ l. W: q        self.pairread(clonedfragmentList)
        readsList_1 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 0]
  c8 f- I9 V0 s- E. ]2 F        readsList_2 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 1]
$ H/ [; \$ x4 t0 k; D        SeqIO.write(readsList_1, sequencingResult_1, "fasta")
9 R) f! W, |9 k$ ?2 I4 J        SeqIO.write(readsList_2, sequencingResult_2, "fasta")
1 m$ p* `. X& x  ?% H' S/ R' s
& V" ^5 H. H/ Z! u  J0 w; o) J5 p    def resultsummary(self):
        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
* {3 E$ D2 q0 c        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
, o+ ^3 f7 L% m6 b# N5 M' V1 o, n* R        print("N值:" + str(self.N))
        print("期望片段长度:" + str(self.averagefragmentlength))
; r9 [* a' X3 A1 d6 L        print("克隆保留率:" + str(self.cloneRetainprobability))
+ i* X- X3 M, j5 {* r+ V        print("片段数量:" + str(len(self.fragmentList)))
6 j1 m/ I% V& U7 E" |' N        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
( G4 m- b# z/ A( i3 R6 `        print("reads总数量:" + str(len(self.readsList)))
        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
6 D! A  I# g8 G& w. s2 a5 [        m = self.allreadslength / self.genomeLength
        print("覆盖度(m值):" + str(round(m, 5)))
3 z7 G9 b7 h: u% v/ n, g5 K! w5 V        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
2 t% t1 w" N# t; p9 f5 @        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
# -------------------------------------------主程序-------------------------------------------
1 C* ~7 M: t) w% S# 模拟单端测序
* J" }' l$ `$ w1 SsequencingObj = Sequencing()
sequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
sequencingObj.resultsummary()

/ y1 D/ @7 j" J' Y# 模拟双端测序
$ K9 F3 L8 r3 y' ~3 {sequencingObj = Sequencing()
/ }) I  Q" C, x: W  ?/ ]# a% ysequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
sequencingObj.resultsummary()

* z. Y! ~1 b0 ^: L. J
' E0 s9 a. R1 ~, t, W

7 z2 N' W- H  P1 E% Y  Q

3477959497

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