pIRS (profile based Illumina pair-end Reads Simulator)

Contents
========
1. Introduction
2. Program framework
3. Usage
4. Examples
5. Output file format
6. Notes

1 Introduction
==============

pIRS is a program for simulating paired-end reads from a reference genome.  It
is optimized for simulating reads similar to those generated from the Illumina
platform.

See `INSTALL' for installation instructions.

There are two subcommands: `pirs simulate' and `pirs diploid'.  See section 3
for more details, or run `pirs simulate -h' or `pirs diploid -h'.

2 Program framework
===================
  2.1 Profile Generator
    Six tools are supplied. SOAP2 or BWA, soap.coverage (http://soap.genomics.org.cn/soapaligner.html) 
  are required. The full process is shown in getprofile.sh.example as an example.
    2.1.1 GC%-Depth Profile Stat.
      a). Run soap and soap.coverage to get .depth single file(s). gzip is OK to over it.
      b). Run gc_coverage_bias on all depth single files. You will get gc-depth stat by 1 GC% and other files.
      c). Run gc_coverage_bias_plot on the gc-depth stat file. You'll get PNG plot and a .gc file by 5 GC%.
      d). Manually check the .gc file for any abnormal levels due to the lower depth on certain GC% windows.
    2.1.2 Base-Calling Profile Stat:
      a). Run soap or bwa to get .{soap,single} or .sam file(s).
      b). Run error_matrix_calculator on those file(s). You will get *.{count,ratio}.matrix .
      c). You can use error_matrix_merger to merge several .{count,ratio}.matrix files.
         However, it is up to you to keep the read length matches.
    2.1.3 InDel Profile Stat:
      a). Choose samples with NO polymorphism InDel, such as the Coliphage samples that shipped with Illumina Sequencers.
      b). Run bwa to get .sam/.bam file.
      c). Run indelstat_sam_bam to get the profile.
    2.1.4 Insert size & mapping ratio stat:
      a). Run soap or bwa to get .{soap,single} or .sam file(s).
      b). Run alignment_stator *.
         * alignment_stator cannot stat. mapping ratio for sam files now.
  2.2 Simulator
    Two commands:
    pirs diploid: use for generating diploid genome sequence. Read the input genome sequence and 
                  then simulate SNP, InDel, SV(structure variation) on it. At last, output the 
                  result genome sequence.
    pirs simulate: use for simulating Illumina data, output PE-read file.
    Note:
      a) If you only want to simulate the diploid genome sequence, "pirs diploid" is enough.
      b) If you want to simulate sequencing data of haploid genome, only you need is "pirs simulate".
      c) If you want to simulate sequencing data of diploid genome, you first need to run "pirs diploid" 
         to get the other diploid genome sequence, and then run "pirs simulate" using both the original 
         genome sequence and the previous output sequence as the input.

3 Usage
=======

 pirs <command> [option]
    diploid     generate diploid genome.
    simulate    simulate Illumina reads.

  3.1 pirs diploid
	Usage: pirs diploid [OPTIONS...] REFERENCE

	Simulate a diploid genome by creating a copy of a haploid genome with
	heterozygosity introduced.  REFERENCE specifies a FASTA file containing
	the reference genome.  It may be compressed (gzip).  It may contain multiple
	sequences (scaffolds or chromosomes), each marked with a separate FASTA tag
	line.  The introduced heterozygosity takes the form of SNPs, indels, and
	large-scale structural variation (insertions, deletions and inversions).
	If REFERENCE is '-', the reference sequence is read from stdin, but it must be
	uncompressed.

	The probabilities of SNPs, indels, and large-scale structural variation can be
	specified with the '-s', '-d', and '-v' options, respectively.  You can also
	set the ratio of transitions to transversions (for SNPs) with the '-R' option.

	Indels are split evenly between insertions and deletions. The length
	distribution of the indels is as follows and is derived from panda
	re-sequencing data:
		1bp	64.82%
		2bp	17.17%
		3bp	7.20%
		4bp	7.29%
		5bp	2.18%
		6bp	1.34%

	Large-scale structural variation is split evenly among large-scale insertions,
	deletions, and inversions.  By default, the length distribution of these
	large-scale features is as follows:
		100bp	70%
		200bp	20%
		500bp	7%
		1000bp	2%
		2000bp	1%

	`pirs diploid' does not use multiple threads, even if pIRS was configured with
	--enable-multiple threads.

	OPTIONS:
	  -s, --snp-rate=RATE    A floating-point number in the interval [0, 1] that
				   specifies the heterozygous SNP rate.  Default: 0.001

	  -d, --indel-rate=RATE  A floating-point number in the interval [0, 1] that
				   specifies the heterozygous indel rate.
				   Default: 0.0001

	  -v, --sv-rate=RATE     A floating-point number in the interval [0, 1] that
				 specifies the large-scale structural variation
				 (insertion, deletion, inversion) rate in the diploid
				   genome. Default: 0.000001

	  -R, --transition-to-transversion-ratio=RATIO
				 In a SNP, a transition is when a purine or pyrimidine
				   is changed to one of the same (A <=> G, C <=> T)
				   while a transversion is when a purine is changed
				   into a pyrimidine or vice-versa.  This option
				   specifies a floating-point number RATIO that gives
				   the ratio of the transition probability to the
				   transversion probability for simulated SNPs.
				   Default: 2.0

	  -o, --output-prefix=PREFIX
				 Use PREFIX as the prefix of the output file and logs.
				   Default: "pirs_diploid"

	  -O, --output-file=FILE
				Use FILE as the name of the output file. Use '-'
				   for standard output; this also moves the
				   informational messages from stdout to stderr.

	  -c, --output-file-type=TYPE
				 The string "text" or "gzip" to specify the type of
				   the output FASTA file containing the diploid copy
				   of the genome, as well as the log files.
				   Default: "text"

	  -n, --no-logs          Do not write the log files.

	  -S, --random-seed=SEED Use SEED as the random seed. Default:
				   time(NULL) * getpid()

	  -q, --quiet            Do not print informational messages.

	  -h, --help             Show this help and exit.

	  -V, --version          Show version information and exit.

	EXAMPLE:
	  ./pirs diploid ref_sequence.fa -s 0.001 -d 0.0001 -v 0.000001\
			 -o ref_sequence >pirs.out 2>pirs.err
      
  3.2 pirs simulate

	Usage: ./pirs simulate [OPTION]... REFERENCE.FASTA...

	pIRS is a program for simulating paired-end reads from a genome.  It is
	optimized for simulating reads from the Illumina platform.  The input to
	pIRS is any number of reference sequences. Typically you would just provide
	one FASTA file containing your reference sequence, but you may provide two
	if you have generated a diploid sequence with `pirs diploid', or if you have
	chromosomes split up into multiple FASTA files.  The output of pIRS is two
	FASTQ files containing the simulated paired-end reads, as well as several log
	files.

	Input sequences are assumed to be haploid.  If you instead want to simulate
	reads from a diploid genome, you must give the --diploid option so that
	the diploidy is taken into account when computing coverage.  If you do
	not do this, you will get twice as many reads as you wanted.

	pIRS simulates a normally-distributed insert (fragment) length using the
	Box-muller method.  Usually you want the insert length standard deviation to
	be 1/20 or 1/10 of the insert length mean (see the -m and -v options).
	This program also simulates Illumina sequencing error, quality score and
	GC bias based on empirical distribution profiles. Users may use the default
	profiles in this package, which are generated by large real sequencing data,
	or they may generate their own profiles.

	OPTIONS:
	  -l LEN, --read-len=LEN
			Generate reads having a length of LEN.  Default: 100

	  -x VAL, --coverage=VAL
			 Set the average sequencing coverage (sometimes called depth).
			 It may be either a floating-point number or an integer.

	  -m LEN, --insert-len-mean=LEN
			 Generate inserts (fragments) having an average length of LEN.
			 Default: 180

	  -v LEN, --insert-len-sd=LEN
			 Set the standard deviation of the insert (fragment) length.
			 Default: 10% of insert length mean.

	  -j, --jumping, --cyclicize
			 Make the paired-end reads face away from either other, as
			 in a jumping library.  Default: the reads face towards each
			 other.

	  -d, --diploid
			 This option asserts that reads are being simulated from a
			 diploid genome.  It causes the program to abort if there
			 are not exactly two reference sequences; in addition, the
			 coverage is divided in half, since the two reference
			 sequences are in reality the same genome.  This option
			 is not required to simulate diploid reads, but you must
			 set the coverage correctly otherwise (it will be half
			 as much as you think).

	  -B FILE, --base-calling-profile=FILE, --subst-error-profile=FILE
			 Use FILE as the base-calling profile.  This profile will be
			 used to simulate substitution errors.  Default:
			 PREFIX/share/pirs/Base-Calling_Profiles/humNew.PE100.matrix.gz

	  -I FILE, --indel-error-profile=FILE, --indel-profile=FILE
			 Use FILE as the indel-error profile.  This profile will be
			 used to simulate insertions and deletions in the reads that
			 are artifacts of the sequencing process.  Default:
			 PREFIX/share/pirs/InDel_Profiles/phixv2.InDel.matrix

	  -G FILE, --gc-bias-profile=FILE, --gc-content-bias-profile=FILE
			 Use FILE as the GC content bias profile.  This profile will
			 adjust the read coverage based on the GC content of
			 fragments.  Defaults:
			 PREFIX/share/pirs/GC-depth_Profiles/humNew.gcdep_100.dat,
			 PREFIX/share/pirs/GC-depth_Profiles/humNew.gcdep_150.dat,
			 PREFIX/share/pirs/GC-depth_Profiles/humNew.gcdep_200.dat,
			 depending on the mean insert length.

	  -e FILE, --error-rate=RATE, --subst-error-rate=RATE
			 Set the substitution error rate.  The base-calling profile
			 will still be used, but the average frequency of errors will
			 be changed to RATE.  Set to 0 to disable substitution errors
			 completely.  In that case, the base-calling profile will not
			 be used.  Default: default error rate of base-calling
			 profile.

			 Note: since pIRS parameterizes the error rate by
			 several parameters, it is very difficult to determine exactly
			 what needs to be done to make the error rate be a given
			 value.  We try to adjust the probabilities of getting each
			 quality score in order to accomodate the user-supplied error
			 rate.  However, depending on your input sequences, the actual
			 error rate simulated by pIRS could be off by 20% or more.
			 Please check the informational output to see the final error
			 rate that was actually simulated.

	  -A ALGO, --substitution-error-algorithm=ALGO, --subst-error-algo=ALGO
			 Set the algorithm used for simulating substitition errors.
			 It may be set to the string "dist" or "qtrans".  The
			 default is "qtrans".

			 The "dist" algorithm looks up the substitution error rate
			 for each base pair based on the current cycle and the true
			 base.  This lookup produces a quality score and a called base
			 that may or may not be the same as the true base.  In the
			 base-calling profile, the matrix we use is marked as the
			 [DistMatrix].

			 The "qtrans" algorithm is a Markov-chain model based on the
			 previous quality score and current cycle.  The next quality
			 score is looked up with a certain probability based on these
			 parameters.  The matrix used for this is marked as
			 [QTransMatrix] in the base-calling profile. Then, the the
			 DistMatrix is used to find a called base for the quality score.
			 The DistMatrix is also used to call the base in the first
			 cycle.

	  -M MODE, --mask=MODE, --eamss=MODE
			 Use the EAMSS algorithm for masking read quality.  MODE may be
			 the string "quality" or "lowercase".  The EAMSS algorithm
			 identifies low-quality, GC-rich regions near the ends of reads.
			 "quality" mode will change the quality scores on these
			 regions to (2 + quality_shift), while "lowercase" mode
			 will change the base pairs to lower case, but not change
			 the quality values.  Default: Do not use EAMSS.

	  -Q VAL, --quality-shift=VAL, --phred-offset=VAL
			 Set the ASCII shift of the quality value (usually 64 or 33 for
			 Illumina data).  Default: 33

	  --no-quality-values
	  --fasta
			 Do not simulate quality values.  The simulated reads will be
			 written as a FASTA file rather than a FASTQ file.
			 Substitution errors may still be done; if you do not want
			 to simulate any substition errors, provide --error-rate=0 or
			 --no-substitution-errors.

	  --no-subst-errors
	  --no-substitution-errors
			 Do not simulate substitution errors.  Equivalent to
			 --error-rate=0.

	  --no-indels
	  --no-indel-errors
			 Do not simulate indels.  The indel error profile will not be
			 used.

	  --no-gc-bias
	  --no-gc-content-bias
			 Do not simulate GC bias.  The GC bias profile will not be
			 used.

	  -o PREFIX, --output-prefix=PREFIX
			 Use PREFIX as the prefix of the output files.  Default:
			 "pirs_reads"

	  -c TYPE, --output-file-type=TYPE
			 The string "text" or "gzip" to specify the type of
			 the output FASTQ files containing the simulated reads
			 of the genome, as well as the log files.  Default: "text"

	  -z, --compress
			 Equivalent to -c gzip.

	  -n, --no-logs, --no-log-files
			 Do not write the log files.

	  -S SEED, --random-seed=SEED
			 Use SEED as the random seed. Default:
			 time(NULL) * getpid().  Note: If pIRS was not compiled with
			 --disable-threads, each thread actually uses its own random
			 number generator that is seeded by this base seed added to
			 the thread number; also, if you need pIRS's output to be
			 exactly reproducible, you must specify the random seed as well
			 as use only 1 simulator thread (--threads=1, or configure
			 with --disable-threads, or run on system with 4 or fewer
			 processors).

	  -t, --threads=NUM_THREADS
			 Use NUM_THREADS threads to simulate reads.  This option is
			 not available if pIRS was compiled with the --disable-threads
			 option.  Default: number of processors minus 2 if writing
			 uncompressed output, or number of processors minus 3 if
			 writing compressed output, or 1 if there are not this many
			 processors

	  -q, --quiet    Do not print informational messages.

	  -h, --help     Show this help and exit.

	  -V, --version  Show version information and exit.


4 Examples
==============
  4.1 Simulating a diploid genome sequence.
    Example command line:
    pirs diploid Human_ref.fa -s 0.001 -R 2 -d 0.00001 -v 0.000001 -o Human >simulate_seq.o 2>simulate_seq.e
    Output files:
      a) Human.snp.indel.invertion.fa:   another diploid genome sequence.
      b) Human_indel.lst:                InDel information list.
      c) Human_snp.lst:                  SNP information list.
      d) Human_invertion.lst:            invertion information list.
      e) simulate_seq.o, simulate_seq.e: records of the program running information.
  4.2 Simulate Illumina paired-end reads from a haploid genome.
    Example command line:
    pirs simulate Human_ref.fa -m 170 -l 90 -x 5 -v 10 -o Human >simulate_170.o 2>simulate_170.e
    Output files:
      a) Human_90_170_1.fq, Human_90_170_2.fq:    the paired-end read files.
      b) Human_90_170.error_rate.distr:           the error distribution file.
      c) Human_90_170.insert_len.distr:           the insert length distribution file.
      d) Human_90_170.read.info:                  information about every simulated reads
      e) simulate_170.o, simulate_170.e:          records of the program running information.
  4.3 Simulate Illumina paired-end reads from a diploid genome.
    Example command line:
    pirs simulate --diploid Human_ref.fa Human.snp.indel.invertion.fa.gz -m 800 -l 70 -x 5 -v 10 
                            -o Human >simulate_800.o 2>simulate_800.e
    Output files:
      a) Human_70_800_1.fq, Human_70_800_2.fq:   the pair-end read files.
      b) Human_70_800.error_rate.distr:          the error distribution file.
      c) Human_70_800.insertsize.distr:          the insert size distribution file.
      d) Human_70_800.read.info.gz:              record the information of every reads.
      e) simulate_800.e, simulate_800.o:         records of the program running information.

5 Output file format
====================
  a) *.fq/*.fq.gz

    FASTQ files containing the simulated reads.  The files are given names
    similar to PREFIX_70_800_1.fq, where PREFIX is the prefix provided by the -o
    option (default: "pirs_reads"), 70 is the read length, 800 is the mean
    insert length, and 1 means the file for read 1 of the read pairs.  These
    files will be written as GZIP files with the .gz suffix if you provide the
    `-c gzip' option.

      @read_800_21/1
      ACGGAAAAGTTACGCTATCGCATGCGTGTAAGAACACTGCTCCTACGCCCATTTTATCGATGGCGCCCAG
      +
      egcggdggfgfgggggfeggggYbcgegfgggbggg^e]egfegggfbSeggdggegg`^eJgggcbEeb
      @read_800_22/1
      CACGGGGGGACTTTATTTAATGAGCGGCTGTAACTTGGTCCGTCGTTTGAGAGGGGACACCTCATATGAT
      +
      gggggegcgeggggcgfcgc_gf_ggfefcgVgegcfcdgf`geggdd[ge`ggafeggggdgdgee^gg
    
  b) *.read.info

    Information about each simulated read.
    
    Column 1: reads ID.
    Column 2: Name of the reference file; this is mainly provided so that you
              tell which chromosome set a read came from if you simulate reads
              from a diploid genome produced with the `pirs diploid' command.
    Column 3: FASTA/FASTQ sequence ID of the contig/scaffold/chromosome.
    Column 4: position(1-based) in chromosome.
    Column 5: "+" forward direction; "-" reverse direction.
    Column 6: real insert size.
    Column 7: length of read-end masking by EAMSS algorithm.
    Column 8: read-position(1-based) of substitution error and raw base(->)error base.
    Column 9: read-position(1-based) of insertion error and the base of insertion.
    Column 10: read-position(1-based) of deletion error and the base of deletion

  c) *_snp.lst
    For example:
    
      I	3	3	G	A 
      I	45342	45355	C	T 
      I	104775	104680	C	T
      .....
    Column 1: chromosome sequence ID.
    Column 2: position(1-based) of SNP in reference chromosome.
    Column 3: position(1-based) of SNP in simulated diploid chromosome.
    Column 4: base of reference chromosome sequence.
    Column 5: the base of SNP.

  d) *_indel.lst
    For example:
      I  - 3 1 C
      IV + 5 2 AC
      .....
    Column 1: chromosomesequence ID.
    Column 2: "-" Deletion; "+" Insertion.
    Column 3: position(1-based) of InDel in the reference chromosome sequence.
              For deletions, it is the position of the first deleted base.  For
              insertions, it is the position of the reference base corresponding
              to the base directly before the insertion. 
    Column 4: position(1-based) of InDel in the diploid sequence.  For
              deletions, it is the position of the diploid sequence base
              corresponding to the reference base directly before the deletion.
              For insertions, it is the position of the first inserted base.
    Column 5: length of InDel (always positive).
    Column 6: sequence of the InDel.

      I - 3 3 1 C
      Ref: a t G c a // 3 is the position of G base in the chromosome base.
         : a t G - a // following the G base is a deletion of "C" base.
      IV + 5 5 2 A C
      Ref: t t g c A - - g t t// 5 is the position of A base in the chromosome base.
         : t t g c A A C g t t// following the A base is a insertion of "AC" bases.

  e) *_inversion.lst
    For example:
      I	50191	50195	100 
      I	948984	948903	200
    Column 1: chromosome sequence ID.
    Column 2: position(1-based) of beginning of inversion in the reference
              chromosome sequence.
    Column 3: position(1-based) of beginning of inversion of the diploid
              chromosome sequence.
    Column 3: length of inversion.

6 Notes
======
  a) pIRS does not simulate reads containing "N" char.  If your input reference
     contains the "N" character, reads generated in this in this region will be
     discarded.
  b) When running a simulation without the --no-subst-errors option, the maximum
     length of the simulated reads depends on the number of cycles recorded in
     the base-calling profile.  The user must set the read length to no more
     than half the number of cycles recorded in the base-calling profile.
  c) When masking quality values, the program uses the same EAMSS algorithm
     from CASAVA v1.8.0.  This is only done if the --eamss option is supplied.
  d) The program parses one chromosome at a time.  Reads are evenly distributed
     among the chromosomes.
  e) To re-do a simulation and get the exact same results, you must set the
     random seed with the -S parameter.  In addition, you must use the
     single-threaded version of pIRS, or else set --threads=4 (for only one
     simulator thread).
  f) pIRS will shift quality values when the substitution-error rate setting by
     user is different from that in profile. The quality score in the output
     data range from 2 to the maximum score of profile.  You can find the real
     substitution-error rate of output data in file *.error_rate.distr.
  
  For update & support, please refer to http://code.google.com/p/pirs/ .