数学建模社区-数学中国

标题: 基因组测序模拟 [打印本页]
作者: 杨利霞    时间: 2019-4-21 14:56
标题: 基因组测序模拟
基因组测序模拟
% O1 I* F# G8 o1 z- e& B基因组测序模拟1 [1 U$ u: `' w- ?4 p3 C; z
' b" `7 @/ F. ~- Y5 Z
一、摘要
# y9 {  j5 Y6 {2 K# i# ]! g, g8 @) ^! T; D- E7 ?6 e
通过熟悉已有的基因组测序模拟和评估程序，加深全基因组鸟枪法测序原理的理解，并且能够编写程序模拟全基因组鸟枪法测序，理解覆盖度、测序深度、拷贝数等概念，设置测序相关参数，生成单端/双端测序结果文件
8 }8 d4 Y0 ~' U: S% i' c! d9 |1 y8 O2 m* x
二、材料和方法: `4 n+ l+ ?  s1 t5 Z6 [
5 K& \+ i% `: n" ]6 r* Z
1、硬件平台
" `: \) o) }" D$ C8 U  C7 k4 W: n$ v/ T4 S! k' A: R
处理器:Intel（R） Core（TM）i7-4710MQ CPU @ 2.50GHz / j% Q0 P7 `! e( p9 @
安装内存(RAM)：16.0GB
' D& Z4 i  W, x
/ z1 N- S7 p8 I. b8 W2、系统平台- M8 m3 n  D- d$ k* q- R" y
Windows 8.1，Ubuntu% s: K# H5 e9 U( A

3 M7 L  P+ S6 x3、软件平台
: N( Y; g5 [& o9 m2 U; K$ p0 H! G1 V/ Y, [/ s
art_454
9 [' x5 v9 M  \  u9 f' A' t# l: |GenomeABC http://crdd.osdd.net/raghava/genomeabc/
( u3 P( p3 ^9 A; M% _- d# aPython3.5+ ~7 y: c' A7 q7 t1 V# z% g- T
Biopython
8 O( o0 j4 J+ e+ n. A" p4、数据库资源# k% u( l9 p0 b4 D
# }, ?& F! g! u' Y0 C7 k, ~7 Y4 |! A
NCBI数据库：https://www.ncbi.nlm.nih.gov/+ q* B/ \+ x# p; u# y5 U

( ^/ t( u/ s# t2 ^* T5、研究对象# W  P9 a% K$ D  F0 N

' T0 n  ]8 m1 w3 y7 {酵母基因组Saccharomyces cerevisiae S288c (assembly R64) ! e6 k$ h: @* U. _5 G& v) Y+ O
ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_genomic.fna.gz
3 a8 n9 r& ]$ H+ h4 o
8 A+ M% x. W2 C6、方法
+ f" E2 H& i( @3 r( A* a
, F; y1 o9 u9 n$ {art_454的使用 ) B7 E& q% H# {( Z9 l* c% ?5 A3 y6 M
首先至art系列软件的官网，下载软件，在ubuntu系统安装，然后阅读相关参数设置的帮助文档，运行程序。
- }8 q8 q7 Z7 C/ v9 YGenomeABC
# z+ O' X" u; m( z% }进入GenomeABC(http://crdd.osdd.net/raghava/genomeabc/)，输入参数，获得模拟测序结果。
" r9 e; ]9 Y" e# J" c" h' c编程模拟测序
, U, a( {3 z- H. C1 v$ |下载安装python，并且安装biopython扩展模块，编写程序，模拟单端/双端测序。" X1 r8 [# s0 ]% w! p# W$ f6 M
三、结果! z! o9 q* E7 G6 A# [
+ D. V( D5 J+ d9 j( z1 M
1、art_454的运行结果$ j0 [1 {! T  ^8 G1 S4 Q2 F

# L: r4 d$ |  n无参数art_454运行，阅读帮助文档
) Q+ V2 o' B3 b. s
* C, Q& m: [9 F8 D$ m图表 1无参数art_454运行 ( s8 L* h. {/ t
对酵母基因组进行基因组单端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20. : A4 {& i1 a# M3 g' K) V
下图为模拟单端测序，程序运行过程及结果
! o5 `8 d5 `# A+ _5 C( M8 Y' s. q) X+ ]' n0 ^& w' b
图表 2 art454单端测序 ; ?5 C4 t) g+ w. t/ l/ s
0 ^7 g4 b7 p+ I; B) O
图表 3 art454单端模拟结果
* o( ^) F. _" {) \双端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20；MEAN_FRAG_LEN设为1500，即平均片段长度为1500；STD_DEV设为20，即长度的标准差为20 : g5 {5 N' V; p8 b
下图为模拟双端测序，程序运行过程及结果
# L) ^' L! p. h% s
6 G" q/ q. H* l图表 4 art454双端测序 ( t- d8 c# P  @( i1 L

* u% ]1 K- u4 [, C图表 5 art454双端模拟结果
+ K0 I/ @- H! m9 u  S- c0 H2、GenomeABC . Z4 S/ B4 s$ M, L1 Q4 ~$ M
下图为设置参数页面 , `5 t# {7 h+ ^/ N  n

( x( ?4 i, {/ [下图为结果下载页面
/ o8 V0 t$ J+ ]. \+ E1 H+ ]
% u; D- U5 x) l图表 6 结果下载页面
$ C  z4 G& a5 K* D. Y7 A) W3、编程模拟测序结果
3 Y* j6 B5 V5 h& _- [+ r! S. I; L拷贝数是这里的N值；覆盖度是m，测序深度是宏观的量，在这里与覆盖度意思相同，就是测序仪10X，20X。
0 t1 T( X0 J% n! f% W单端测序 , t/ _) z" r6 `4 O6 P# I$ b
" k6 q' X3 w( m4 U0 u/ B& S+ P
图表 7 程序模拟单端测序 ) n( _. l7 d& [# ?6 F8 Y: K
双端测序 0 @) d# a( t* r3 o2 h. L7 c
- ]/ i9 A. T# w2 o- h1 l9 K/ H
图表 8 程序模拟双端测序 1 D/ t$ f( L8 J" @2 J1 i4 G
测序结果
* q5 \( Z( G5 d  D
& I5 h( z) w8 l6 C- o& _图表 9 结果文件  g( Q. z2 v4 {" j+ h9 T

9 V" e9 N8 h7 I% i$ u" e. f; W6 B因为期望片段长度是600bp，在片段长度区间200-1000bp内，所以大部分的片段都没有删除。
0 Q( {; ^7 [: |: }$ M测序结果统计表
& [; G1 }# x7 I5 l2 C; p9 p, Z( M: \0 H
测序方式        基因组大小(bp)        片段长度区间 (bp)        N值        期望片段长度        克隆保留率        片段数量        Reads长度范围（bp）        Reads总数量        Reads总长度        覆盖度(m值)        理论丢失率(e-m)        覆盖率(1-e-m)
  a" \9 _4 p( g' I单端        12157kb        200-1000        10        600        0.95        107378        50-100        101968        7645.541kb        0.62889        0.53318        0.46682
# x; l0 j$ B0 G* C2 E" Q单端        12157kb        200-1000        20        600        0.95        213722        50-100        202996        15227.882kb        1.25259        0.28576        0.71424- J+ H. ?- u7 i; e' ]
双端        12157kb        200-1000        10        600        0.95        106704        50-100        202770        15212.662kb        1.25134        0.28612        0.71388
  L; B! |( d' F) h0 q9 I: I2 N双端        12157kb        200-1000        20        600        0.95        214212        50-100        407186        30534.265kb        2.51164        0.08114        0.918862 O" m% f/ S6 h" O
四、讨论和结论
) Q4 R' G: G0 e8 d7 j6 O) o
7 A- T- F# V* F6 q; t. @+ Y( [4 R程序运行方法: y  V- w3 h  I! j* s4 h
% Q5 [% x: W: d
在类的构造方法init()中，调整参数。 & m8 w8 n1 G* V$ w: G; F
Averagefragmentlength为片段平均的长度； 5 T7 n8 O, H% Q; `
minfragmentlength和maxfragmentlength是保留片段的范围； 1 c0 H& J1 R* e; E1 M4 N& q/ ~; X
cloneRetainprobability是克隆的保留率；
! g, S7 g  u. H2 |0 S- _: f5 Uminreadslength和maxreadslength是测序reads的长度范围
$ y6 j( Q0 N+ @  v! N* R2 o1 \, @+ Q9 [( l% I& j+ U
模拟测序的诸多方法都封装成了Sequencing类，只需要创建类，并调用singlereadsequencing()和pairreadsequencing()方法，传入文件名的参数即可。
; \9 q; D* b; {- D: ~# F% u+ b
& G! e- ~/ K; Q附录
! O  U1 h9 k' j# w6 O5 m5 ~4 w
: }' o1 S5 ?6 ^% Jfrom Bio import SeqIO
5 Y+ ~/ S/ B' `from math import exp# ^0 D5 B: }0 Y* R3 i: D
import random
4 [8 N/ f; r) @8 j( d* U- A
4 _4 i& j$ W) Q( g! M  q# _class Sequencing:
* }6 d& u  i" u5 d    # N代表拷贝份数6 r6 k7 k) J7 h* ~- w; p
    def __init__(self)
6 C- h. g* p& J        self.fragmentList = []7 b6 g: n; U! i
        self.readsID = 10 s! `* M# I1 L2 t
        self.readsList = []
+ Y7 D8 r. V9 w* v& g5 @" L        self.averagefragmentlength = 650
" ~$ ?* C; Q% y1 C( [$ J1 @* m8 J8 i        self.minfragmentlength = 500
7 B. _6 x$ F9 _, t) J" c0 ^        self.maxfragmentlength = 800
' Y: ?* ?% q2 X$ Q        self.cloneRetainprobability = 12 M. C& {/ D' q  P3 j& z! h
        self.minreadslength = 50
" Y0 T: W6 N4 H6 e" E' s        self.maxreadslength = 1506 [+ p& R! Z( |$ f- P/ O
        self.N = 10, Q3 ]6 }/ [0 e( S
        self.genomeLength = 0
1 ~3 F/ [+ J7 k        self.allreadslength = 0! X/ H* w- D+ L# R- E, U2 m

0 Z% X2 W6 G; }    # 生成断裂点, o. [) G5 c( [1 b% O* g, c/ L
    def generatebreakpoint(self, seqlen, averageLength):& t# ?! }9 I( w2 J* f9 o) c
        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
" r& R; q( N" @1 W/ c        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]9 S$ z/ s, ~7 O' w0 z
        breakpoint.append(seqlen)
9 g6 ?: G4 K. W  X1 g( J        breakpoint.append(0)
# H2 o- _. t/ v        # 把随机断裂点从小到大排序
8 Z; Y7 \7 i3 P% y        breakpoint.sort()9 ]$ X' o8 r. `8 _: ~
        return breakpoint
, ?5 {9 f& b% N) S6 b& K% h9 @- g) b5 S3 D' \) {
    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp+ c$ k' s" F4 I
    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
1 V. N$ T- v# N. f5 v4 n        for i in range(len(breakpoint) - 1):
7 t: x* Y  U) F5 D            fragment = seq[breakpoint:breakpoint[i + 1]]
( H+ @; x' w9 F5 {            if maxfragmentlength > len(fragment) > minfragmentlength:
* r( {$ E0 G" D8 J. O+ V                self.fragmentList.append(fragment)
, V' Y$ p, {/ y- b        return self.fragmentList- o) h/ Y( O, r
" Z2 l. k! l1 P. T8 S: g
    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率7 A: Q. U4 W- y' R2 R& s
    def clonefragment(self, fragmentList, cloneRetainprobability):
+ ?2 p$ r- w0 `. ]$ X        clonedfragmentList = []
) J. W- P" v5 Z* b        Lossprobability = [random.random() for _ in range(len(fragmentList))]8 U! [6 |7 }" \' K  f
        for i in range(len(fragmentList)):2 m% }) J' P9 n
            if Lossprobability <= cloneRetainprobability:, f6 s0 W8 n. `) N8 m
                clonedfragmentList.append(fragmentList)
  M( z0 [) @7 A        return clonedfragmentList& e0 [- E( G& S6 I4 O; F" J; H
: v5 ?1 t& d2 @" H" Q2 c6 x% _
    # 模拟单端测序,并修改reads的ID号
% M7 Q, e: Y  J& P2 J% b3 L    def singleread(self, clonedfragmentList):% H; O; A# l- D+ T9 y0 ^
        for fragment in clonedfragmentList:
# K! Q% h5 l5 M5 G$ S            fragment.id = ""
( o/ e! C% c! m$ ^1 O            fragment.name = ""( G9 _/ |; Z% U2 P  A/ ~4 P
            fragment.description = fragment.description[12:].split(",")[0]
; c, v7 l  T- D% u9 W- N$ Q6 N            fragment.description = str(self.readsID) + "." + fragment.description
  S  m8 u* ?/ k: |- ]3 I            self.readsID += 1
1 v8 e0 R% R7 S6 \! {, P& l            readslength = random.randint(self.minreadslength, self.maxreadslength)$ F4 ~# s$ D  f1 p& m, e3 h* c1 u
            self.allreadslength += readslength
1 l5 ^- v1 q% z" z/ t% D! Z' W            self.readsList.append(fragment[:readslength])
; U3 t0 n8 t8 p5 X$ x$ B5 }7 L& b. L: X* e# Y; ]
    def singlereadsequencing(self, genomedata, sequencingResult):
. w. a1 @6 w5 i        for seq_record in SeqIO.parse(genomedata, "fasta"):
6 o; e( I+ |, {; F            seqlen = len(seq_record)
/ a+ A2 _# M3 v# z+ J9 {            self.genomeLength += seqlen
/ B7 O. f: H  C/ [0 N. z8 G$ s            for i in range(self.N):
7 i+ o: b% }- q9 U! ~$ A                # 生成断裂点( Z2 y7 y$ H0 b8 B
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
, ~3 r- J! {) Z: q% b7 R; |) T7 }                # 沿断裂点打断基因组
! Y& ]: V& @6 F* j! S' W                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
8 u9 V0 D7 P& u        # 模拟克隆时的随机丢失情况1 k3 a% `& A# h  G+ b: ]2 b8 b
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability): m2 z4 C  q5 e$ U% X$ L
        # 模拟单端测序; P  W: A0 o# v. ]
        self.singleread(clonedfragmentList)+ L' H$ E8 ?) c# Q7 S9 X- ^0 w
        SeqIO.write(self.readsList, sequencingResult, "fasta")) }5 P3 ?+ l" |: q

! ^$ D; J8 v+ M    def pairread(self, clonedfragmentList):0 I' o$ B7 ^$ E( t/ e
        for fragment in clonedfragmentList:8 K/ n1 \0 \; J  n$ S; ]
            fragment.id = ""2 p9 i) b2 Z# M. T9 z' }
            fragment.name = "". c* Q. |  ~* b1 S
            description = fragment.description[12:].split(",")[0]
! }* S0 M! F0 R* T- ]            fragment.description = str(self.readsID) + "." + description( _& J" j9 A- s
            readslength = random.randint(self.minreadslength, self.maxreadslength)
9 m5 t7 Q/ v- N& a. Q# @& x            self.allreadslength += readslength
+ A, t  P* ?; A$ R            self.readsList.append(fragment[:readslength])
' A  B9 w+ w8 K' l/ H, V* M: e! G/ x) d
            readslength = random.randint(self.minreadslength, self.maxreadslength); J$ f2 [( {0 }5 K* K3 G, t8 B1 m
            self.allreadslength += readslength0 C8 ?* m4 i4 B0 a- l
3 c* s9 M6 r  G7 s1 ~
            fragmentcomplement = fragment.reverse_complement()
$ F+ b8 W; I& l( I* h% R+ a8 N            fragmentcomplement.id = ""4 a0 S) Z, g8 X% @3 d, i8 z
            fragmentcomplement.name = ""/ P) ^- N3 v  ~4 |' `7 A" M% _
            fragmentcomplement.description = str(self.readsID) + "." + description- `; z1 Z. r3 q3 O: F9 y
            self.readsList.append(fragmentcomplement[:readslength])
% ~% P% }6 D0 I2 m
  Y2 Q- h9 M% I/ {3 x, F2 X            self.readsID += 12 i+ Z% C+ E( d; ~9 v

: s  e8 S& z  p$ U; ~    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
, l  w% l# L0 L7 W6 u2 }1 k% k        for seq_record in SeqIO.parse(genomedata, "fasta"):
! U6 }8 P+ o, D$ r3 X7 t0 ^            seqlen = len(seq_record)  n# p! ^+ n2 Q; M
            self.genomeLength += seqlen6 E' u1 |9 H# m' r) K. a. [2 D  S
            for i in range(self.N):' v( d3 \* I% w9 P2 e
                # 生成断裂点2 r6 O# N4 c. i2 K' q
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
$ C3 C9 q" s+ u( k, g+ ~                # 沿断裂点打断基因组* a7 o5 ^: [; W
                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)- E& P: t8 O, ]9 r# E
        # 模拟克隆时的随机丢失情况7 d) u0 H4 u8 j; m9 p6 k
        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
9 R/ C. ]1 P/ J* h" X/ }, X2 ~        # 模拟双端测序
7 @4 S: G5 P# M1 V        self.pairread(clonedfragmentList)1 X$ G: l) v& l+ f5 y
        readsList_1 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 0]* H. B( H4 M. i2 f: k! B4 T' g
        readsList_2 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 1]+ z* |4 w) P% J# W% V5 E
        SeqIO.write(readsList_1, sequencingResult_1, "fasta")$ Z# j5 T: Y' T3 C4 K0 G
        SeqIO.write(readsList_2, sequencingResult_2, "fasta")9 u+ c6 P8 h) b: b! L

  F2 ]" c- B9 {' }    def resultsummary(self):3 ?" o3 Z) J9 @, g" V9 W5 e7 ~, g
        print("基因组长度:" + str(self.genomeLength / 1000) + "kb"): ?, }/ r2 v% v4 Y3 ]$ K( @9 M2 [
        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
3 r8 b4 u" W8 Y/ z        print("N值:" + str(self.N))2 M2 Q( @! a" G/ G; q, T
        print("期望片段长度:" + str(self.averagefragmentlength))
# p  K+ ?+ w/ w9 k/ K        print("克隆保留率:" + str(self.cloneRetainprobability))+ f+ m" g+ w7 X+ b& L
        print("片段数量:" + str(len(self.fragmentList)))6 x% _7 o9 v# w/ N$ t
        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
) m( W$ W; q% J        print("reads总数量:" + str(len(self.readsList)))
) P' h8 d4 X6 C$ n) j        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
3 \& c% v$ k3 V4 l4 q6 u" V' E        m = self.allreadslength / self.genomeLength! T4 H# ?6 Q% _. h
        print("覆盖度(m值):" + str(round(m, 5)))  A3 r- e* t) q/ v0 Y5 f$ O
        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))3 L' r  U+ ~, R2 m7 e' l( t5 C
        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5))), Y& |. \. ^: P: P9 p, L
# -------------------------------------------主程序-------------------------------------------
- A' V: `2 e+ o. J; j! |# 模拟单端测序3 \# L  |1 {+ J
sequencingObj = Sequencing()
6 [# I1 |4 h' z, i% L% jsequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
2 W) _  ~  D2 @7 i$ S0 \sequencingObj.resultsummary()
0 p. U: G0 k2 G& m
0 N/ }6 I0 [/ m9 X6 N2 f# 模拟双端测序
& p' w2 H) e3 T5 u) g. N/ Q! CsequencingObj = Sequencing()
" E- S& `: F# E  q, m4 esequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
  H$ `# T: R* Q3 ?0 A2 c" y/ z% M( [sequencingObj.resultsummary()
3 h( d! f4 l% o3 w+ w! f- ^( O  E9 Gfrom Bio import SeqIO; _& x+ h/ p3 H" q) C- c7 k
from math import exp5 I! \# x9 q6 Y: u! \$ ]+ ]2 a
import random
6 m6 ?9 r/ E: [# t: P7 u0 K* e1 q0 c- P0 ]5 T' Q# Q& h; h3 ?/ p
class Sequencing:* K1 X% l% B0 l* {* |! _
    # N代表拷贝份数( @) ?/ Z! f& _8 \6 m/ _. n2 j' M
    def __init__(self):: v+ q: E- z1 Z$ z; W% d' D
        self.fragmentList = []7 ]. j3 m( n! s. ]9 }
        self.readsID = 13 u$ W8 z1 ^# j. T' U- o
        self.readsList = []
+ Z& W: \. \; [- `' k& @. F        self.averagefragmentlength = 650
& s% {4 o( R0 y8 |1 C$ t$ g# h        self.minfragmentlength = 500/ I9 b& r# P6 x/ e: ^4 i0 L
        self.maxfragmentlength = 800
0 h0 ]7 \5 ^7 V- f" @* i' z        self.cloneRetainprobability = 1
! G* t3 w! @/ U8 F( U: V/ Z! z        self.minreadslength = 50
4 H( Y# G7 s! A% s& p- Y$ [        self.maxreadslength = 1505 D$ ^  V& d+ d1 S* H' _3 E
        self.N = 10
# r0 \% H0 S% q        self.genomeLength = 02 ?5 {9 n  U8 `5 R% [5 w0 [. y
        self.allreadslength = 0) ?$ d/ a4 s( v+ |! X1 L

2 r9 U: K: ]3 ^4 t    # 生成断裂点
7 W; {. A, e2 v( s    def generatebreakpoint(self, seqlen, averageLength):
3 T& ?" d8 `6 A4 R. ]) O        # 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)4 M4 J- t8 Q* c" ]
        breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]1 u3 D$ ?4 H  x7 c3 {: @6 f
        breakpoint.append(seqlen)1 i2 X! ~; g  N
        breakpoint.append(0)8 ]% p' G2 {: B3 r: a* ?
        # 把随机断裂点从小到大排序+ o; G- n/ Q. d2 D3 U
        breakpoint.sort()
* t: `; k6 H$ }4 x, Z; E# r        return breakpoint
# x- \% N2 ?9 A- X: r
! C0 R2 d/ x, {& Z5 `8 `+ l) R    # 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
. n" q- z3 A8 _  T    def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):" @: X5 e: I: v/ h6 E6 P9 q
        for i in range(len(breakpoint) - 1):! l* d6 x' J. _9 v" Z
            fragment = seq[breakpoint:breakpoint[i + 1]]0 y. _- D# ]- G: @. l
            if maxfragmentlength > len(fragment) > minfragmentlength:. g7 R" o* {" B+ C# S
                self.fragmentList.append(fragment)) A8 _2 p& K' f$ E( [7 ]7 T
        return self.fragmentList
- \1 E- l' r0 Y1 Y. n* A& h, R$ g4 N1 y( T, k+ X
    # 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率+ J& r3 w, g$ }. |$ l' U7 K9 F
    def clonefragment(self, fragmentList, cloneRetainprobability):
$ m( g: V! y& {  U        clonedfragmentList = []
+ l% O3 s5 z7 B! _        Lossprobability = [random.random() for _ in range(len(fragmentList))]
+ R6 m$ P1 ^8 E8 h4 B; s        for i in range(len(fragmentList)):: `+ g2 g3 X) ?5 r( ~8 s% @9 }/ b0 [* b
            if Lossprobability <= cloneRetainprobability:
  ]- p0 u& ^) m) {' m% h$ P  i' g& Z                clonedfragmentList.append(fragmentList)5 c& a) l" e' p' \! `$ y: N
        return clonedfragmentList) Y# P, n; g. C! x8 H4 t6 o8 }

8 T, L; N$ f4 K  [( O% J" N  m; ]% Y    # 模拟单端测序,并修改reads的ID号  _, O( Z; N% t3 F! z+ Q
    def singleread(self, clonedfragmentList):
( f, E# Z) T. ^4 y$ T; ]        for fragment in clonedfragmentList:3 B1 _: ?3 x9 a8 A9 [
            fragment.id = ""
: G: c' n" }, w6 H: A2 E            fragment.name = ""
& L9 D) F5 H( I* C% X3 [! G            fragment.description = fragment.description[12:].split(",")[0]8 Q& g9 S% ?! m4 |. f7 D! Q2 T9 {
            fragment.description = str(self.readsID) + "." + fragment.description
) P4 b% m! F& X; }7 J/ p/ B( U1 m8 f) `            self.readsID += 11 S6 f' Z* @  X
            readslength = random.randint(self.minreadslength, self.maxreadslength)
. L: p+ I# }2 W( d3 t            self.allreadslength += readslength) \0 }& j( l- q5 i* ?# `- t- }
            self.readsList.append(fragment[:readslength])
% S* x0 B4 c4 h$ N8 p8 j( B( e9 ]: [1 J9 i5 S, i8 `1 C
    def singlereadsequencing(self, genomedata, sequencingResult):
9 K" W/ K% M/ f5 W6 G        for seq_record in SeqIO.parse(genomedata, "fasta"):
4 k* a# z* y/ b# }            seqlen = len(seq_record)
1 [4 r$ }+ m- H* v' Q            self.genomeLength += seqlen
8 s' o" e. K5 m/ e& b& ^  W            for i in range(self.N):- P$ n+ S5 A! N5 Y0 \
                # 生成断裂点
1 p) q1 k" r; v; U, Z# I                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength); m6 _( O6 X3 Y: V
                # 沿断裂点打断基因组
+ ]2 V1 S# s5 d) h                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)& Q1 s$ z, U1 _1 v
        # 模拟克隆时的随机丢失情况
' W) e  i; x$ R5 h% |7 N        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)3 z7 _2 c( `, B# d$ X' D3 U! A
        # 模拟单端测序2 |$ n; C0 c* \  g  }
        self.singleread(clonedfragmentList)
  a) P( W" m5 ?6 \        SeqIO.write(self.readsList, sequencingResult, "fasta")
6 W& ?% m5 e% J2 J
/ z* Y6 U! D- u  o- N    def pairread(self, clonedfragmentList):3 _+ ~6 {- g, {& w  a) c
        for fragment in clonedfragmentList:
- k; B( G3 A! a            fragment.id = ""
6 T8 T! V1 |( N5 l6 p! `. {            fragment.name = ""
+ M/ v( I9 ?, X4 r            description = fragment.description[12:].split(",")[0]% E, e9 e4 t. C: h5 m  @4 ~
            fragment.description = str(self.readsID) + "." + description+ \& I  ~2 I* n  y
            readslength = random.randint(self.minreadslength, self.maxreadslength)8 E% K0 L1 o! x
            self.allreadslength += readslength8 X; q. k2 @% d7 A
            self.readsList.append(fragment[:readslength])
6 s2 ~( q* q/ J5 S. B) S0 Q6 B# D0 T1 x, E  M" e
            readslength = random.randint(self.minreadslength, self.maxreadslength)
8 ?* ^% |( @+ O+ q& q            self.allreadslength += readslength
. @, d: P& M6 k/ e( ]
# Z" ~& ?0 K# B4 \/ R- Y1 `$ f            fragmentcomplement = fragment.reverse_complement()4 A2 {3 A: }) k2 V/ {; t  E; Y
            fragmentcomplement.id = ""2 m" q$ j' B& v" s% N
            fragmentcomplement.name = ""
! o9 ?7 g* G+ K3 E8 A3 t# |/ ~7 H. u            fragmentcomplement.description = str(self.readsID) + "." + description
6 g" E- q: E: H            self.readsList.append(fragmentcomplement[:readslength]); z9 U. D. ~2 o# G8 ^

# Q0 ^4 r: U' z8 M4 }$ x- X" c: z            self.readsID += 1/ L# C8 _5 v9 N8 s/ Y
5 w0 d! I6 B3 |  Z; j2 }6 {$ x, q
    def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):! ]/ X& ?8 Y; ]8 |
        for seq_record in SeqIO.parse(genomedata, "fasta"):
" y9 z' M' r$ x" P6 j            seqlen = len(seq_record)! b4 T7 u- ^2 D- g
            self.genomeLength += seqlen
7 y4 O8 `& }/ \: `: C            for i in range(self.N):( L( z" N/ n3 l
                # 生成断裂点# m  \  L, w! j, p! O
                breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
# ~. z0 I2 z  h, i. W: w                # 沿断裂点打断基因组
% r( g3 p" r% d, N) M" W                self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)1 ?: |2 Y" }' V5 }! z3 k( x
        # 模拟克隆时的随机丢失情况
$ i0 z- j3 ^4 }8 D4 u4 D9 t        clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)6 q) w- |5 T' O" V% R9 U& D
        # 模拟双端测序) e5 m/ ?1 c( X0 x
        self.pairread(clonedfragmentList)
7 V; l  h* X1 T$ E2 \        readsList_1 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 0]
5 b6 K1 S# D; l& a3 }7 |; b        readsList_2 = [self.readsList for i in range(len(self.readsList)) if i % 2 == 1]
1 K+ M) O$ R; |        SeqIO.write(readsList_1, sequencingResult_1, "fasta")
+ q) e! S* k$ S! i2 w, r7 U$ ]        SeqIO.write(readsList_2, sequencingResult_2, "fasta")
9 G5 h2 W- J5 N
  W( F  F' {9 d/ r    def resultsummary(self):
" G' v8 l$ Y' Q5 J7 U        print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
6 ]# B* l2 a; V6 o$ ^, `        print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
; W; ~% o4 E& {/ d, q; q        print("N值:" + str(self.N))9 T1 i; C% ~3 v
        print("期望片段长度:" + str(self.averagefragmentlength))
2 e4 c1 _* |, v% \        print("克隆保留率:" + str(self.cloneRetainprobability))' S, G/ {0 g" F" T& u, r
        print("片段数量:" + str(len(self.fragmentList)))
8 x$ M1 T. q& U9 @. b        print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
& B  \# P: b6 k9 h9 m0 q2 P        print("reads总数量:" + str(len(self.readsList)))3 ]! @1 {# `7 q( d5 w9 k
        print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
' y" R, ?7 E; i' W2 n. j2 X  y6 L7 h        m = self.allreadslength / self.genomeLength1 {( i$ ~) `' o9 d# q
        print("覆盖度(m值):" + str(round(m, 5)))# U3 D2 U: Y! A: ?" D: S0 x0 D. ?6 q3 d
        print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
+ X* W: p: U/ ?, o  j6 }        print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5))). e3 O8 B3 h, Q+ t) Q2 P5 Z
# -------------------------------------------主程序-------------------------------------------6 ~3 R; |( g- u% a2 g
# 模拟单端测序: V0 H* Y4 l. K, e5 p/ Q6 H" a' w
sequencingObj = Sequencing()
1 z% u% F/ I5 i4 m; M& bsequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa"), l, H* d9 J6 x5 N, G* v; e
sequencingObj.resultsummary(), D7 t! f% m% G

5 n: _  B1 m; M1 K: \# 模拟双端测序4 Q( `4 d& r2 Z( J$ r' ]% \0 L/ q
sequencingObj = Sequencing()& I; b$ M( G* E. L3 n/ V
sequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
! d" \4 h3 J8 T& W3 {sequencingObj.resultsummary()
6 W1 q; w: z6 d* a. q
) K% r/ V& s* C2 _2 s1 {, i
0 z5 a. ^" [2 K! f' g  m5 j
$ k4 [: }' Q2 H: v2 g 8 c+ g' ]( C5 Z

数学建模解题思路与方法.pptx

117.69 KB, 下载次数: 4, 下载积分: 体力 -2 点
作者: 3477959497    时间: 2019-4-22 10:49
不错。。。。。。。。。。。。。。。。。3 t" B5 c( U  @: n




欢迎光临 数学建模社区-数学中国 (http://www.madio.net/) Powered by Discuz! X2.5