it's a simulation for wes , aimed at assessing tumor heterogeneity analysis software
该脚本用于模拟wes/wgs的fastq格式的测序数据.脚本最初的目的是为了模拟肿瘤测序,期望用来比较肿瘤异质性分析软件的准确性.因为肿瘤测序需要很高的深度,常常使用外显子组测序,所以一开始以wes测序为设计出发点,后来加入wgs.
这个工具考虑了wes/wgs两种测序技术,模拟近乎所有种类变异型(以公式的方式表示),并能按照不同细胞比例组织样模拟测序数据.
下表列出了脚本的特点
| 功能 | 有无 | 定义方式 |
|---|---|---|
| 倍型 | √ | 配置 |
| 突变型 | √ | 配置 |
| 宏基因组 | × | × |
| read质量分布 | √ | 文件 |
| 文库片段长度 | √ | 配置 |
| read长度 | √ | 文件 |
| 插入,缺失,替换错误 | √ | 配置 |
| 错配 | √ | 文件 |
| 捕获区域 | √ | 文件 |
| 各种突变 | √ | 文件 |
| 质量错误率转换 | √ | 配置+文件 |
| PE/SE | √ | 配置 |
| ME | × | × |
| illumina | √ | × |
| 其它平台 | × | × |
python 3.7(3.5可能支持), matplotlib
共有6个文件: simu.py, read.py, mutation.py, sequence.py, basic.py, profile.ini
simu.py是主操作工程文件,profile.ini是配置文件,在其中设置主要参数,其它几个文件顾名可知功能,是部分功能组件.
- 获取各种输入文件,
- 设置配置文件
- 命令行输入执行
- 参考基因组文件(fna)
- 捕获区域文件(bed)
- 原始实验数据(fastq) (从已有的fastq文件中获取一部分测序平台信息)
- 突变设置文件
1-3文件都是现成的数据,4文件不是vcf等常用格式文件,因为这些文件大多数从结果推原因,不是原因推结果,定义突变的能力有限,所以这里采用新定义的格式,需要使用者按照自己的目的进行突变设置.格式说明请见(突变设置)
配置采用INI格式文件,详情说明见文件中的注释和模拟过程,配置主要设置了输入文件的路径和各种参数值,以及模式切换
脚本主要的设置都在配置文件中,命令行的操作比较简便易懂,一个命令就是一个操作,没有过多复杂的可选参数,分为5个主要部分,2个辅助部分,如下表
| 命令行 | 功能 |
|---|---|
| ref | 格式化参考基因组 |
| dep | 从depth文件中获得捕获区域 |
| reg | 格式化捕获区域 |
| qph | 格式化参考fastq |
| mut | 格式化突变设置文件 |
| read | 输出fastq下机文件 |
| seq | 依照格式化的参考基因组和捕获 区域,获得可捕获区域的序列 |
| view | 参考基因组浏览(需要格式化)和显示深度与外显子区间 |
ref,reg,qph,mut,read是4个主要部分,依次执行即可.seq可以查看正常基因组的外显子格式化后的序列(加上了侧翼序列的),view是辅助设计突变文件或查看深度与区间的匹配情况的,dep是在没有现成区域文件时,用测序深度文件获得区间的.
python simu.py -ref
[python simu.py -dep depthfile avg_depth]
python simu.py -reg
python simu.py -qph
[python mutation.py VCFfile] # VCF转换格式
python simu.py -mut
python simu.py -read [-R label]注:[]是可选的
内容主要分为三列: 基因型, 缺失位点, 插入序列,三者用制表符分隔开.其中缺失位点再分为染色体位置,缺失开始位点,缺失结束位点三列.再加上可选的描述部分一共7列.
-
基因型
*表示突变,.表示没突变,各个单倍型直接用/分开,比如genotype description */. 一个变异,另一个正常 */* 纯合突变
-
缺失位点
染色体,开始,结尾:表示染色体的开始位点到结尾位点的区间缺失(包括两点),之间用制表符分隔开.如果位点没缺失只是插入,开始/结尾其中一个用.替代即可,不可全都为.,比如chr start end description 1 100 100 1号染色体100号位点缺失 1 100 200 1号染色体100-200缺失 1 100 . 1号染色体100位点后插入 1 . 100 1号染色体100位点后插入(两种写法) 1(.-.) 无效 -
插入序列
染色体(开始-(indel开始-indel结尾)indel序列-结尾)*拷贝数: 表示插入了染色体开始到结束,并且这个区间的indel开始到结尾发生了indel变异,插入的区间是有多个(拷贝数个)这样的序列拷贝串联的.如果倒位,拷贝数为负数即可.如果有多段拼接,用,区分开.比如:insertion description 1(100-200)*2 插入了串联的两条1号染色体100-200序列 1(100-(150,.)AA-200)*1 插入了1号染色体100-200区间,其中的150位点发生了indel插入突变 1(100-200)*-2 插入了策略的两条倒转的200-100+200-100序列 (1(100-200)*2,2(100-200)*1)*1 插入了2个1号的100-200和1个二号的100-200 -
完整的一行
genotype chr begin end insertion description */. 1 100 200 2(100-200)*2 杂合突变,1号染色体100-200缺失,替换为2号染色体100-200串联2个的序列