wgsim说明

wgsim 就是这样一个软件，它是由开发了BWA等软件大牛 Li heng 写的基因组 转成  短序列的软件。

安装wgsim

~/app$ git clone https://github.com/lh3/wgsim.git && cd wgsim
~/app$ gcc -g -O2 -Wall -o wgsim wgsim.c -lz -lm
~/app$ sudo ln -s `pwd`/wgsim /usr/local/sbin
~/app$ wsgim -h

Program: wgsim (short read simulator)
Version: 0.3.1-r13
Contact: Heng Li <lh3@sanger.ac.uk>

Usage:   wgsim [options] <in.ref.fa> <out.read1.fq> <out.read2.fq>

Options: -e FLOAT      base error rate [0.020]
         -d INT        outer distance between the two ends [500]
         -s INT        standard deviation [50]
         -N INT        number of read pairs [1000000]
         -1 INT        length of the first read [70]
         -2 INT        length of the second read [70]
         -r FLOAT      rate of mutations [0.0010]
         -R FLOAT      fraction of indels [0.15]
         -X FLOAT      probability an indel is extended [0.30]
         -S INT        seed for random generator [-1]
         -A FLOAT      disgard if the fraction of ambiguous bases higher than FLOAT [0.05]
         -h            haplotype mode

模拟基因组短序列数据

使用所有参数为默认值，将大肠杆菌基因组数据转换为PE250 fastq 格式数据。

~/data$ wget https://raw.github.com/ecerami/samtools_primer/master/tutorial/genomes/NC_008253.fna
~/data$ wgsim -S11 -1250 -2250 NC_008253.fna reads_1.fastq reads_2.fastq
~/data$ head -8 reads_1.fastq

@gi|110640213|ref|NC_008253.1|_3151106_3151623_4:0:0_5:1:0_0/1
AACATAAGTGGTATTAATTCCCCAAGATTCAAGATTACGAATAGTATTATCCGCAAAAATATCATCACCTACTTTCGTCAGCATCAGGACTTTTGAATTCAATTTAGCCGCCGCCACCGCTTGATTAGCACCTTTGCCACCACATCCGATTTTGAAGGCAGGTGCTTCCAGAGTTTCTCCTTCTTTAGGCATCTGATAAGTGTAAGTAATGAGATCCACCATATTGGAACCAATAAGTGCAATGTCCATT
+
2222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222
@gi|110640213|ref|NC_008253.1|_2429297_2429873_8:0:0_7:0:0_1/1
ACGGTATCACCATGATTGCCCGTTCCGTCAACAGCATGGGGCTGGTCATTATAGGCGGTGGATCGCTTGAAGAAGCGTTAACTGAACTGGAAACCGGACGCGGCGACGCGGTGGTGGTGCTGGAAAACGAACTGCATCGTCACGCTTGCGCTACCCGCGTGAATGCTGCGCTGGCTAAAGCGCCGCTGGTGATTGTGGTTGACCATCAACGCACAGCGATTATGGAAAACGCTCATCTGGTACTATCCGC
+
222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222

 

 

下载安装：

git clone  https://github.com/lh3/wgsim.git
gcc -g -O2 -Wall -o wgsim wgsim.c -lz -lm

五、软件使用：

软件比较简单，输入文件为基因组序列，fasta格式，输出为illumina的fastq格式，这些格式我们在前面都介绍过。然后是一些选项。
-e 是错误率，默认是0.02
-d reads两头的距离，也就是插入片段长度，默认250bp，注意插入片段本身是包含reads长度的，而不是reads之间的距离
-s 是-d插入片段的偏差，默认是20，也就是-d的值加减20，我们知道插入片段长度并不是固定的，而是一个范围
-N 是测序的层数，控制输出数据量
-1 是reads1长度，默认70bp
-2 是reads2长度，默认70bp
-r 突变率
-R -X 都是调整indels的
-h 是单倍体模式
下面我们来运行一下
wgsim 参考序列 reads1 reads2 这里插入片段我们选择500bp，偏差-s在50，reads长度-1 -2为100bp，二者可以不一样，其余默认。

六、使用案例：

wgsim ref.fna reads1.fq reads2.fq -d 500 -s 50 -1 90 -2 90

七、注意事项：

1、模拟出的reads质量值是无法更改的，都是“I”，如果程序用到reads的质量值模拟数据就会有问题。
2、不支持Mate-pair文库，就是即使把-d设置微6K，那么它是不能像实际过程中发生环化的，两条reads的方向和小片段还是一样的。

介绍
============

Wgsim是从参照基因组中模拟序列的小工具。
它能够模拟二倍体基因组与SNP和插入/缺失（INDEL）
多态性，并能够模拟均匀替代测序错误的reads。
它不产生INDEL测序错误，但是这可能是部分地
通过模拟INDEL多态性补偿。

Wgsim输出是模拟多态性，并写入真正的reads坐标
以及在reads名称的多态性和测序错误的数量。
我们可以wgsim_eval.pl自带的包评估映射的准确性或SNP caller。


编译
===========

GCC -g -O2 -Wall -o wgsim wgsim.c -lz -lm


评估
==========

仿真与评估
-------------------------

仿真命令行：

  wgs​​im -Nxxx -1yyy -d0 -S11 -e0 -rzzz hs37m.fa YYY-zzz.fq的/ dev / null的

其中yyy是read长，zzz为错误率和$ XXX * $ YYY =10000000。
默认情况下，多态性的15％以上是插入缺失和它们的长度是从绘制
几何分布密度0.7 * 0.3 ^ {L-1}。

评估命令行：

  wgs​​im_eval.pl独特aln.sam | wgs​​im_eval.pl alneval -g 20

在'-g'选项可以和映射器来改变。


系统
------

GCC: 4.1.2
CPU: AMD Opteron 8350 @ 2.0GHz
Mem: 128GB


实际操作：

使用过的命令：

../seqtk_concat   out_1.fq  out_2.fq | head -n 8

../seqtk_concat   out_1.fq  out_2.fq  >  ../wgsim.fq   #双端序列合并

seqtk_names  wgsim.fq  | cut -f1,2 -d'_'| cut -f3 -d'|' | perl -ane 's/chr\S+/9606/; print' | perl -ne 'print qq{$.\t$_}' | tabtk_bins - 1 >tabtk_stats 

seqtk_names  wgsim.fq  | cut -f1,2 -d'_'| cut -f3 -d'|' | perl -ane 's/chr\S+/9606/; print' | perl -ne 'print qq{$.\t$_}' | tabtk_bins - 1 >tabtk_stats.tsv

cut -f1,2 tabtk_stats.tsv | tabtk_decorate  /biostack/database/taxonomy/node-name.tsv   0  -  > tab.xls
/*

解释：
perl -ane 's/chr\S+/9606/; print'    将每一行中包含chr并且其后非空格部分（比如chrab和chrbc）替换为9606
perl -ne 'print qq{$.\t$_}'  打印每一行的行号和内容;

可简化为：

|perl -ane 's/chr\S+/9606/; print qq{$.\t$_}'  或 perl -ane 's/chr\S+/9606/; print qq{$.\t$_}'

*/

cut说明：

cut 不就是『切』吗？没错啦！这个指令可以将一段讯息的某一段给他『切』出来～ 处理的讯息是以
『行』为单位喔！

常用的参数：