README file for Leucippus software distribution.

Software Goal
=============
Leucippus is a Java program aimed at determining mosaic SNVs and estimating
their allelic frequency at specified sites of interest from amplicon and
capture sequencing data. Mosaic SNV is one that is not uniformly present in
the DNA of all cells sequenced. To make a decision about an SNV being mosaic,
Leucippus estimates background sequencing noise from the other sites
surrounding specified sites of interest. Confidence in conclusion for each
site is reflected in a p-value, which is calculated based on the background
noise. A dedicated functionality displays noise in the data graphically, to
enable data visual inspection and quality control. The software also allows 
constructing long read-fragments from pairs of reads with overlapping 3'-ends.

Leucippus Five Tasks
=====================
 a.  Constructs long reads
        -construction of long read by sliding paired reads and then choosing
          the long read with the best matching window.
           (the sliding is performed from left to right and 
           from right to left(the paired reads are in fastq format, thus 
           one read must be reversed and complemented before the sliding 
           process).
			 
 b.  Creates noise tables that provide the following information 
       for each reference position around the sites of interest:
         - expected base counts
         - specific mismatch counts
         - next position(s) deletion  (number counts)
         - next position(s) insertion (sequence-number counts)
         - total position counts

 c.  Creates graphs (single or comparison):
         - p-value and mutation rate graphs by reading the above table(s).
 
 d.  Making decision about provided sites, whether they are mosaics.

 e.  Create pvalue or mutation rate graphs for a particular site, by reading
	noise table files

1. Compilation
==============
Leucippus is a Java application. It runs on UNIX operating system with java 
installed on it. First the Java extension files (programs) must be compiled
into byte-code using the following directions.

To compile Leucippus users must do the following:

$ cd src
$ javac Leucippus.java

To test "Leucippus.class" while they are in Leucippus root directory, users
must run:

$ java Leucippus

It must print on the screen the following general info about Leucippus
commands:

Program: Leucippus (Mosaics: Site Noise-Estimation, Validation, Identification)
 Version: 0.1.0
 Usage:   Leucippus <command> [options]

  Command: frag           create long reads
           noisetab       create noise table
           graph          create graph
           decide         decide [somatic | germline  | unknown | omitted]
           extract        extract reads that present alternative allele

           config         -cf+suffix <value> [none]



2. Leucippus Commands
=====================
 Leucippus has 4 commands: frag, noisetab, graph, and decide. Each command is
associated with a particular task:
	   
 a.  'frag'     => Constructs long reads from short paired reads.
 b.  'noisetab' => Creates noise tables for sites of interest.
 c.  'graph'    => Creates graphs from noise tables.
 d.  'decide'   => Makes decision whether sites of interest are mosaics.
 e.  'posgraph' => Creates position graphs from noise tables.


Cut offs at each step:

a. 'frag'
min [50] bp overalp
min 75% of matching bases in the overlap (hardcoded)

b. 'noisetab'
min read mapping quality 30 (hardcoded)
min base quality [20]
distance from interval boundaries [-1]

c. 'graph'
min coverage [100]
max range [0.005]

d. 'decide'
min germline AF [0.35]
min coverage [100]
max p-value [0.05]

e. 'posgraph'
min coverage [100]
max range [0.005]

A. Command 'frag'
-----------------

Leucippus constructs long fragments from paired reads which overlap. To find
overlap, reads are sled against each other by one base pair at a time.

    Beginning of read sliding 
    read1 5'----------3'            slide to the right >
    read2          3'----------5'   slide to the left  <

    End of read sliding 
    read1          5'----------3'
    read2 3'----------5'

Matched and mismatched nucleotides in the overlap are counted. Then a
statistical test, using binomial coefficients, is performed to calculate
the probability that such or a larger count of matches for the overlap length
can occur by chance. Here a random chance for a base match is set to 0.25. The
best overlap is the one with the smallest such probability. The pair of reads
with the best overlap of at least 50 nucleotides in length and at least 75% of
matched nucleotides is used to construct a long fragment. Nucleotides in the
overlapping part of both reads are constructed in the following way:
* if nucleotides match, then this nucleotide is assigned at the current
  position, and its quality is equal to the sum of the nucleotide qualities at
  this position in each read.
* if nucleotides mismatch, then the one with the higher quality is assigned at
  the current position, while its quality is reduced by the quality of the
  other nucleotide. If the qualities are equal, one of the two nucleotides is
  chosen randomly and its quality is assigned zero.

Leucippus finds pairs of reads from their names. It recognizes the following
naming conventions for read pairs.
	@name type
	@name:type
	@name/type
	@name
where 'type' is one character (1-9 or a letter) to specify a read in a pair.
For example,

First  read  @HWI-7001311 1
Second read  @HWI-7001311 2

First  read  @read_pair_1:1
Second read  @read_pair_1:3

First  read  @HWI-7001311:50:H9V9KADXX:1:2116:15471:92324/1
Second read  @HWI-7001311:50:H9V9KADXX:1:2116:15471:92324/2

First  read  @HWI-7001311:50:H9V9KADXX:1:2116:15471:92324
Second read  @HWI-7001311:50:H9V9KADXX:1:2116:15471:92324

If other separator (other than 'space', ':', or '/') is used then this can be
stated using a configuration argument -cfseparator.

The usage of the command is as follows:

java Leucippus frag [options] [file1.fq.gz file2.fq.gz]

Options: -o  FILE        output file in gzip (e.g., longreads.fq.gz)
         -o1 FILE        output file in gzip (e.g., remainedreads1.fq.gz)
         -o2 FILE        output file in gzip (e.g., remainedreads2.fq.gz)
         -l              maximum long reads length [500]
Input:   file1.fq.gz     input gzipped file with reads in fastq format
                         (if not provided, input is taken from standard input)
         file2.fq.gz     input gzipped file with reads in fastq format


First fastq file may contain both reads. In such a case there is no need to
provide second fastq file. If no input files are provided, Leucippus expects 
reads in fastq format from standard input (e.g., STDIN). This feature allows
to generate long fragments from aligned .bam file by piping the output 
of 'samtools bam2fq' command to Leucippus (see example below).

Examples:

java Leucippus frag -o long.fq.gz \
     pathtoreads/reads1.fastq.gz \
     pathtoreads/reads2.fastq.gz

java Leucippus frag -o long.fq.gz -o1 remained1.fq.gz -o2 remained2.fq.gz \
     pathtoreads/reads1.fastq.gz \
     pathtoreads/reads2.fastq.gz

java Leucippus frag -o out.fq.gz -l 300 pathtoreads/both_reads.fastq.gz

gzip -cd pathtoreads/both_reads.fastq.gz | java Leucippus frag -o out.fq.gz

samtools bam2fq pathtobam/file.bam | java Leucippus frag -o out.fq.gz

To skip secondary aligments in BAM use the following

samtools view -b -F 0x900 pathtobam/file.bam | samtools bam2fq - | \
     java Leucippus frag -o out.fq.gz


B. Command 'noisetab'
---------------------

Leucippus tabulates sequencing error (e.g., "noise") for sites targeted by
amplicon or capture sequencing. This information is used later to call mosaic
sites. In essence, the information tabulated is a comprehensive count of
nucleotides and small "indels" in each read, covering a particular position.
Therefore, the required inputs are a list of targeted sites, and .bam files
with aligned reads. The table includes all necessary information that is used
to visualize noise and make a decision about sites of interest being mosaic.
	
The tab-delimited interval file specifying the sites of interest is of the
following format:

CHROM	POS	REF	ALT	LEFT	RIGHT
1	11242672	T	A	11242499	11242922
1	72605472	C	T	72605349	72605579
2	15636892	G	C	15636830	15637227

where  
 CHROM is chromosome
 POS   is site of interest
 REF   is expected nucleotide (same as reference nucleotide)
 ALT   is alternative nucleotide (different from REF)
 LEFT  is start of a targeted interval
 RIGHT is end of a targeted interval
                    
All of these columns are required, but additional columns (e.g., with other
information relevant to a particular experiment) are allowed. The distance cut
off is applied to reads covering a particular site of interest and specifies:
  i) the maximum distance from the first aligned position in a long read to the
LEFT  position of the corresponding interval;
 ii) the maximum distance from the last aligned position in a long read to the
RIGHT position of the corresponding interval.
If distance cut off is not specified or negative, then all reads will be 
considered as data source (the default value is -1). The output table is 
of the following format:

chr	pos	x	y	Nuc_ref	Nuc_exp	tot	As	Cs	Ts  Gs
2	5627738	261	184	A	A	8003	7987	9	3   4
2	5627739	262	183	C	C	8003	10	7988	2   3
2	5627740	263	182	C	G	8003	14	7980	3   6
2	5627741	264	181	A	A	8004	7994	3	1   6
2	5627742	265	180	C	C	8004	13	7980	6   5

where  
 chr       is chromosome
 pos       is position
 x         is the distance from the LEFT of the targeted interval
 y         is the distance from the RIGHT of the targeted interval
 Nuc_ref   is the reference nucleotide for this position
 Nuc_exp   is the expected nucleotide (only different from Nuc_ref at
            position for the sites of interest)
 tot       is the total number of by good bases (i.e., passign base quality
 	    and padding cutoffs) at this position
            tot = #A + #C + #T + #G + #Ds, where D represents deletion
 As        is the coverage by As passing all cut-offs
 Cs        is the coverage by Cs passing all cut-offs
 Ts        is the coverage by Ts passing all cut-offs
 Gs        is the coverage by Gs passing all cut-offs

 Noise tables also contain detailed information about mutations, 
deletions (start, length), and insertions (start, length, and sequence) 
of all positions in intervals of targeted site.
Example of additional columns:  
Ds	Ps	cov	D1	D2... D10   I1  I2... I10  DelInfo	InsInfo
0	0	8053	2	0 ... 0	    0	1 ... 0	    C,C         TA,
2	0	8053	0	0 ... 0	    2	0 ... 0	    0	        A,A
0	0	8053	0	0 ... 0	    0	0 ... 0	    0	        0
0	0	8053	0	0 ... 0	    0	0 ... 0     0	        0
0	0	8053	0	0 ... 0	    0	0 ... 0	    0	        0

where 
Ds      is the count of deletions found for this position
Ps	is the count of Ps when padding has been applied
cov     is the total number of reads covering the base; this is used for
	 calculation of indel allele frequency
D1      is the count of reads with one base deletion starting at next position
D2      is the count of reads with two base deletion starting at next position
...
I1      is the count of reads with one base insertion at the next position
I2      is the count of reads with two base insertion at the next position
...
DelInfo is detailed deletion(s) information, for example, "0_1D_5627738_1" is
        read with one base deletion after position 5627738 (at 5627739)
InsInfo is detailed insertion(s) information including the insertion sequence


The usage of the "noisetab" command is as follows:

java Leucippus noisetab [options] file1.bam [file2.bam, ...]

Options:
        -interval FILE  file with intervals of targeted site
        -ref      FILE  FASTA file with reference sequence (could be .gz file)
        -q              base quality cut off 0-100 [20]
        -d              distance cut off [-1]
        -o        FILE  output file
Input:
        file1.bam       a file with aligned reads
        file2.bam       a file with aligned reads


Examples:

java Leucippus noisetab \
     -interval pathtosites/sites.txt \
     -ref      pathtoreference/ref.fa.gz \
     -o        pathtooutput/noisetable.tsv \
     pathtobams/file01_LONG.bam pathtobams/file02_LONG.bam

(example with custom 'samtools' directory)
java Leucippus noisetab \
     -interval pathtosites/sites.txt \
     -ref      pathtoreference/ref.fa.gz \
     -o        pathtooutput/noisetable.tsv \
     -cfsam    pathtosamtools/bin \
     pathtobams/file01_LONG.bam pathtobams/file02_LONG.bam

(example with custom quality cut off, and distance cut off)
java Leucippus noisetab \
     -interval pathtosites/sites.txt \
     -ref      pathtoreference/ref.fa.gz \
     -o        pathtooutput/noisetable.tsv \
     -q 20 -d 10 \
     pathtobams/file01_LONG.bam pathtobams/file02_LONG.bam


C. Command 'graph'
------------------

Leucippus can create graphs for substitution rates and p-values of random
chance, for each nucleotide substitution type (i.e., C to T, etc.). The
graphs are created from the noise table generated at the previous step and
exclude the sites of interest (where Nuc_exp and Nuc_alt are not the same).
Mismatches to the reference base in aligned long reads are assumed to be
noise, either from sample preparation (e.g., amplification) or sequencing. A
mutation rate graph shows the frequency of sites versus rate of a particular
substitution type (e.g., C to T). A p-value graph shows the proportion of
sites having substitution rate larger than a value. For example, for the
expected nucleotide C, substitution type C to T, and substitution frequency
of 0.01 the p-value is the fraction of C-sites having T substitution frequency
larger than 0.01, i.e., #C_sites_with_CtoT_more_0.01/#C_sites. To make the
graphs Leucippus calls an R script passing the necessary parameters to it.

The usage of the "graph" command is as follows:

Leucippus graph [options] -o <prefix> table1.file [table2.file]

 Options:    -type            graph type: pvalue|mutrate [pvalue]
             -coverage        minimum site/position coverage [100]
             -range    DOUBLE maximum range for error [0.005]
             -overlap  INT    use positions in overlapping 3'-ends of reads.
                              The number specifies read length. Only useful 
			      for analysis of amplicon-seq data
             -o               prefix for output files: prefix.pdf,
                                              prefix.fpvtb[1,2].tsv


Second file with noise table is optional. If it is given, graphs from the two
files will be superimposed, with colours representing substitution type and
line style (solid and dashed) representing data from different files. Such
visualization allows for instant comparison of noise level in different data.


Examples:

java Leucippus graph \
	-o pathtooutput/graphortable \
	-type pvalue \
	pathtoinputnoisetable/MosaicNoiseTable.tsv
	
Example of Leucippus command that provides a path for R (graph):
java Leucippus graph \
	-o pathtooutput/graphortable \
	-type pvalue \
	-cfR ~/pathtoR \
	pathtoinputnoisetable/MosaicNoiseTable.tsv

Example of Leucippus command that provides range and coverage arguments:
java Leucippus graph \
        -o pathtooutput/graphortable \
        -coverage 1000 \
        -range 0.03 \
        -type pvalue \
        -cfR pathtoR \
        pathtoinputnoisetable/MosaicNoiseTable.tsv

Example of Leucippus command for comparison (graph):
java Leucippus graph \
	-o pathtooutput/graphortable \
	-type pvalue \
	pathtoinputnoisetable/MosaicNoiseTable1.tsv \
	pathtoinputnoisetable/MosaicNoiseTable2.tsv


D. Command 'decide'
-------------------

Leucippus can decide whether sites of interest are germline or somatic. By
default sites with at least 35% non reference nucleotide coverage are termed
germline.  The decision to call a variant as somatic is done by comparing
the frequency of expected non reference nucleotide with the background
substitution rate of reference to expected nucleotide. The background rate is
derived from all sites in the noise table, excluding sites of interest and
likely germline variants that have more than 10% of non reference nucleotide
coverage. By default, sites with p-value of less than 0.05 are deemed somatic.
Also, by default, sites with read coverage less than 100 reads are deemed as
having insufficient data and omitted from decision process. The remaining sites
are called undetermined, i.e., they are not germline, but frequency of 
expected nucleotide is consistent with background.

The input is a noise table described above for the section for 'noisetab'
command. The output is a tab delimited file in the following format:

CHR	POS	REF	EXP	COV	FREQ(#EXP/COV)	P-VALUE	CONCLUSION
2	5627740	C	G	8107	0.00333045516	0.00131	somatic
2	6912245	T	C	3305	0.00121028744	0.08378	unknown
2	16266496	G	T	19691	0.47036	0.0	germline
2	40177413	C	T	52788	0.05804	0.0	somatic
2	47883415	A	G	152619	0.47365	0.0	germline
2	60074304	A	T	27167	0.07026	0.0	somatic

where:

CHR is chromosome
POS is position
REF is the reference nucleotide
EXP is the expected nucleotide (different than reference)
COV is the site coverage(number of reads found, which include the site)
FREQ is the number of expected nucleotides found, divided by site coverage
P-VALUE is the p-value of the site being somatic
CONCLUSION is the decision made for the site:
         - germline       close to 50% allele frequency
         - omitted        insufficient data (low coverage)
         - somatic        p-value lower than threshold
         - undetermined   not germline and with p-value higher than threshold.

The usage of the "decide" command is as follows:

 Leucippus decide [options] table.file

   Options:  -o           FILE    results
             -coverage            minimum site/position coverage [100]
             -germlineAF  DOUBLE  minimum AF to call variant as germline [0.35]
             -pvalue      DOUBLE  p-value for somatic call[0.05]

Examples:

Example of complete decide command (all parameters have been stated)
java Leucippus decide \
     -coverage 200 \
     -germlineAF 0.35 \
     -o pathtoresults/resultstable.tsv \
     -pvalue 0.03
     pathtoinput/noisetable.tsv

Example of simplest Leucippus decide command with default values for 
"-coverage[100], -germlineAF[0.35], -pvalue[0.05]":

java Leucippus decide \
     -o pathtoresults/resultstable.tsv \
     pathtoinput/noisetable.tsv

E. Command 'posgraph'
--------------------

Leucippus can create graphs for visualizing  nucleotide substitution rates 
and p-values for a site-position or any other position present in multiple noise 
tables. The position graph depicts the behaviour of Allele Frequency of the site
alternative nucleotide found in noise tables. Sinse the site alternative 
nucleotide can be more than one, the posgraph command generates three lines for 
all possible substitutions. For example if the reference in site-position under 
investigation is A, then posgraph (with -type p-value command) will generate 
three pvalue lines for substitutions A to T, A to C, and A to G. Each one of 
the three lines in the plot has different color. The color information is 
provided with the x-title.


Usage: Leucippus posgraph [options] -o <prefix> [table1 table2 .... table(n)]

Options:    -type             graph type: pvalue|mutrate [pvalue]
            -pos              chr(x):position (1-22, X, Y) [chr1:1000000]
            -coverage INT     minimum site/position coverage [100]
            -range    DOUBLE  maximum range for error [0.05]
            -overlap  INT     use positions in overlapping 3'-ends of reads
                               (the number specifies read length and it is
                               only useful for analysis of amplicon-seq data)
            -o                prefix for output files: prefix.pdf

Example:

Example of posgraph command
java Leucippus posgraph \
     -type pvalue \
     -pos 4:153253817 \
     -cfR /pathtoRscript/Rscript \
     -coverage 200 \
     -range 0.05 \
     -o pathtooutput/outposgraph \
     pathtoinput/noisetable1.tsv pathtoinput/noisetable2.tsv ...

F. Command 'extract'
--------------------
Leucippus provides 'extract' command for read retrieval from bam files. User
have to specify position and nucleotide in question. The returned reads will
be those that present the nucleotide in question.
 
The usage of the "extract" command is as follows:

java Leucippus extract [options] file.bam

Options: -o           FILE   results
         -a                  alelle_X [X could be A,C,T,G]
         -p                  chr:position in format according to bam
                               (for examples, chr3:234234, 1:5454635)
example:
java Leucippus extract \
     -o outfile.tsv \
     -a T \
     -p chr15:34563 \
     /pathtobam/file.bam



3. Configuration
================
Leucippus uses softwares samtools and R at various steps of calculations.
To run these softwares it assumes that path to their binaries exists in the
system PATH variable and that bash shell is run by /bin/bash. Additionally,
it generates temporary folder 'Temp' in the current directory, to store
temporary files during execution of 'noisetab' command. Path to the software,
bash shell, and temporary folder can be redefined through configuration
parameters. These parameters can be used with any command.

Configuration Options:
	-cfsam   <path to samtools executable [samtools]>
	-cfR     <path to R executable [R]>
	-cftmp   <path to temporary folder [./Temp]>
	-cfbash  <path to bash shell [/bin/bash]>
	-cfseparator  character.

Examples:

java Leucippus noisetab -interval pathtosites/intervalsites.txt \
    -ref /pathtoref/GenomeRef.tsv \
    -o pathtooutput/Table1.tsv \
    -q 20 -d 0 \
    -cfsam  pathtosamtools/my_samtools \
    -cfbash pathtobash/my_bash \
    -cftmp  pathtotmpfolder/my_tmp \
     pathtoinputbam/longreads.bam

java Leucippus graph \
	-o pathtooutput/graphortable \
	-type   pvalue \
	-cfR    pathtoR/my_R \
	-cfbash pathtobash/my_bash \
	pathtoinputnoisetable/MosaicNoiseTable.tsv

java Leucippus frag -o long.fq.gz \
     -cfseparator #
     pathtoreads/reads1.fastq.gz \
     pathtoreads/reads2.fastq.gz


Please send your comments and suggestions to Abyzov.Alexej@mayo.edu