数学建模社区-数学中国
标题:
基因组测序模拟
[打印本页]
作者:
杨利霞
时间:
2019-4-21 14:56
标题:
基因组测序模拟
基因组测序模拟
+ [/ [1 I' u# f: B+ _! ]
基因组测序模拟
8 {, [8 ~ w5 w" P
' r: y/ a# K2 J/ C
一、摘要
" S) h9 s7 a/ \* v: i3 V
. n3 c: |8 m8 L9 u! l) c
通过熟悉已有的基因组测序模拟和评估程序，加深全基因组鸟枪法测序原理的理解，并且能够编写程序模拟全基因组鸟枪法测序，理解覆盖度、测序深度、拷贝数等概念，设置测序相关参数，生成单端/双端测序结果文件
% \! V+ M- C) E; N/ I
6 I# L9 d; i$ x
二、材料和方法
$ u+ a$ L/ l$ p) ?# m
: `% V, F& v: W
1、硬件平台
! {: U& f9 o1 u4 j4 L3 X4 H, q
$ y* G: N# f$ G
处理器:Intel（R） Core（TM）i7-4710MQ CPU
@
2.50GHz
0 n8 ]' Y9 ?& U: y* J p
安装内存(RAM)：16.0GB
8 X0 \% Z: R! ?0 m0 l7 h: H
3 L% Q6 w3 j5 u5 v! V
2、系统平台
; |- R5 [- }2 B1 p4 G- _* _6 a0 E
Windows 8.1，Ubuntu
$ ^: d T/ m$ Y8 h7 m: W0 U
% E: l$ V8 O8 }5 V$ B
3、软件平台
. {; ^$ N- s) l. d* I; ~" ], r
; Z5 E T% G; U& B
art_454
$ o5 f' C$ q6 G' e- Y1 g& o
GenomeABC http://crdd.osdd.net/raghava/genomeabc/
( ^/ V7 `8 X' ?9 q) d8 o: [
Python3.5
5 s# J: j& C& g
Biopython
( F/ N# C4 ^. S9 I0 N( Z& @
4、数据库资源
/ A2 `4 Q# d6 s% m! |
- P# w3 u# [5 u' h( V) D
NCBI数据库：https://www.ncbi.nlm.nih.gov/
' w, O3 z; u* v& K5 {/ T) B4 i5 o
7 F! Y# |1 G, @9 m
5、研究对象
" E0 n5 d. u& ]5 V
1 q& R5 h& n0 K$ m
酵母基因组Saccharomyces cerevisiae S288c (assembly R64)
& A- S/ }* F, I8 F* h! d
ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/146/045/GCF_000146045.2_R64/GCF_000146045.2_R64_genomic.fna.gz
1 K7 ?) a' }& g, t
* v$ X6 M- n0 L3 W* ~( ?5 o; J
6、方法
+ e+ O% M6 p, J% s6 I# V$ W
Z( c/ h5 k7 ]. k6 k
art_454的使用
$ \, h$ g V* I N
首先至art系列软件的官网，下载软件，在ubuntu系统安装，然后阅读相关参数设置的帮助文档，运行程序。
& }- s; ?2 J' ]( A7 {- l
GenomeABC
% \' c7 ?- d8 _4 O
进入GenomeABC(http://crdd.osdd.net/raghava/genomeabc/)，输入参数，获得模拟测序结果。
2 y3 j1 s! X. I b1 }5 l- N
编程模拟测序
" b4 {, D# V$ E$ |6 V5 h6 i
下载安装python，并且安装biopython扩展模块，编写程序，模拟单端/双端测序。
2 f4 }; _3 ?6 R3 A3 O- g. X8 V
三、结果
: F! D3 F& R* s/ ]' Y7 I+ [
3 q! w/ ?! u( Y. F$ L9 d H5 a
1、art_454的运行结果
9 w: i1 Q+ k2 U2 V
% ?6 F" q4 i7 U# H8 S
无参数art_454运行，阅读帮助文档
' [0 I. g( w! Z0 k0 s& _, @
7 r( h+ s1 T. _; ?( s( r& N
图表 1无参数art_454运行
2 a/ j: L h2 O6 ?: y
对酵母基因组进行基因组单端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20.
" n2 w% V2 O: L( _- ]% ?0 }* r
下图为模拟单端测序，程序运行过程及结果
9 Y4 t) E+ } @- n _
3 ^$ M3 P7 o/ H* J
图表 2 art454单端测序
7 w# K0 p6 E7 F, L. C; f
: I2 U$ G# E/ u
图表 3 art454单端模拟结果
8 e% ^& v5 t1 j% A( J8 }' N
双端测序模拟，FOLD_COVERAGE设为20，即覆盖度为20；MEAN_FRAG_LEN设为1500，即平均片段长度为1500；STD_DEV设为20，即长度的标准差为20
0 @; v/ C0 @5 w
下图为模拟双端测序，程序运行过程及结果
8 X) e) I# [) j% @8 D9 G
. N. E8 Z- J7 ]% Z' M7 Z6 a1 p }; W
图表 4 art454双端测序
% k/ P7 A3 l& r5 y
9 I6 L2 Z3 {( u7 }8 w5 {% s2 F+ @
图表 5 art454双端模拟结果
- W/ C" f0 i s" B0 Z
2、GenomeABC
% {5 Z3 x: |. Y$ H2 W4 K5 T
下图为设置参数页面
& S( f2 b/ y1 _9 ~' B8 |# x
3 S$ a) O6 `) q5 \: z2 N
下图为结果下载页面
9 c0 r- i& v* Y7 F% i9 H
0 z, i \7 ? w6 U
图表 6 结果下载页面
* _/ f9 o. T) ?( E0 s8 J g
3、编程模拟测序结果
. [ a/ e5 J2 f+ P2 _4 C
拷贝数是这里的N值；覆盖度是m，测序深度是宏观的量，在这里与覆盖度意思相同，就是测序仪10X，20X。
3 i1 L* m0 w: R
单端测序
7 x, c' ^$ \* B3 n6 h
8 c* r/ E" ]7 h- I
图表 7 程序模拟单端测序
: i( f1 H1 {- v6 w
双端测序
; V( M! U- H7 `
. l. i7 F9 [# l2 v, \, ?
图表 8 程序模拟双端测序
{( `: Y6 `, `" |; ]4 M
测序结果
. t9 ~$ `& P: W4 ?5 i
5 ~; b3 [# K8 ~1 G4 l w
图表 9 结果文件
& ^7 s3 W* ?3 h3 k9 s0 f0 P$ d
% Z( z! _8 H7 ]* j6 F" l5 j
因为期望片段长度是600bp，在片段长度区间200-1000bp内，所以大部分的片段都没有删除。
+ ?# H9 J2 Y1 Y4 A! d$ k8 u7 I
测序结果统计表
- k% s1 k6 z5 b! _8 v
3 C" `( o7 m# {4 i
测序方式 基因组大小(bp) 片段长度区间 (bp) N值 期望片段长度 克隆保留率 片段数量 Reads长度范围（bp） Reads总数量 Reads总长度 覆盖度(m值) 理论丢失率(e-m) 覆盖率(1-e-m)
, M1 C9 t' G8 `3 n% }8 n6 q! {, E/ p
单端 12157kb 200-1000 10 600 0.95 107378 50-100 101968 7645.541kb 0.62889 0.53318 0.46682
+ L4 S5 C+ x5 h( l9 |; I
单端 12157kb 200-1000 20 600 0.95 213722 50-100 202996 15227.882kb 1.25259 0.28576 0.71424
/ |3 w) V0 m d: e0 x5 Q/ G8 g
双端 12157kb 200-1000 10 600 0.95 106704 50-100 202770 15212.662kb 1.25134 0.28612 0.71388
5 A3 \0 ]0 ?. F4 O" x9 v7 k0 E# ]
双端 12157kb 200-1000 20 600 0.95 214212 50-100 407186 30534.265kb 2.51164 0.08114 0.91886
" s/ Q0 \( `2 G' C
四、讨论和结论
3 a1 A! s% e# `2 k+ J
, }0 I' e3 ?: T! L
程序运行方法
0 S5 t1 c8 k3 |; V5 A. P8 E
; e# E7 \: V8 u$ `" V& x6 H
在类的构造方法init()中，调整参数。
( X+ S9 S' Q! s" @9 d& M. I
Averagefragmentlength为片段平均的长度；
4 }; T+ @ V( ~! h$ i
minfragmentlength和maxfragmentlength是保留片段的范围；
. s1 [- n( L/ R3 r" _" f l+ I
cloneRetainprobability是克隆的保留率；
4 [6 c9 f5 _7 D
minreadslength和maxreadslength是测序reads的长度范围
& ?6 J" ~! b4 w6 W* r; e/ ?; B
* g2 E5 E& O1 l* o' s
模拟测序的诸多方法都封装成了Sequencing类，只需要创建类，并调用singlereadsequencing()和pairreadsequencing()方法，传入文件名的参数即可。
( }" z+ c! W) c2 N, _
- Y/ c4 b# q5 P% a9 J* W
附录
- ^2 w. k& d' E& Y
9 o( w& d# D- r9 w/ {9 @7 P% c: f
from Bio import SeqIO
- M- F+ a# T6 d3 E1 P: Q1 P
from math import exp
5 p3 x& ^! c! J. A$ X% L! E
import random
$ B/ q6 X9 d0 N9 t
3 T: X& v8 \# ^1 [+ F! M0 _
class Sequencing:
4 r! \& c8 I( [1 c! X. w+ s9 o
# N代表拷贝份数
; V+ B' U/ c3 g' o' F
def __init__(self)
) H7 A( b+ {- L2 d) @
self.fragmentList = []
- W' Y) @5 J) Z4 e
self.readsID = 1
& k" ^8 u1 C, f$ Y4 k
self.readsList = []
2 V3 ]2 d% a8 g% M) @
self.averagefragmentlength = 650
* C) c5 L* R' p7 v
self.minfragmentlength = 500
1 }( ?" o7 f& s2 ^/ |( Z
self.maxfragmentlength = 800
6 s/ h& O( d+ I! O8 R
self.cloneRetainprobability = 1
! Z/ f N* C& g
self.minreadslength = 50
/ ]+ o% E$ D2 z! T4 X
self.maxreadslength = 150
% E: e# n: `" u6 s' N3 y% A
self.N = 10
3 n/ s5 g& e7 D! U0 x
self.genomeLength = 0
& s5 ]5 ^$ h7 p
self.allreadslength = 0
- x+ Q/ c5 a [4 n
8 G- z* B/ E7 e8 \5 @) B
# 生成断裂点
" t. ~5 j) I, A* q
def generatebreakpoint(self, seqlen, averageLength):
+ s( l) l5 E! D3 U! a, H7 k
# 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
1 q# L+ D9 P1 k. `
breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
0 B% ^- ?6 X; O+ \/ q4 l% m1 X) ]) ?
breakpoint.append(seqlen)
2 Z5 C, L8 D; J! L
breakpoint.append(0)
; r* S# b2 p2 k: g0 U+ Z7 T* I
# 把随机断裂点从小到大排序
% Q2 U: g1 ~* T* q% m
breakpoint.sort()
& K, s: O' \7 ?+ X6 D
return breakpoint
# m$ }+ l. i9 k' D6 I$ m' e
! L; i% M7 ^ ?9 K$ J( ~% k d! Z
# 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
+ F5 c4 ~* `0 y0 K9 Q$ b
def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
& s% x! O' ~! j1 [6 w
for i in range(len(breakpoint) - 1):
, S; p6 `) W5 Y2 v$ s
fragment = seq[breakpoint
:breakpoint[i + 1]]
2 T; Y( D0 ?$ |6 f
if maxfragmentlength > len(fragment) > minfragmentlength:
4 @, T A" Y9 v0 g' C& u
self.fragmentList.append(fragment)
* k& K! K3 d F( u- L% E! i; y
return self.fragmentList
* V) L& n. S1 {
8 c; B& N* N4 W
# 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
4 f$ t' X. A6 L# u; o9 z+ o' d
def clonefragment(self, fragmentList, cloneRetainprobability):
( t# P8 }1 e4 f; U- `) t
clonedfragmentList = []
" F+ C' @' G G8 h9 x; l
Lossprobability = [random.random() for _ in range(len(fragmentList))]
^7 p/ K6 I; ]. ]1 M
for i in range(len(fragmentList)):
. S9 m4 x4 \+ [% l
if Lossprobability
<= cloneRetainprobability:
6 h% T2 T! Y0 {3 Q, I
clonedfragmentList.append(fragmentList
)
2 Y. W9 b1 n2 r6 O+ W
return clonedfragmentList
: y3 a: m- Y# v+ V. \* Z1 \
' \1 [1 B3 k% X. ]& r! N B
# 模拟单端测序,并修改reads的ID号
: o& W2 l& M' r8 [9 J
def singleread(self, clonedfragmentList):
/ l4 @) X3 l2 a) D- j
for fragment in clonedfragmentList:
$ Y' k( ~/ n4 M3 ?% n" [
fragment.id = ""
$ |5 V. m3 w! j; A/ A
fragment.name = ""
' p2 _- _! ~+ \" ?8 E
fragment.description = fragment.description[12:].split(",")[0]
6 Z1 ^; U4 z( a2 ?1 O/ a
fragment.description = str(self.readsID) + "." + fragment.description
& l, P& v% s6 |: R( L
self.readsID += 1
! @5 b4 m! Q) c4 M: M: B
readslength = random.randint(self.minreadslength, self.maxreadslength)
0 L# V/ V6 y E
self.allreadslength += readslength
* G" {) b; a$ ]5 {- }8 ?, l
self.readsList.append(fragment[:readslength])
2 U) K+ M& E* {% x
' d$ M/ {7 y# L
def singlereadsequencing(self, genomedata, sequencingResult):
5 S! b: F3 k# Q" j/ s
for seq_record in SeqIO.parse(genomedata, "fasta"):
@+ l# w b; R+ i
seqlen = len(seq_record)
" ~0 l6 T& h' P* [
self.genomeLength += seqlen
4 M3 C" N/ R2 q# K# a
for i in range(self.N):
2 H" v. k5 T6 F$ @/ { j
# 生成断裂点
. ]6 I& | W- h0 M! S1 y/ O
breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
$ J5 ~; S* A: }) C6 R& J
# 沿断裂点打断基因组
3 r3 A5 p2 g- L' H L
self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
9 x9 W/ K) ^9 s$ v; q
# 模拟克隆时的随机丢失情况
) r4 D2 k: C6 e7 {/ L. c5 R$ O2 |% j
clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
A5 [3 K' x Z. G6 X
# 模拟单端测序
. h. C/ H: ?/ \2 d
self.singleread(clonedfragmentList)
" g% D" a8 A, z P( f3 P
SeqIO.write(self.readsList, sequencingResult, "fasta")
& r! t/ H+ a0 a8 L ~( ^( A* ^
. I! ]8 ^% Q# o9 g! y# ?5 _
def pairread(self, clonedfragmentList):
, P) d$ o7 d: Z+ X
for fragment in clonedfragmentList:
- `3 t' W2 y/ U. ?4 W7 ^4 Y' b2 N
fragment.id = ""
( J' B) U, d) p/ i$ i: D0 M
fragment.name = ""
/ e! K6 p" u L9 b" `& T
description = fragment.description[12:].split(",")[0]
* N; c. {) i r- G7 b( G/ Z
fragment.description = str(self.readsID) + "." + description
9 ?9 I7 q5 h3 l( L7 F0 R
readslength = random.randint(self.minreadslength, self.maxreadslength)
2 w' w( e' w+ J
self.allreadslength += readslength
5 a9 p* d9 g' B% ]' a f8 V/ ]
self.readsList.append(fragment[:readslength])
) g/ {* u' I) ^- m6 o( ^. s
" H. u/ b/ a' O! l) h
readslength = random.randint(self.minreadslength, self.maxreadslength)
4 n( z5 d# x( Y# P" x; J/ q
self.allreadslength += readslength
6 e& X8 Q# {1 p+ b7 ~. t+ \( I
+ y8 N& K. }/ v" S
fragmentcomplement = fragment.reverse_complement()
9 K- ?' [& x) u$ o8 {; x2 f
fragmentcomplement.id = ""
( ~4 v4 |! d( M7 s3 j# Y: [
fragmentcomplement.name = ""
+ K3 {6 m0 H, X' q+ n
fragmentcomplement.description = str(self.readsID) + "." + description
1 H( g. h5 F" l l0 {, G8 K: O3 T1 B
self.readsList.append(fragmentcomplement[:readslength])
3 T4 M8 G; Y; }, g5 ^
6 X5 s% N' D- R) f4 J
self.readsID += 1
5 E- w) W. p; [2 x" E$ r' l
6 B6 }- Y; Q* i' J- F
def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
- ~$ E6 f* l, _4 ^$ M7 a' d
for seq_record in SeqIO.parse(genomedata, "fasta"):
0 w6 ~7 v8 e" c9 a
seqlen = len(seq_record)
- P. \) B$ [$ F- {
self.genomeLength += seqlen
6 k$ o4 K3 ]- h1 |9 ?" T
for i in range(self.N):
& \" m( F) j8 n7 W6 y
# 生成断裂点
. F1 X/ u r- l3 A
breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
% ~" C# D9 g# e: l# [! q
# 沿断裂点打断基因组
: s$ x, U3 W. ?4 ^
self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
) B" @6 H* ~, P8 n0 {7 K
# 模拟克隆时的随机丢失情况
$ _& N3 ~# [0 m# Q( Q
clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
4 q0 ?* @5 r) f! G
# 模拟双端测序
& `* `5 S3 d* F/ w5 B9 C
self.pairread(clonedfragmentList)
9 Z+ Z( h3 k3 n' _
readsList_1 = [self.readsList
for i in range(len(self.readsList)) if i % 2 == 0]
& _1 X2 [7 g8 K' s1 _( X |
readsList_2 = [self.readsList
for i in range(len(self.readsList)) if i % 2 == 1]
# [& a; z: I8 a2 N/ J3 }- {1 g
SeqIO.write(readsList_1, sequencingResult_1, "fasta")
0 M' [+ }/ V. \, P$ e
SeqIO.write(readsList_2, sequencingResult_2, "fasta")
) d4 q" H3 ~6 [# M5 ~& r+ D
7 J0 b- V9 C7 ?# a
def resultsummary(self):
" X: U! m# _5 o1 c. n; A" q8 y' ?! k
print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
0 ?1 }. b* M3 x, d1 N
print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
( [# G0 W" z9 z8 m! }% d
print("N值:" + str(self.N))
9 l5 _ y) h% J' U& ?3 {
print("期望片段长度:" + str(self.averagefragmentlength))
+ e k8 m T5 ]3 P- P2 c$ ?9 d+ k9 A3 d
print("克隆保留率:" + str(self.cloneRetainprobability))
$ f( {3 O; f+ W! |# [5 ^
print("片段数量:" + str(len(self.fragmentList)))
1 k9 p2 I/ j0 d4 _
print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
" _" u' V o6 C4 s9 k: K
print("reads总数量:" + str(len(self.readsList)))
, `& v, }9 x! [
print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
4 o" _* @! w1 ^* ^
m = self.allreadslength / self.genomeLength
! o) t( N, e- i8 O/ t1 N4 {
print("覆盖度(m值):" + str(round(m, 5)))
- C) I4 {; D1 D! G2 \: i5 h- o
print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
s( P/ R* |8 C- [$ o
print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
- Q# t8 E7 m( X. @& N) w' T1 o* C
# -------------------------------------------主程序-------------------------------------------
, l& n: |$ ^4 w6 d8 w1 S
# 模拟单端测序
4 E5 O# W' o$ Z9 }% i' u
sequencingObj = Sequencing()
A/ X! F1 [2 F# m$ l: F
sequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
- z8 \9 B0 w" E* K0 ]9 m2 ~% ]
sequencingObj.resultsummary()
+ H2 U- u! H. P* {$ x& Q4 s, e6 p
$ G+ w( ]9 |# r. v8 @
# 模拟双端测序
6 N* Z( P4 a' n* X/ Y0 x
sequencingObj = Sequencing()
. q z3 o& U) ~
sequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
+ t. i: l. m7 U. r0 Z
sequencingObj.resultsummary()
, [# W0 F- g' ]' v. [
from Bio import SeqIO
f6 l4 w0 W. ^8 Y9 h$ ?
from math import exp
4 u4 w. f" t/ J( u1 \; Y4 J# ^
import random
8 z# I6 y W3 q3 b& X! K5 B, P
3 `% W) t& H) N! a+ n4 @( ?
class Sequencing:
8 l5 r# O. d/ i. O# h) ]
# N代表拷贝份数
2 ?3 A V; H: L4 C' V
def __init__(self):
5 I9 w* d. p5 Q# E# \
self.fragmentList = []
' X7 X% Y1 p* k* R! [# r) |) f
self.readsID = 1
: N8 W' [ @8 C) P, y, `) I
self.readsList = []
0 a, d5 C+ ^2 p6 m) F+ D
self.averagefragmentlength = 650
2 w* Y* `! s( o3 d' n8 O! R$ p- K* J6 p
self.minfragmentlength = 500
4 s$ i! Y1 ]0 H+ A
self.maxfragmentlength = 800
# M3 l7 s4 V' S
self.cloneRetainprobability = 1
6 U k. t/ u: C& T0 h6 U
self.minreadslength = 50
( e$ X* z/ @* k* w( z3 r% H! L, M
self.maxreadslength = 150
0 q, V$ o% v7 C/ j9 p2 I; |2 A2 v0 G% f
self.N = 10
/ Q# F+ M! H+ c" V: j' Z% R% K4 E
self.genomeLength = 0
; P3 V: }2 G: H0 x
self.allreadslength = 0
, N4 [4 T0 l4 j3 O: J3 I& j) N
! a6 n- s4 Y6 N" H
# 生成断裂点
# z! _2 c3 f) ?$ }# `
def generatebreakpoint(self, seqlen, averageLength):
* u5 z" Y- X. o
# 假设平均每500bp 产生一个断裂点(averageLength = 500),通过随机函数生成seqlen/500个随机断裂点(1到seqlen之间的随机整数)
' ]- t, C. c4 \
breakpoint = [random.randint(0, seqlen) for _ in range(int(seqlen / averageLength))]
3 {7 e4 Q! B7 a j I) Z
breakpoint.append(seqlen)
5 n( }1 ^( v6 E3 M# e
breakpoint.append(0)
% [ u& }( T+ o) s# s2 X4 E- r
# 把随机断裂点从小到大排序
9 j5 p: G6 b. N
breakpoint.sort()
' V8 d. p& j7 N- X2 L4 C0 K
return breakpoint
, x B# p8 o- e+ n D( q
* |6 x) Q: [7 y5 o/ ^8 h
# 沿断裂点打断基因组,并删除不符合长度要求的序列片段,定义片段范围:200-1000bp
3 U4 C0 Y# K3 M4 f2 @, I4 ~; d5 L
def breakgenome(self, seq, breakpoint, minfragmentlength, maxfragmentlength):
+ i( b& A2 V! @
for i in range(len(breakpoint) - 1):
$ z0 S$ U" @# O! C# q Q. B
fragment = seq[breakpoint
:breakpoint[i + 1]]
7 M8 x. ]; n3 D$ c% `) j
if maxfragmentlength > len(fragment) > minfragmentlength:
3 y! W! c" P; i* X9 a. H9 E$ |0 }
self.fragmentList.append(fragment)
* C4 D8 K0 |5 |
return self.fragmentList
_' z# `* m" C6 Z5 \1 B6 W0 [' D) a
) ^3 }* k" z9 O' m3 F( p6 r# Q
# 模拟克隆时的随机丢失情况,random.random()生成0-1的保留概率
, s: K' H9 \0 _0 d6 C( ~8 H! S; k
def clonefragment(self, fragmentList, cloneRetainprobability):
V4 o$ D0 S: j/ @8 r1 N
clonedfragmentList = []
1 w9 ?! J4 ?* v1 T8 C
Lossprobability = [random.random() for _ in range(len(fragmentList))]
' t. f$ B. M. M" b+ E
for i in range(len(fragmentList)):
" y7 Q' e3 |7 O" F8 J
if Lossprobability
<= cloneRetainprobability:
; w2 d) X8 m( F) X) V
clonedfragmentList.append(fragmentList
)
5 K- w, | T6 Y3 L( s
return clonedfragmentList
{: ^% W/ J6 S- {, V0 ?
' o9 k5 y3 j3 @3 k% S4 j: O& j& B1 ?
# 模拟单端测序,并修改reads的ID号
8 [; N) [) o) S4 U$ j5 S
def singleread(self, clonedfragmentList):
) S( D- G1 c& y$ K! z" ^3 e X8 W
for fragment in clonedfragmentList:
3 ^8 s8 u$ O6 a: N1 p
fragment.id = ""
9 e p! U5 o T# s: t( X
fragment.name = ""
; q3 x# @4 j) k- |6 Y
fragment.description = fragment.description[12:].split(",")[0]
$ `- ^/ _- }4 r- E% T6 {. m/ V9 o( a
fragment.description = str(self.readsID) + "." + fragment.description
5 f3 b1 W! B$ h0 E6 m. v
self.readsID += 1
9 }9 _7 m6 X* ~/ K( A- }' F6 o# m+ @
readslength = random.randint(self.minreadslength, self.maxreadslength)
) ^0 q4 t2 y% f. T/ g
self.allreadslength += readslength
) y, X$ ~' c& S2 b" @8 v& E6 f% ]- f0 x
self.readsList.append(fragment[:readslength])
) [. X7 Z0 K# u" u0 }3 X: h
" k0 j) g8 |! s" m+ \2 |2 K
def singlereadsequencing(self, genomedata, sequencingResult):
7 w; l m9 ^) c
for seq_record in SeqIO.parse(genomedata, "fasta"):
) w9 f! s/ K4 |6 t( q6 w
seqlen = len(seq_record)
' B1 X0 A% L0 V4 H/ P7 m+ ]1 K e
self.genomeLength += seqlen
7 k: j e0 Z7 M" T
for i in range(self.N):
; r: H- x: G7 q( A9 l
# 生成断裂点
& j0 ]4 f+ ~4 d/ ?+ P6 A
breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
' S8 e9 n$ N; ?7 f5 D8 \4 z
# 沿断裂点打断基因组
/ @' ^1 K5 @% ~7 h
self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
8 q8 E( K% X l* i0 [ s. ^
# 模拟克隆时的随机丢失情况
3 A8 a. c% k/ r6 C9 g
clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
4 g) w Y7 |8 D% @& g) d
# 模拟单端测序
4 _5 y) }! s. ]$ B
self.singleread(clonedfragmentList)
+ v" o8 y" U6 q. K7 L( d; g
SeqIO.write(self.readsList, sequencingResult, "fasta")
- X7 w2 @$ b+ {, K, T' T( n: ?" S
. D, H' D, Z9 Q
def pairread(self, clonedfragmentList):
$ g1 a$ {9 }9 h5 D0 b* }
for fragment in clonedfragmentList:
* l6 ?' \0 r3 I) ^3 c
fragment.id = ""
" J! D( a& i- `( W4 a% g
fragment.name = ""
7 n1 ?- j/ o& e4 p/ t! p& R
description = fragment.description[12:].split(",")[0]
# L- d. k4 B$ [7 r3 P
fragment.description = str(self.readsID) + "." + description
2 r; h" B. W* `) ?3 Y) z
readslength = random.randint(self.minreadslength, self.maxreadslength)
! O: \) u. c' w- `. R5 b
self.allreadslength += readslength
. u4 f+ a$ [) m' G
self.readsList.append(fragment[:readslength])
# E$ [, j0 \( v
4 A- I' d/ |% y0 d/ t
readslength = random.randint(self.minreadslength, self.maxreadslength)
8 S# Z2 Q. u- V9 l. Y
self.allreadslength += readslength
2 ~ V% {# y& i" r4 O3 K
2 D/ `; r) M8 U/ S( i: }
fragmentcomplement = fragment.reverse_complement()
" j8 o* e7 {0 r! [" Y. P6 y- k2 H3 x* b
fragmentcomplement.id = ""
* t* S3 h; {, ]/ b. f3 Q0 D4 G6 `3 q& c
fragmentcomplement.name = ""
; x3 V# p9 f0 U
fragmentcomplement.description = str(self.readsID) + "." + description
- }, G% f/ F9 V' L- K: e& b% v
self.readsList.append(fragmentcomplement[:readslength])
4 g3 l& _1 F5 J; b) _8 T
% Y; P p7 Z% c3 K5 K9 b2 ~
self.readsID += 1
* f3 ]" ~1 j+ Z$ [' f) y+ D' v
9 V" _4 `" _; j0 x( c
def pairreadsequencing(self,genomedata, sequencingResult_1, sequencingResult_2):
' p6 e) y A& B" I6 {
for seq_record in SeqIO.parse(genomedata, "fasta"):
5 P! f# B% g& n$ C r& E% `
seqlen = len(seq_record)
( K) H0 B$ D$ K; u! Y
self.genomeLength += seqlen
& x: k" {- d* D/ t5 F1 |9 M0 }; Y" z
for i in range(self.N):
2 o7 F, n$ C' [- O
# 生成断裂点
9 H2 Z- Z8 v8 \) k o! L
breakpoint = self.generatebreakpoint(seqlen, self.averagefragmentlength)
[: I3 [# @, x. ~! q
# 沿断裂点打断基因组
& `- R. O9 F- n
self.breakgenome(seq_record, breakpoint, self.minfragmentlength, self.maxfragmentlength)
; Z% _' @6 C8 l" A% N. |8 H- d! D
# 模拟克隆时的随机丢失情况
4 V4 Y$ _" [) h2 A
clonedfragmentList = self.clonefragment(self.fragmentList, self.cloneRetainprobability)
4 Y% N* [! K7 H; }3 J1 }
# 模拟双端测序
, b8 e" V$ j! B; o' e
self.pairread(clonedfragmentList)
2 I9 O6 R, g& W% I% |" W' X
readsList_1 = [self.readsList
for i in range(len(self.readsList)) if i % 2 == 0]
. `+ J6 u. D# ^$ r
readsList_2 = [self.readsList
for i in range(len(self.readsList)) if i % 2 == 1]
3 Y& A7 l1 t& e# g( C: j4 f6 |7 F
SeqIO.write(readsList_1, sequencingResult_1, "fasta")
% _% [$ u+ t. v" T
SeqIO.write(readsList_2, sequencingResult_2, "fasta")
5 B& Q- c. ?5 N% M% I3 ~
6 h. w: W8 I! m. ]
def resultsummary(self):
( ]& n! g( m3 p4 W" K
print("基因组长度:" + str(self.genomeLength / 1000) + "kb")
; f4 n2 _* {4 C6 A; s
print("片段长度区间:" + str(self.minfragmentlength) + "-" + str(self.maxfragmentlength))
( O6 j, {2 ^# ^: J- [5 m& }* L
print("N值:" + str(self.N))
6 a# }% C# o& x! Q- q1 _2 I! W8 j
print("期望片段长度:" + str(self.averagefragmentlength))
2 T" O5 V% w! M; [: S
print("克隆保留率:" + str(self.cloneRetainprobability))
2 V$ _8 q; w4 B& x5 [, F2 t
print("片段数量:" + str(len(self.fragmentList)))
% b$ q3 i4 b7 T+ t; \) D7 \
print("reads长度:" + str(self.minreadslength) + "-" + str(self.maxreadslength))
$ l' V' j& C$ x/ B' I1 T
print("reads总数量:" + str(len(self.readsList)))
2 V4 h$ R' [( B% `7 D
print("reads总长度:" + str(self.allreadslength / 1000) + "kb")
6 _. `. c% j E6 V
m = self.allreadslength / self.genomeLength
% p& {2 d# H" N
print("覆盖度(m值):" + str(round(m, 5)))
, Y: V+ U) i9 F* Y$ j* {. v( }1 v" g
print("理论丢失率(e^-m):" + str(round(exp(-m), 5)))
. |( e' f! o7 Z: G+ P
print("覆盖率(1-e^-m):" + str(round(1 - exp(-m), 5)))
" }$ O+ n+ j! x. d5 [
# -------------------------------------------主程序-------------------------------------------
9 n X6 a9 X/ Q' x7 S
# 模拟单端测序
8 y" Z& s5 z" F/ p7 s8 v
sequencingObj = Sequencing()
% u# p( X f5 r
sequencingObj.singlereadsequencing("data/NC_025452.fasta", "result/virusSingleRead.fa")
' p& S- c7 [: k2 `& Q$ N$ q* s
sequencingObj.resultsummary()
( H7 D, I, ]* y3 t* J: n9 i
& F7 s. H& v- o8 A; N0 X( t
# 模拟双端测序
* G# k: L. q+ O+ S+ o3 `6 J
sequencingObj = Sequencing()
+ b6 g* m4 ~" r& f2 V. { w' z4 C
sequencingObj.pairreadsequencing("data/GCF_000146045.2_R64_genomic.fna", "result/yeastPairRead_1.fa", "result/yeastPairRead_2.fa")
) y Q5 Q8 P" |: s
sequencingObj.resultsummary()
6 i# g- o8 ?7 t$ ]; ^
; Z; i+ d) l6 r5 a6 ^6 ]5 X
3 o+ Q/ \: O( t! H3 Y
, X& X& K( X8 P7 P6 N _- f
5 A( d1 P/ g1 j- g
数学建模解题思路与方法.pptx
2019-4-21 14:56 上传
点击文件名下载附件
下载积分: 体力 -2 点
117.69 KB, 下载次数: 4, 下载积分: 体力 -2 点
作者:
3477959497
时间:
2019-4-22 10:49
不错。。。。。。。。。。。。。。。。。
4 Z$ D4 h' q/ x* }6 G/ {
欢迎光临 数学建模社区-数学中国 (http://www.madio.net/)
Powered by Discuz! X2.5
