Join GitHub today
GitHub is home to over 28 million developers working together to host and review code, manage projects, and build software together.Sign up
A tool for confirmation of mosaic variants from targeted sequencing
Fetching latest commit…
Cannot retrieve the latest commit at this time.
|Type||Name||Latest commit message||Commit time|
|Failed to load latest commit information.|
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  bp overalp min 75% of matching bases in the overlap (hardcoded) b. 'noisetab' min read mapping quality 30 (hardcoded) min base quality  distance from interval boundaries [-1] c. 'graph' min coverage  max range [0.005] d. 'decide' min germline AF [0.35] min coverage  max p-value [0.05] e. 'posgraph' min coverage  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  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  -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  -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  -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, -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  -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