pywgsim 0.5.2

pywgsim

pywgsim是wgsim短读模拟器的修改版本。

• https://github.com/lh3/wgsim

对wgsim的代码进行了修改，以便将模拟的突变可视化为GFF文件。

该软件包提供了Python包装器和适用于Linux和MacOS的独立编译可执行文件。

安装

pip install pywgsim

PyPI页面：https://pypi.ac.cn/project/pywgsim/

用法

$ pywgsim -h

打印

usage: pywgsim [-h] [-e 0.02] [-D 500] [-s 50] [-N 1000] [-1 70] [-2 70]
               [-r 0.001] [-R 0.15] [-X 0.25] [-S 0] [-A 0.05] [-f]
               genome [read1] [read2]

positional arguments:
  genome                FASTA reference sequence
  read1                 FASTQ file for first in pair (read1.fq)
  read2                 FASTQ file for second in pair (read2.fq)

optional arguments:
  -h, --help            show this help message and exit
  -e 0.02, --err 0.02   the base error rate
  -D 500, --dist 500    outer distance between the two ends
  -s 50, --stdev 50     standard deviation
  -N 1000, --num 1000   number of read pairs
  -1 70, --L1 70        length of the first read
  -2 70, --L2 70        length of the second read
  -r 0.001, --mut 0.001
                        rate of mutations
  -R 0.15, --frac 0.15  fraction of indels
  -X 0.25, --ext 0.25   probability an indel is extended
  -S 0, --seed 0        seed for the random generator
  -A 0.05, --amb 0.05   disregard if the fraction of ambiguous bases higher
                        than FLOAT
  -f, --fixed           each chromosome gets N sequences

与wgsim相比的变化

wgsim原始代码的修改如下

1. 描述wgsim引入的突变的输出以GFF格式生成。
2. 读名中的分隔符字符已从_更改为|。
3. 有一个名为--fixed的新标志，可以为每个染色体生成相同数量的N个读数。

读名命名

读名现在遵循更广泛接受的规范（即NCBI），并允许使用下划线的contigs。此外，更容易进行读名的视觉检查

@NC_002945.4|1768156|1768694|0:0:0|4:0:0|4

固定模式

在wgsim的默认操作中，N个读数被分配以创建所有染色体上的均匀覆盖（较长的染色体获得N的更大比例）。

当启用--fixed模式时，将为每个染色体生成N个读数。引入--fixed模式是为了简化分类器的评估。由于每个输入序列生成相同数量的读数，因此评估分类的质量（即有多少个被正确分类）变得更容易

突变输出

该工具模拟双倍体基因组中的突变。由pywgsim生成的输出将如下所示

##gff-version 3
#
# N=10000 err_rate=0 mut_rate=0.001 indel_frac=0.15000001 indel_ext=0.25 size=500 std=50 len1=70 len2=70 seed=1607013056
#
NC_001416.1    wgsim    snp    89      89      .    +    .    Name=A/R;Ref=A;Alt=R;Type=het
NC_001416.1    wgsim    snp    2825    2825    .    +    .    Name=-/A;Ref=-;Alt=A;Type=het
NC_001416.1    wgsim    snp    3712    3712    .    +    .    Name=G/A;Ref=G;Alt=A;Type=hom
NC_001416.1    wgsim    snp    4622    4622    .    +    .    Name=G/-;Ref=G;Alt=-;Type=hom

解释

• A/R表示具有A/A和A/G等位的杂合突变。
• -/A 表示相对于参考的插入一个 A，类型字段表示杂合突变。
• G/A 表示两个拷贝中都有 G/A 等位的同质突变。
• G/- 表示从参考序列中删除一个 G，类型字段表示同质突变。

模糊代码

A           A           T
C           C           G
G           G           C
T/U         T           A
M         A or C        K
R         A or G        Y
W         A or T        W
S         C or G        S
Y         C or T        R
K         G or T        M
V      A or C or G      B
H      A or C or T      D
D      A or G or T      H
B      C or G or T      V
N    G or A or T or C   N

读取名称约定

读取名称现在是以下形式

   @NC_002945.4|1768156|1768694|0:0:0|4:0:0|4

其中

• NC_002945.4 是生成片段的序列片段名称。
• 1768156 是片段最左侧的位置。
• 1768694 是片段最右侧的位置。
• 0:0:0 是这对中左侧读取的错误、替换和插入/删除的数量。
• 4:0:0 是这对中右侧读取的错误、替换和插入/删除的数量。
• 4 是读取对编号，对每个序列片段是唯一的。

API

到 wgsim 的 C 接口可以通过单个函数调用访问。

from pywgsim import wgsim

wgsim.core(r1="read1.fq", r2="read2.fq", ref="genome.fa", err_rate=0.02, mut_rate=0.001, indel_frac=0.15, indel_ext=0.25, max_n=0.05, is_hap=0, N=100000,  dist=500, stdev=50, size_l=100, size_r=100, is_fixed=0, seed=0)

该函数创建 r1 和 r2 文件。