Genome Analysis Wiki - User contributions [en]

BamUtil: convert

2011-10-04T18:47:48Z

Pha: /* BAM File Recovery */ better grammar and more text

[[Category:BamUtil|convert]]
[[Category:BAM Software]]
[[Category:Software]]

= Overview of the <code>convert</code> function of <code>bamUtil</code> =
The <code>convert</code> option on the [[bamUtil]] executable reads a SAM/BAM file and writes it as a SAM/BAM file.

The executable converts the input file into the format of the output file.

It has options to allow for the conversion of the sequence to/from '=' from/to the actual bases by using the reference sequence.

If you want to convert a BAM file to a SAM file, just call:
<pathToExe>/bam --in <bamFile>.bam --out <newSamFile>.sam
Don't forget to put in the paths to the executable and your test files.

= Parameters =
<pre>
Required Parameters:
--in : the SAM/BAM file to be read
--out : the SAM/BAM file to be written
Optional Parameters:
--refFile : reference file name
--noeof : do not expect an EOF block on a bam file.
--params : print the parameter settings
--recover : attempt to recover the input bam file.
Optional Sequence Parameters (only specify one):
--seqOrig : Leave the sequence as is (default & used if reference is not specified).
--seqBases : Convert any '=' in the sequence to the appropriate base using the reference (requires --ref).
--seqEquals : Convert any bases that match the reference to '=' (requires --ref).
</pre>

== Sequence Representation Parameters ==
The sequence parameters options specify how to represent the sequence if the reference is specified (refFile option). If the reference is not specified or seqOrig is specified, no modifications are made to the sequence. If the reference and seqBases is specified, any matches between the sequence and the reference are represented in the sequence as the appropriate base. If the reference and seqEquals is specified, any matches between the sequence and the reference are represented in the sequence as '='.

=== Examples ===
ExtendedCigar: SSMMMDDMMMIMNNNMPMSSS
Sequence: AATAA CTAGA T AGGG
Reference: TAACCCTA ACCCT A
Sequence with Orig: AATAACTAGATAGGG
Sequence with Bases: AATAACTAGATAGGG
Sequence with Equals: AA======G===GGG

ExtendedCigar: SSMMMDDMMMIMNNNMPMSSS
Sequence: AATGA CTGGA T AGGG
Reference: TAACCCTA ACCCT A
Sequence with Orig: AATGACTGGATAGGG
Sequence with Bases: AATGACTGGATAGGG
Sequence with Equals: AA=G===GG===GGG

ExtendedCigar: SSMMMDDMMMIMNNNMPMSSS
Sequence: AAT=A CT=GA T AGGG
Reference: TAACCCTA ACCCT A
Sequence with Orig: AAT=ACT=GATAGGG
Sequence with Bases: AATGACTGGATAGGG
Sequence with Equals: AA======G===GGG

ExtendedCigar: SSMMMDDMMMIMNNNMPMSSS
Sequence: AA=== ===G= = =GGG
Reference: TAACCCTA ACCCT A
Sequence with Orig: AA======G===GGG
Sequence with Bases: AATAACTAGATAGGG
Sequence with Equals: AA======G===GGG

= BAM File Recovery =

A BAM file that has been corrupted or truncated due to a copy or disk problem can often be partially recovered.

Both the BGZF format and binary BAM format have enough information to scan forward and resynchronize the input data. While some data will be lost, substantial recovery can often be done.

When a file has bad blocks in it, normal copy commands (cp) will truncate the file at the point of disk read failure. To recover the maximum amount of data possible, use the dd command with the conv=noerror option.

So a normal use case for recovery would look this this:

<pre>
# dd if=brokenbamfile.bam of=/tmp/brokenbamfile1.bam conv=noerror bs=4k
# bam convert --recover --in /tmp/brokenbamfile1.bam --out /tmp/brokenbamfilerecovered.bam
</pre>

Note, you will of course need to output the result file to a known good filesystem.

Currently, no statistics are printed as far as how many BAM records are recovered, but subsequent tests can readily be done on the resulting file to determine the quality of recovery.

In real cases, we have recovered better than 94% of reads from a set of severely damaged files (numerous 64K chunks of a RAID were lost), and better than 99.9% recovery from a moderately damaged file (3 disk pages were corrupt).

= Usage =
./bam convert --in <inputFile> --out <outputFile.sam/bam/ubam (ubam is uncompressed bam)> [--refFile <reference filename>] [--seqBases|--seqEquals|--seqOrig] [--noeof] [--params]

= Return Value =
Returns the SamStatus for the reads/writes.

= Example Output =
<pre>
Number of records read = 10
Number of records written = 10
</pre>

BamUtil: convert

2011-10-04T18:44:16Z

Pha: add recovery section

[[Category:BamUtil|convert]]
[[Category:BAM Software]]
[[Category:Software]]

= Overview of the <code>convert</code> function of <code>bamUtil</code> =
The <code>convert</code> option on the [[bamUtil]] executable reads a SAM/BAM file and writes it as a SAM/BAM file.

The executable converts the input file into the format of the output file.

It has options to allow for the conversion of the sequence to/from '=' from/to the actual bases by using the reference sequence.

If you want to convert a BAM file to a SAM file, just call:
<pathToExe>/bam --in <bamFile>.bam --out <newSamFile>.sam
Don't forget to put in the paths to the executable and your test files.

= Parameters =
<pre>
Required Parameters:
--in : the SAM/BAM file to be read
--out : the SAM/BAM file to be written
Optional Parameters:
--refFile : reference file name
--noeof : do not expect an EOF block on a bam file.
--params : print the parameter settings
--recover : attempt to recover the input bam file.
Optional Sequence Parameters (only specify one):
--seqOrig : Leave the sequence as is (default & used if reference is not specified).
--seqBases : Convert any '=' in the sequence to the appropriate base using the reference (requires --ref).
--seqEquals : Convert any bases that match the reference to '=' (requires --ref).
</pre>

== Sequence Representation Parameters ==
The sequence parameters options specify how to represent the sequence if the reference is specified (refFile option). If the reference is not specified or seqOrig is specified, no modifications are made to the sequence. If the reference and seqBases is specified, any matches between the sequence and the reference are represented in the sequence as the appropriate base. If the reference and seqEquals is specified, any matches between the sequence and the reference are represented in the sequence as '='.

=== Examples ===
ExtendedCigar: SSMMMDDMMMIMNNNMPMSSS
Sequence: AATAA CTAGA T AGGG
Reference: TAACCCTA ACCCT A
Sequence with Orig: AATAACTAGATAGGG
Sequence with Bases: AATAACTAGATAGGG
Sequence with Equals: AA======G===GGG

ExtendedCigar: SSMMMDDMMMIMNNNMPMSSS
Sequence: AATGA CTGGA T AGGG
Reference: TAACCCTA ACCCT A
Sequence with Orig: AATGACTGGATAGGG
Sequence with Bases: AATGACTGGATAGGG
Sequence with Equals: AA=G===GG===GGG

ExtendedCigar: SSMMMDDMMMIMNNNMPMSSS
Sequence: AAT=A CT=GA T AGGG
Reference: TAACCCTA ACCCT A
Sequence with Orig: AAT=ACT=GATAGGG
Sequence with Bases: AATGACTGGATAGGG
Sequence with Equals: AA======G===GGG

ExtendedCigar: SSMMMDDMMMIMNNNMPMSSS
Sequence: AA=== ===G= = =GGG
Reference: TAACCCTA ACCCT A
Sequence with Orig: AA======G===GGG
Sequence with Bases: AATAACTAGATAGGG
Sequence with Equals: AA======G===GGG

= BAM File Recovery =

A BAM file that has been corrupted, or truncated due to a copy or disk problem can often be partially recovered.

Both the BGZF format and binary BAM format have enough information to scan forward and resynchronize the input data. While some data will be lost, substantial recovery can often be done.

When a file has bad blocks in it, normal copy commands (cp) will truncate the file at the point of disk read failure. To recover the maximum amount of data possible, use the dd command with the conv=noerror option.

So a normal use case for recovery would look this this:

<pre>
# dd if=brokenbamfile.bam of=/tmp/brokenbamfile1.bam conv=noerror bs=4k
# bam convert --recover --in /tmp/brokenbamfile1.bam --out /tmp/brokenbamfilerecovered.bam
</pre>

Note, you will of course need to output the result file to a known good filesystem.

Currently, no statistics are printed as far as how many BAM records are recovered, but subsequent tests can readily be done on the resulting file to determine the quality of recovery.

= Usage =
./bam convert --in <inputFile> --out <outputFile.sam/bam/ubam (ubam is uncompressed bam)> [--refFile <reference filename>] [--seqBases|--seqEquals|--seqOrig] [--noeof] [--params]

= Return Value =
Returns the SamStatus for the reads/writes.

= Example Output =
<pre>
Number of records read = 10
Number of records written = 10
</pre>

BamUtil: convert

2011-10-04T18:32:26Z

Pha: /* Parameters */ add recover option

[[Category:BamUtil|convert]]
[[Category:BAM Software]]
[[Category:Software]]

= Overview of the <code>convert</code> function of <code>bamUtil</code> =
The <code>convert</code> option on the [[bamUtil]] executable reads a SAM/BAM file and writes it as a SAM/BAM file.

The executable converts the input file into the format of the output file.

It has options to allow for the conversion of the sequence to/from '=' from/to the actual bases by using the reference sequence.

If you want to convert a BAM file to a SAM file, just call:
<pathToExe>/bam --in <bamFile>.bam --out <newSamFile>.sam
Don't forget to put in the paths to the executable and your test files.

= Parameters =
<pre>
Required Parameters:
--in : the SAM/BAM file to be read
--out : the SAM/BAM file to be written
Optional Parameters:
--refFile : reference file name
--noeof : do not expect an EOF block on a bam file.
--params : print the parameter settings
--recover : attempt to recover the input bam file.
Optional Sequence Parameters (only specify one):
--seqOrig : Leave the sequence as is (default & used if reference is not specified).
--seqBases : Convert any '=' in the sequence to the appropriate base using the reference (requires --ref).
--seqEquals : Convert any bases that match the reference to '=' (requires --ref).
</pre>

== Sequence Representation Parameters ==
The sequence parameters options specify how to represent the sequence if the reference is specified (refFile option). If the reference is not specified or seqOrig is specified, no modifications are made to the sequence. If the reference and seqBases is specified, any matches between the sequence and the reference are represented in the sequence as the appropriate base. If the reference and seqEquals is specified, any matches between the sequence and the reference are represented in the sequence as '='.

=== Examples ===
ExtendedCigar: SSMMMDDMMMIMNNNMPMSSS
Sequence: AATAA CTAGA T AGGG
Reference: TAACCCTA ACCCT A
Sequence with Orig: AATAACTAGATAGGG
Sequence with Bases: AATAACTAGATAGGG
Sequence with Equals: AA======G===GGG

ExtendedCigar: SSMMMDDMMMIMNNNMPMSSS
Sequence: AATGA CTGGA T AGGG
Reference: TAACCCTA ACCCT A
Sequence with Orig: AATGACTGGATAGGG
Sequence with Bases: AATGACTGGATAGGG
Sequence with Equals: AA=G===GG===GGG

ExtendedCigar: SSMMMDDMMMIMNNNMPMSSS
Sequence: AAT=A CT=GA T AGGG
Reference: TAACCCTA ACCCT A
Sequence with Orig: AAT=ACT=GATAGGG
Sequence with Bases: AATGACTGGATAGGG
Sequence with Equals: AA======G===GGG

ExtendedCigar: SSMMMDDMMMIMNNNMPMSSS
Sequence: AA=== ===G= = =GGG
Reference: TAACCCTA ACCCT A
Sequence with Orig: AA======G===GGG
Sequence with Bases: AATAACTAGATAGGG
Sequence with Equals: AA======G===GGG

= Usage =
./bam convert --in <inputFile> --out <outputFile.sam/bam/ubam (ubam is uncompressed bam)> [--refFile <reference filename>] [--seqBases|--seqEquals|--seqOrig] [--noeof] [--params]

= Return Value =
Returns the SamStatus for the reads/writes.

= Example Output =
<pre>
Number of records read = 10
Number of records written = 10
</pre>

QPLOT

2010-11-02T20:33:28Z

Pha: initial copy from internal wiki - needs more edits

== Introduction ==
The qplot program is to calculate various summary statistics some of which will be plotted in a pdf file which can be used to assess the sequencing quality for illumina sequencing after mapping reads to the reference genome. The main statistics are empirical Phred scores which was calculated based on the background mismatch rate. By background mismatch rate, it means the rate that sequenced bases are different from the reference genome, EXCLUDING dbSNP positions. Other statistics include GC biases, insert size distribution, depth distribution, genome coverage, empirical Q20 count and so on. An example plot and summary text will follow at the end

== Usage ==
The code is now deposited into the pipeline. If you have not seen it in your pipeline, you shall update you pipeline. See [[pipeline and git]].

In the pipeline folder, go to qplot and make will create an executable qplot in the qplot sub-directory. Here is the qplot help page

wonderland > ./qplot
The following parameters are available. Ones with "[]" are in effect:
References : --reference [/data/local/ref/karma.ref/human.g1k.v37.umfa],
--dbsnp [/home/bingshan/data/db/dbSNP/dbSNP130.UCSC.coordinates.tbl],
--gccontent [/home/bingshan/data/db/gcContent/gcContent.hg37.w250.out]
Create gcContent file : --create_gc [], --winsize [100]
Flag filters : --read1_skip, --read2_skip, --paired_skip,
--unpaired_skip
Dup and QCFail : --dup_keep, --qcfail_keep
Mapping filters : --minMapQuality [0.00]
Records to process : --first_n_record [-1]
Lanes to process : --lanes []
Output files : --plot [], --stats [], --Rcode []
Plot labels : --label [], --bamLabel []

== Input files ==

Three (3) precomputed files are required. Multiple bam/sam files should be appended after all other parameters.

* --reference
** The reference genome is the same as karma reference genome. If the index files do not exist, qplot will create the index files using the input reference fasta file.
** Or use the latest karma to generate reference index files using the command:
karma create -i reference.fa

* --dbsnp
This file has two columns. First column is the chromosome name wich have to be consistent with the reference created above

* --gccontent

Although GC content can be calculated on fly each time, it is much more efficient to load a precomputed GC content from a file. To generate the file, use the following command
/home/bingshan/code/calcGCContent/genGCContent -r referenceByKarma -w windowSize -o gccontent.out

''Note'': Before running the qplot, it is critical to check how the chromosome numbers are coded. Some bam file use just numbers, others use chr + numbers. '''You need to make sure that the chromosome numbers from reference and dbsnp are consistent with the bam file.'''

== Parameters ==

Most of command line parameters are self explanatory and some of them are described here

*Flag filter
** By default all reads are processed. If it is desired to check only the first read of a pair, use --read2_skip to ignore the second read. And so on.

*Duplication and QCFile
** By default reads marked as dup and QCFile are ignored but can be retained by
--dup_keep or --qcfail_keep

*Records to process is to try the first n reads to test the bam files and check whether it works.

* Lanes to process
** If the input bam files have more than one lane and only some of them need to be checked, they can be specified by --lanes 1,3,5 whatever the number of lanes needed.
** In order for this to work, the lane info has to be encoded in the read name such that lane number is the second field with the delimit of ":".

* Plot filters
** GC bias plot requires large memory (genome size which is ~3Gb) and if it is not desired to generate it to save memory it can be disabled by
--gc_plot_skip

* Region list
** If the interest of QA is a list of regions, e.g. exons, this can be achieved by providing a list of region. The regions should be in the form of "chr start end label" each line in the file. In order for this option to work, within each chromosome (contig) the regions have to be sorted by starting position, and also the input bam files have to be sorted.
1 100 500 region_A
1 600 800 region_B
2 100 300 region_C
...

* Plot labels
** Two kinds of labels are enabled. First one is the label for the plot (default is empty), e.g. label on the title of each subplot. Second one is a set of labels for each input bam files, e.g. sample ID (default is numbers 1, 2, ... until the number of input bam files.
--label Run100 --bamLabels s1,s2,s3,s4,s5,s6,s7,s8

* Multiple threading
** Number of concurrent threads running for the input bam files. One bam file will be processed by one thread. Therefore using a number which is dividable by the number of input bam files will make it efficient. One extra thread requires memory about 375Mb on top of around 4Gb memory used to hold reference and GC content file.

== Output files ==

There are three (optional) output files.
* --plot qa.pdf
** This is a pdf file containing 2 pages each with 4 figures. If --pages 1 is specified, only page 1 is output. The plot is generated using Rscript.
* --stats qa.stats
** This is a text file containing various summary statistics for each input bam/sam file
* --Rcode qa.R
** This is the R code used for plotting the figures in qa.pdf file. If Rscript is not installed in the system the qplot is run, you can use the qa.R to generate the figures in other machines, or extract plotting data from each run and combine multiple runs together to generate more comprehensive plots.

== Example output ==

* Figures
https://statgen.sph.umich.edu/w/images/5/53/Sardinia_Run_84_QA.pdf

* Summary statistics text file
TotalReads(e6) 75.74 76.37 76.33 76.47 76.76 76.33 76.34 75.30
MappingRate(%) 95.27 95.21 95.68 95.64 95.49 95.99 95.72 95.54
MapRate_pass(%) 95.27 95.21 95.68 95.64 95.49 95.99 95.72 95.54
ZeroMapQual(%) 3.91 3.90 3.86 3.95 3.65 3.44 3.47 3.46
MapQual<10(%) 4.15 4.16 4.12 4.22 3.91 3.69 3.72 3.72
TargetMapping(%) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
PairedReads(%) 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00
ProperPaired(%) 91.74 90.92 92.53 92.19 91.74 92.64 91.99 92.10
MappedBases(e9) 7.62 7.51 7.56 7.52 7.55 7.60 7.62 7.41
Q20Bases(e9) 7.41 7.35 7.42 7.37 7.40 7.44 7.47 7.28
Q20BasesPct(%) 97.27 97.96 98.16 98.02 98.00 97.99 98.12 98.16
MeanDepth 3.01 2.98 3.01 2.99 2.99 3.00 2.99 2.93
GenomeCover(%) 88.45 87.92 87.74 87.88 88.00 88.37 88.78 88.21
GCBiasMSE 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02
ISize_mode 242 267 284 282 268 258 260 262
ISize_medium 241 268 285 281 269 258 259 262
DupRate(%) 1.58 1.55 1.77 1.71 1.46 1.51 1.73 1.95
QCFailRate(%) 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
BaseComp_A(%) 30.3 30.1 30.6 30.3 30.0 30.4 30.4 30.5
BaseComp_C(%) 19.9 20.1 19.6 19.9 20.2 19.8 19.8 19.7
BaseComp_G(%) 19.9 20.1 19.7 20.0 20.2 19.8 19.9 19.8
BaseComp_T(%) 29.9 29.7 30.1 29.8 29.6 30.0 29.9 30.0
BaseComp_O(%) 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0

SAM

2010-07-30T17:40:32Z

Pha: /* What is SAM */ second test

== What is SAM ==
The '''SAM Format''' is a text format for storing sequence data in a series of tab delimited ASCII columns.

Most often it is generated as a human readable version of its sister [[BAM]] format, which stores the same data in a compressed, indexed, binary form.

Currently, most SAM format data is output from aligners that read FASTQ files and assign the sequences to a position with respect to a known reference genome. In the future, SAM will also be used to archive unaligned sequence data generated directly from sequencing machines.

The current definition of the format is at [[http://samtools.sourceforge.net/SAM1.pdf BAM/SAM Specification]].

If you are writing software to read SAM or BAM data, our C++ [[C++ Library: libbam|BamFile]] is a good resource to use.

This is a test to see how merges work - will be reverted. Or Will It?

This is a second test line to see if merges work.

== What Information is in SAM & BAM ==
SAM files and BAM files contain the same information, but in a different format. Refer to the specs to see a format description.

Both SAM & BAM files contain a header section and an alignment section.
The header section may contain information about the entire file and additional information for alignments. The alignments then associate themselves with specific header information.

The alignment section contains the information for each sequence about where/how it aligns to the reference genome.

=== What Information Does SAM/BAM Have for an Alignment ===
Each Alignment has:
* query name, QNAME (SAM)/read_name (BAM). It is used to group/identify alignments that are together, like paired alignments or a read that appears in multiple alignments.
* a bitwise set of information describing the alignment, FLAG. Provides the following information:
** are there multiple fragments?
** are all fragments properly aligned?
** is this fragment unmapped?
** is the next fragment unmapped?
** is this query the reverse strand?
** is the next fragment the reverse strand?
** is this the 1st fragment?
** is this the last fragment?
** is this a secondary alignment?
** did this read fail quality controls?
** is this read a PCR or optical duplicate?

Not all alignments contain The rest of the alignment fields may be set to default values if the information is unknown.
* reference sequence name, RNAME, often contains the Chromosome name.
* leftmost position of where this alignment maps to the reference, POS. For SAM, the reference starts at 1, so this value is 1-based, while for BAM the reference starts at 0,so this value is 0-based. Beware to always use the correct base when referencing positions.
* mapping quality, MAPQ, which contains the "phred-scaled posterior probability that the mapping position" is wrong. (from SAM-1.pdf)
* string indicating alignment information that allows the storing of clipped, [[SAM#What is a CIGAR?|CIGAR]]
* the reference sequence name of the next alignment in this group, MRNM or RNEXT. In paired alignments, it is the mate's reference sequence name. (A group is alignments with the same query name.)
* leftmost position of where the next alignment in this group maps to the reference, MPOS or PNEXT. For SAM, the reference starts at 1, so this value is 1-based, while for BAM the reference starts at 0,so this value is 0-based. Beware to always use the correct base when referencing positions.
* length of this group from the leftmost position to the rightmost position, ISIZE or TLEN
* the query sequence for this alignment, SEQ
* the query quality for this alignment, [[SAM#What is QUAL?|QUAL]], one for each base in the query sequence.
* Additional optional information is also contained within the alignment, [[SAM#What are TAGs?|TAGs]]. A bunch of different information can be stored here and they appear as key/value pairs. See the spec for a detailed list of commonly used tags and what they mean.

==== What is a CIGAR? ====
You may have heard the term CIGAR, but wondered what it means. Hopefully this section will help clarify it.

The sequence being aligned to a reference may have additional bases that are not in the reference or may be missing bases that are in the reference. The CIGAR string is a sequence of of base lengths and the associated operation. They are used to indicate things like which bases align (either a match/mismatch) with the reference, are deleted from the reference, and are insertions that are not in the reference.

For example:
RefPos: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
Reference: C C A T A C T G A A C T G A C T A A C
Read: ACTAGAATGGCT
Aligning these two:
RefPos: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
Reference: C C A T A C T G A A C T G A C T A A C
Read: A C T A G A A T G G C T
With the alignment above, you get:
POS: 5
CIGAR: 3M1I3M1D5M

The POS indicates that the read aligns starting at position 5 on the reference.
The CIGAR says that the first 3 bases in the read sequence align with the reference. The next base in the read does not exist in the reference. Then 3 bases align with the reference. The next reference base does not exist in the read sequence, then 5 more bases align with the reference. Note that at position 14, the base in the read is different than the reference, but it still counts as an M since it aligns to that position.

==== What is QUAL? ====
QUAL stands for query quality. It is an indicator for how accurate each base in the query sequence (SEQ) is. If QUAL is specified, there is a quality value for each base in SEQ.

Quality is calculated based on the probability that a base is wrong, p, using the following formula:
<math>quality = -10 \log_{10}p</math>
This quality is called the [http://en.wikipedia.org/wiki/Phred_quality_score Phred Quality Score].

Since a human readable format is desired for SAM, 33 is added to the calculated quality in order to make it a printable character ranging from ! - ~.

So, for SAM, the QUAL field is:
<math>QUAL = (-10 \log_{10}p) + 33</math>

Phred Quality is also found in a FASTQ file, described here: http://en.wikipedia.org/wiki/FASTQ_format#Quality

==== What are TAGs? ====
TAGs are optional fields on a SAM/BAM Alignment.
A TAG is comprised of a two character TAG key, they type of the value, and the value:
[A-Za-z][A-za-z]:[AifZH]:.*

There are a set of predefined tags that are general used in Alignments.

A user can also use any additional tags to store any information they want.

== Example SAM ==
=== Example Alignments ===
This is what the alignment section of a SAM file looks like:

1:497:R:-272+13M17D24M 113 1 497 37 37M 15 100338662 0 CGGGTCTGACCTGAGGAGAACTGTGCTCCGCCTTCAG 0;==-==9;>>>>>=>>>>>>>>>>>=>>>>>>>>>> XT:A:U NM:i:0 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:37
19:20389:F:275+18M2D19M 99 1 17644 0 37M = 17919 314 TATGACTGCTAATAATACCTACACATGTTAGAACCAT >>>>>>>>>>>>>>>>>>>><<>>><<>>4::>>:<9 XT:A:R NM:i:0 SM:i:0 AM:i:0 X0:i:4 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:37
19:20389:F:275+18M2D19M 147 1 17919 0 18M2D19M = 17644 -314 GTAGTACCAACTGTAAGTCCTTATCTTCATACTTTGT ;44999;499<8<8<<<8<<><<<<><7<;<<<>><< XT:A:R NM:i:2 SM:i:0 AM:i:0 X0:i:4 X1:i:0 XM:i:0 XO:i:1 XG:i:2 MD:Z:18^CA19
9:21597+10M2I25M:R:-209 83 1 21678 0 8M2I27M = 21469 -244 CACCACATCACATATACCAAGCCTGGCTGTGTCTTCT <;9<<5><<<<><<<>><<><>><9>><>>>9>>><> XT:A:R NM:i:2 SM:i:0 AM:i:0 X0:i:5 X1:i:0 XM:i:0 XO:i:1 XG:i:2 MD:Z:35

In this example, the fields are:
{| style="margin: 1em 1em 1em 0; background-color: #f9f9f9; border: 1px #aaa solid; border-collapse: collapse;" border="1"
|-style="background: #f2f2f2; text-align: center;"
! '''Field''' !! '''Alignment 1''' !! '''Alignment 2''' !! '''Alignment 3''' !! '''Alignment 4'''
|-
|QNAME
|1:497:R:-272+13M17D24M
|19:20389:F:275+18M2D19M
|19:20389:F:275+18M2D19M
|9:21597+10M2I25M:R:-209
|-
|FLAG
|113
|99
|147
|83
|-
|RNAME
|1
|1
|1
|1
|-
|POS
|497
|17644
|17919
|21678
|-
|MAPQ
|37
|0
|0
|0
|-
|CIGAR
|37M
|37M
|18M2D19M
|8M2I27M
|-
|MRNM/RNEXT
|15
|=
|=
|=
|-
|MPOS/PNEXT
|100338662
|17919
|17644
|21469
|-
|ISIZE/TLEN
|0
|314
|-314
|-244
|-
|SEQ
|CGGGTCTGACCTGAGGAGAACTGTGCTCCGCCTTCAG
|TATGACTGCTAATAATACCTACACATGTTAGAACCAT
|GTAGTACCAACTGTAAGTCCTTATCTTCATACTTTGT
|CACCACATCACATATACCAAGCCTGGCTGTGTCTTCT
|-
|QUAL
|0;==-==9;>>>>>=>>>>>>>>>>>=>>>>>>>>>>
|>>>>>>>>>>>>>>>>>>>><<>>><<>>4::>>:<9
|;44999;499<8<8<<<8<<><<<<><7<;<<<>><<
|<;9<<5><<<<><<<>><<><>><9>><>>>9>>><>
|-
|TAGs
|XT:A:U NM:i:0 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:37
|XT:A:R NM:i:0 SM:i:0 AM:i:0 X0:i:4 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:37
|XT:A:R NM:i:2 SM:i:0 AM:i:0 X0:i:4 X1:i:0 XM:i:0 XO:i:1 XG:i:2 MD:Z:18^CA19
|XT:A:R NM:i:2 SM:i:0 AM:i:0 X0:i:5 X1:i:0 XM:i:0 XO:i:1 XG:i:2 MD:Z:35
|}

You should now be a SAM expert ;-)

SAM

2010-07-30T17:38:12Z

Pha: /* What is SAM */ test edit - revert

== What is SAM ==
The '''SAM Format''' is a text format for storing sequence data in a series of tab delimited ASCII columns.

Most often it is generated as a human readable version of its sister [[BAM]] format, which stores the same data in a compressed, indexed, binary form.

Currently, most SAM format data is output from aligners that read FASTQ files and assign the sequences to a position with respect to a known reference genome. In the future, SAM will also be used to archive unaligned sequence data generated directly from sequencing machines.

The current definition of the format is at [[http://samtools.sourceforge.net/SAM1.pdf BAM/SAM Specification]].

If you are writing software to read SAM or BAM data, our C++ [[C++ Library: libbam|BamFile]] is a good resource to use.

This is a test to see how merges work - will be reverted.

== What Information is in SAM & BAM ==
SAM files and BAM files contain the same information, but in a different format. Refer to the specs to see a format description.

Both SAM & BAM files contain a header section and an alignment section.
The header section may contain information about the entire file and additional information for alignments. The alignments then associate themselves with specific header information.

The alignment section contains the information for each sequence about where/how it aligns to the reference genome.

=== What Information Does SAM/BAM Have for an Alignment ===
Each Alignment has:
* query name, QNAME (SAM)/read_name (BAM). It is used to group/identify alignments that are together, like paired alignments or a read that appears in multiple alignments.
* a bitwise set of information describing the alignment, FLAG. Provides the following information:
** are there multiple fragments?
** are all fragments properly aligned?
** is this fragment unmapped?
** is the next fragment unmapped?
** is this query the reverse strand?
** is the next fragment the reverse strand?
** is this the 1st fragment?
** is this the last fragment?
** is this a secondary alignment?
** did this read fail quality controls?
** is this read a PCR or optical duplicate?

Not all alignments contain The rest of the alignment fields may be set to default values if the information is unknown.
* reference sequence name, RNAME, often contains the Chromosome name.
* leftmost position of where this alignment maps to the reference, POS. For SAM, the reference starts at 1, so this value is 1-based, while for BAM the reference starts at 0,so this value is 0-based. Beware to always use the correct base when referencing positions.
* mapping quality, MAPQ, which contains the "phred-scaled posterior probability that the mapping position" is wrong. (from SAM-1.pdf)
* string indicating alignment information that allows the storing of clipped, [[SAM#What is a CIGAR?|CIGAR]]
* the reference sequence name of the next alignment in this group, MRNM or RNEXT. In paired alignments, it is the mate's reference sequence name. (A group is alignments with the same query name.)
* leftmost position of where the next alignment in this group maps to the reference, MPOS or PNEXT. For SAM, the reference starts at 1, so this value is 1-based, while for BAM the reference starts at 0,so this value is 0-based. Beware to always use the correct base when referencing positions.
* length of this group from the leftmost position to the rightmost position, ISIZE or TLEN
* the query sequence for this alignment, SEQ
* the query quality for this alignment, [[SAM#What is QUAL?|QUAL]], one for each base in the query sequence.
* Additional optional information is also contained within the alignment, [[SAM#What are TAGs?|TAGs]]. A bunch of different information can be stored here and they appear as key/value pairs. See the spec for a detailed list of commonly used tags and what they mean.

==== What is a CIGAR? ====
You may have heard the term CIGAR, but wondered what it means. Hopefully this section will help clarify it.

The sequence being aligned to a reference may have additional bases that are not in the reference or may be missing bases that are in the reference. The CIGAR string is a sequence of of base lengths and the associated operation. They are used to indicate things like which bases align (either a match/mismatch) with the reference, are deleted from the reference, and are insertions that are not in the reference.

For example:
RefPos: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
Reference: C C A T A C T G A A C T G A C T A A C
Read: ACTAGAATGGCT
Aligning these two:
RefPos: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
Reference: C C A T A C T G A A C T G A C T A A C
Read: A C T A G A A T G G C T
With the alignment above, you get:
POS: 5
CIGAR: 3M1I3M1D5M

The POS indicates that the read aligns starting at position 5 on the reference.
The CIGAR says that the first 3 bases in the read sequence align with the reference. The next base in the read does not exist in the reference. Then 3 bases align with the reference. The next reference base does not exist in the read sequence, then 5 more bases align with the reference. Note that at position 14, the base in the read is different than the reference, but it still counts as an M since it aligns to that position.

==== What is QUAL? ====
QUAL stands for query quality. It is an indicator for how accurate each base in the query sequence (SEQ) is. If QUAL is specified, there is a quality value for each base in SEQ.

Quality is calculated based on the probability that a base is wrong, p, using the following formula:
<math>quality = -10 \log_{10}p</math>
This quality is called the [http://en.wikipedia.org/wiki/Phred_quality_score Phred Quality Score].

Since a human readable format is desired for SAM, 33 is added to the calculated quality in order to make it a printable character ranging from ! - ~.

So, for SAM, the QUAL field is:
<math>QUAL = (-10 \log_{10}p) + 33</math>

Phred Quality is also found in a FASTQ file, described here: http://en.wikipedia.org/wiki/FASTQ_format#Quality

==== What are TAGs? ====
TAGs are optional fields on a SAM/BAM Alignment.
A TAG is comprised of a two character TAG key, they type of the value, and the value:
[A-Za-z][A-za-z]:[AifZH]:

There are a set of predefined tags that are general used in Alignments.

A user can also use any additional tags to store any information they want.

== Example SAM ==
=== Example Alignments ===
This is what the alignment section of a SAM file looks like:

1:497:R:-272+13M17D24M 113 1 497 37 37M 15 100338662 0 CGGGTCTGACCTGAGGAGAACTGTGCTCCGCCTTCAG 0;==-==9;>>>>>=>>>>>>>>>>>=>>>>>>>>>> XT:A:U NM:i:0 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:37
19:20389:F:275+18M2D19M 99 1 17644 0 37M = 17919 314 TATGACTGCTAATAATACCTACACATGTTAGAACCAT >>>>>>>>>>>>>>>>>>>><<>>><<>>4::>>:<9 XT:A:R NM:i:0 SM:i:0 AM:i:0 X0:i:4 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:37
19:20389:F:275+18M2D19M 147 1 17919 0 18M2D19M = 17644 -314 GTAGTACCAACTGTAAGTCCTTATCTTCATACTTTGT ;44999;499<8<8<<<8<<><<<<><7<;<<<>><< XT:A:R NM:i:2 SM:i:0 AM:i:0 X0:i:4 X1:i:0 XM:i:0 XO:i:1 XG:i:2 MD:Z:18^CA19
9:21597+10M2I25M:R:-209 83 1 21678 0 8M2I27M = 21469 -244 CACCACATCACATATACCAAGCCTGGCTGTGTCTTCT <;9<<5><<<<><<<>><<><>><9>><>>>9>>><> XT:A:R NM:i:2 SM:i:0 AM:i:0 X0:i:5 X1:i:0 XM:i:0 XO:i:1 XG:i:2 MD:Z:35

In this example, the fields are:
{| style="margin: 1em 1em 1em 0; background-color: #f9f9f9; border: 1px #aaa solid; border-collapse: collapse;" border="1"
|-style="background: #f2f2f2; text-align: center;"
! '''Field''' !! '''Alignment 1''' !! '''Alignment 2''' !! '''Alignment 3''' !! '''Alignment 4'''
|-
|QNAME
|1:497:R:-272+13M17D24M
|19:20389:F:275+18M2D19M
|19:20389:F:275+18M2D19M
|9:21597+10M2I25M:R:-209
|-
|FLAG
|113
|99
|147
|83
|-
|RNAME
|1
|1
|1
|1
|-
|POS
|497
|17644
|17919
|21678
|-
|MAPQ
|37
|0
|0
|0
|-
|CIGAR
|37M
|37M
|18M2D19M
|8M2I27M
|-
|MRNM/RNEXT
|15
|=
|=
|=
|-
|MPOS/PNEXT
|100338662
|17919
|17644
|21469
|-
|ISIZE/TLEN
|0
|314
|-314
|-244
|-
|SEQ
|CGGGTCTGACCTGAGGAGAACTGTGCTCCGCCTTCAG
|TATGACTGCTAATAATACCTACACATGTTAGAACCAT
|GTAGTACCAACTGTAAGTCCTTATCTTCATACTTTGT
|CACCACATCACATATACCAAGCCTGGCTGTGTCTTCT
|-
|QUAL
|0;==-==9;>>>>>=>>>>>>>>>>>=>>>>>>>>>>
|>>>>>>>>>>>>>>>>>>>><<>>><<>>4::>>:<9
|;44999;499<8<8<<<8<<><<<<><7<;<<<>><<
|<;9<<5><<<<><<<>><<><>><9>><>>>9>>><>
|-
|TAGs
|XT:A:U NM:i:0 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:37
|XT:A:R NM:i:0 SM:i:0 AM:i:0 X0:i:4 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:37
|XT:A:R NM:i:2 SM:i:0 AM:i:0 X0:i:4 X1:i:0 XM:i:0 XO:i:1 XG:i:2 MD:Z:18^CA19
|XT:A:R NM:i:2 SM:i:0 AM:i:0 X0:i:5 X1:i:0 XM:i:0 XO:i:1 XG:i:2 MD:Z:35
|}

SAM

2010-07-30T17:30:02Z

Pha: /* What is SAM */ refer to BamFile c++ library code

== What is SAM ==
The '''SAM Format''' is a text format for storing sequence data in a series of tab delimited ASCII columns.

Most often it is generated as a human readable version of its sister [[BAM]] format, which stores the same data in a compressed, indexed, binary form.

Currently, most SAM format data is output from aligners that read FASTQ files and assign the sequences to a position with respect to a known reference genome. In the future, SAM will also be used to archive unaligned sequence data generated directly from sequencing machines.

The current definition of the format is at [[http://samtools.sourceforge.net/SAM1.pdf BAM/SAM Specification]].

If you are writing software to read SAM or BAM data, our C++ [[C++ Library: libbam|BamFile]] is a good resource to use.

== What Information is in SAM & BAM ==
SAM files and BAM files contain the same information, but in a different format. Refer to the specs to see a format description.

Both SAM & BAM files contain a header section and an alignment section.
The header section may contain information about the entire file and additional information for alignments. The alignments then associate themselves with specific header information.

The alignment section contains the information for each sequence about where/how it aligns to the reference genome.

=== What Information Does SAM/BAM Have for an Alignment ===
Each Alignment has:
* query name, QNAME (SAM)/read_name (BAM). It is used to group/identify alignments that are together, like paired alignments or a read that appears in multiple alignments.
* a bitwise set of information describing the alignment, FLAG. Provides the following information:
** are there multiple fragments?
** are all fragments properly aligned?
** is this fragment unmapped?
** is the next fragment unmapped?
** is this query the reverse strand?
** is the next fragment the reverse strand?
** is this the 1st fragment?
** is this the last fragment?
** is this a secondary alignment?
** did this read fail quality controls?
** is this read a PCR or optical duplicate?

Not all alignments contain The rest of the alignment fields may be set to default values if the information is unknown.
* reference sequence name, RNAME, often contains the Chromosome name.
* leftmost position of where this alignment maps to the reference, POS. For SAM, the reference starts at 1, so this value is 1-based, while for BAM the reference starts at 0,so this value is 0-based. Beware to always use the correct base when referencing positions.
* mapping quality, MAPQ, which contains the "phred-scaled posterior probability that the mapping position" is wrong. (from SAM-1.pdf)
* string indicating alignment information that allows the storing of clipped, [[SAM#What is a CIGAR?|CIGAR]]
* the reference sequence name of the next alignment in this group, MRNM or RNEXT. In paired alignments, it is the mate's reference sequence name. (A group is alignments with the same query name.)
* leftmost position of where the next alignment in this group maps to the reference, MPOS or PNEXT. For SAM, the reference starts at 1, so this value is 1-based, while for BAM the reference starts at 0,so this value is 0-based. Beware to always use the correct base when referencing positions.
* length of this group from the leftmost position to the rightmost position, ISIZE or TLEN
* the query sequence for this alignment, SEQ
* the query quality for this alignment, [[SAM#What is QUAL?|QUAL]], one for each base in the query sequence.
* Additional optional information is also contained within the alignment, [[SAM#What are TAGs?|TAGs]]. A bunch of different information can be stored here and they appear as key/value pairs. See the spec for a detailed list of commonly used tags and what they mean.

==== What is a CIGAR? ====
You may have heard the term CIGAR, but wondered what it means. Hopefully this section will help clarify it.

The sequence being aligned to a reference may have additional bases that are not in the reference or may be missing bases that are in the reference. The CIGAR string is a sequence of of base lengths and the associated operation. They are used to indicate things like which bases align (either a match/mismatch) with the reference, are deleted from the reference, and are insertions that are not in the reference.

For example:
RefPos: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
Reference: C C A T A C T G A A C T G A C T A A C
Read: ACTAGAATGGCT
Aligning these two:
RefPos: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
Reference: C C A T A C T G A A C T G A C T A A C
Read: A C T A G A A T G G C T
With the alignment above, you get:
POS: 5
CIGAR: 3M1I3M1D5M

The POS indicates that the read aligns starting at position 5 on the reference.
The CIGAR says that the first 3 bases in the read sequence align with the reference. The next base in the read does not exist in the reference. Then 3 bases align with the reference. The next reference base does not exist in the read sequence, then 5 more bases align with the reference. Note that at position 14, the base in the read is different than the reference, but it still counts as an M since it aligns to that position.

==== What is QUAL? ====
QUAL stands for query quality. It is an indicator for how accurate each base in the query sequence (SEQ) is. If QUAL is specified, there is a quality value for each base in SEQ.

Quality is calculated based on the probability that a base is wrong, p, using the following formula:
<math>quality = -10 \log_{10}p</math>
This quality is called the [http://en.wikipedia.org/wiki/Phred_quality_score Phred Quality Score].

Since a human readable format is desired for SAM, 33 is added to the calculated quality in order to make it a printable character ranging from ! - ~.

So, for SAM, the QUAL field is:
<math>QUAL = (-10 \log_{10}p) + 33</math>

Phred Quality is also found in a FASTQ file, described here: http://en.wikipedia.org/wiki/FASTQ_format#Quality

==== What are TAGs? ====

== Example SAM ==
=== Example Alignments ===
This is what the alignment section of a SAM file looks like:

1:497:R:-272+13M17D24M 113 1 497 37 37M 15 100338662 0 CGGGTCTGACCTGAGGAGAACTGTGCTCCGCCTTCAG 0;==-==9;>>>>>=>>>>>>>>>>>=>>>>>>>>>> XT:A:U NM:i:0 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:37
19:20389:F:275+18M2D19M 99 1 17644 0 37M = 17919 314 TATGACTGCTAATAATACCTACACATGTTAGAACCAT >>>>>>>>>>>>>>>>>>>><<>>><<>>4::>>:<9 XT:A:R NM:i:0 SM:i:0 AM:i:0 X0:i:4 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:37
19:20389:F:275+18M2D19M 147 1 17919 0 18M2D19M = 17644 -314 GTAGTACCAACTGTAAGTCCTTATCTTCATACTTTGT ;44999;499<8<8<<<8<<><<<<><7<;<<<>><< XT:A:R NM:i:2 SM:i:0 AM:i:0 X0:i:4 X1:i:0 XM:i:0 XO:i:1 XG:i:2 MD:Z:18^CA19
9:21597+10M2I25M:R:-209 83 1 21678 0 8M2I27M = 21469 -244 CACCACATCACATATACCAAGCCTGGCTGTGTCTTCT <;9<<5><<<<><<<>><<><>><9>><>>>9>>><> XT:A:R NM:i:2 SM:i:0 AM:i:0 X0:i:5 X1:i:0 XM:i:0 XO:i:1 XG:i:2 MD:Z:35

In this example, the fields are:
{| style="margin: 1em 1em 1em 0; background-color: #f9f9f9; border: 1px #aaa solid; border-collapse: collapse;" border="1"
|-style="background: #f2f2f2; text-align: center;"
! '''Field''' !! '''Alignment 1''' !! '''Alignment 2''' !! '''Alignment 3''' !! '''Alignment 4'''
|-
|QNAME
|1:497:R:-272+13M17D24M
|19:20389:F:275+18M2D19M
|19:20389:F:275+18M2D19M
|9:21597+10M2I25M:R:-209
|-
|FLAG
|113
|99
|147
|83
|-
|RNAME
|1
|1
|1
|1
|-
|POS
|497
|17644
|17919
|21678
|-
|MAPQ
|37
|0
|0
|0
|-
|CIGAR
|37M
|37M
|18M2D19M
|8M2I27M
|-
|MRNM/RNEXT
|15
|=
|=
|=
|-
|MPOS/PNEXT
|100338662
|17919
|17644
|21469
|-
|ISIZE/TLEN
|0
|314
|-314
|-244
|-
|SEQ
|CGGGTCTGACCTGAGGAGAACTGTGCTCCGCCTTCAG
|TATGACTGCTAATAATACCTACACATGTTAGAACCAT
|GTAGTACCAACTGTAAGTCCTTATCTTCATACTTTGT
|CACCACATCACATATACCAAGCCTGGCTGTGTCTTCT
|-
|QUAL
|0;==-==9;>>>>>=>>>>>>>>>>>=>>>>>>>>>>
|>>>>>>>>>>>>>>>>>>>><<>>><<>>4::>>:<9
|;44999;499<8<8<<<8<<><<<<><7<;<<<>><<
|<;9<<5><<<<><<<>><<><>><9>><>>>9>>><>
|-
|TAGs
|XT:A:U NM:i:0 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:37
|XT:A:R NM:i:0 SM:i:0 AM:i:0 X0:i:4 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:37
|XT:A:R NM:i:2 SM:i:0 AM:i:0 X0:i:4 X1:i:0 XM:i:0 XO:i:1 XG:i:2 MD:Z:18^CA19
|XT:A:R NM:i:2 SM:i:0 AM:i:0 X0:i:5 X1:i:0 XM:i:0 XO:i:1 XG:i:2 MD:Z:35
|}

BAM

2010-07-30T17:28:27Z

Pha: elaborate on details

The '''BAM Format''' is a binary format for storing sequence data.

The current definition of the format is at [[http://samtools.sourceforge.net/SAM1.pdf BAM/SAM Specification]].

The corresponding [[SAM Format]] can be used to store sequence data, both aligned as well as unaligned, in a human readable format.

BAM and SAM formats are designed to contain exactly the same information. The SAM format is more human readable, and easier to process by conventional text based processing programs, such as awk, sed, python, cut and so on. The BAM format provides binary versions of most of the same data, and is designed to compress reasonably well.

Our [[C++ Library: libbam|BamFile]] library reads both SAM and BAM format files.

SAM

2010-07-29T20:30:44Z

Pha: /* What is SAM */ simplify and elaborate

== What is SAM ==
The '''SAM Format''' is a text format for storing sequence data in a series of tab delimited ASCII columns.

Most often it is generated as a human readable version of its sister [[BAM]] format, which stores the same data in a compressed, indexed, binary form.

Currently, most SAM format data is output from aligners that read FASTQ files and assign the sequences to a position with respect to a genome. In the future, SAM will also be used to archive unaligned sequence data generated directly from sequencing machines.

The current definition of the format is at [[http://samtools.sourceforge.net/SAM1.pdf BAM/SAM Specification]].

== What Information is in SAM & BAM ==
SAM files and BAM files contain the same information, but in a different format. Refer to the specs to see a format description.

Both SAM & BAM files contain a header section and an alignment section.
The header section may contain information about the entire file and additional information for alignments. The alignments then associate themselves with specific header information.

The alignment section contains the information for each sequence about where/how it aligns to the reference genome.

=== What Information Does SAM/BAM Have for an Alignment ===
Each Alignment has:
* query name, QNAME (SAM)/read_name (BAM). It is used to group/identify alignments that are together, like paired alignments or a read that appears in multiple alignments.
* a bitwise set of information describing the alignment, FLAG. Provides the following information:
** are there multiple fragments?
** are all fragments properly aligned?
** is this fragment unmapped?
** is the next fragment unmapped?
** is this query the reverse strand?
** is the next fragment the reverse strand?
** is this the 1st fragment?
** is this the last fragment?
** is this a secondary alignment?
** did this read fail quality controls?
** is this read a PCR or optical duplicate?

Not all alignments contain The rest of the alignment fields may be set to default values if the information is unknown.
* reference sequence name, RNAME, often contains the Chromosome name.
* leftmost position of where this alignment maps to the reference, POS. For SAM, the reference starts at 1, so this value is 1-based, while for BAM the reference starts at 0,so this value is 0-based. Beware to always use the correct base when referencing positions.
* mapping quality, MAPQ, which contains the "phred-scaled posterior probability that the mapping position" is wrong. (from SAM-1.pdf)
* string indicating alignment information that allows the storing of clipped, [[SAM#What is a CIGAR?|CIGAR]]
* the reference sequence name of the next alignment in this group, MRNM or RNEXT. In paired alignments, it is the mate's reference sequence name. (A group is alignments with the same query name.)
* leftmost position of where the next alignment in this group maps to the reference, MPOS or PNEXT. For SAM, the reference starts at 1, so this value is 1-based, while for BAM the reference starts at 0,so this value is 0-based. Beware to always use the correct base when referencing positions.
* length of this group from the leftmost position to the rightmost position, ISIZE or TLEN
* the query sequence for this alignment, SEQ
* the query quality for this alignment, QUAL, one for each base in the query sequence.
* Additional optional information is also contained within the alignment, TAGS. A bunch of different information can be stored here and they appear as key/value pairs. See the spec for a detailed list of commonly used tags and what they mean.

==== What is a CIGAR? ====
You may have heard the term CIGAR, but wondered what it means. Hopefully this section will help clarify it.

The sequence being aligned to a reference may have additional bases that are not in the reference or may be missing bases that are in the reference. The CIGAR string is a sequence of of base lengths and the associated operation. They are used to indicate things like which bases align (either a match/mismatch) with the reference, are deleted from the reference, and are insertions that are not in the reference.

For example:
RefPos: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
Reference: C C A T A C T G A A C T G A C T A A C
Read: ACTAGAATGGCT
Aligning these two:
RefPos: 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
Reference: C C A T A C T G A A C T G A C T A A C
Read: A C T A G A A T G G C T
With the alignment above, you get:
POS: 5
CIGAR: 3M1I3M1D5M

The POS indicates that the read aligns starting at position 5 on the reference.
The CIGAR says that the first 3 bases in the read sequence align with the reference. The next base in the read does not exist in the reference. Then 3 bases align with the reference. The next reference base does not exist in the read sequence, then 5 more bases align with the reference. Note that at position 14, the base in the read is different than the reference, but it still counts as an M since it aligns to that position.

== Example SAM ==
=== Example Alignments ===
This is what the alignment section of a SAM file looks like:

1:497:R:-272+13M17D24M 113 1 497 37 37M 15 100338662 0 CGGGTCTGACCTGAGGAGAACTGTGCTCCGCCTTCAG 0;==-==9;>>>>>=>>>>>>>>>>>=>>>>>>>>>> XT:A:U NM:i:0 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:37
19:20389:F:275+18M2D19M 99 1 17644 0 37M = 17919 314 TATGACTGCTAATAATACCTACACATGTTAGAACCAT >>>>>>>>>>>>>>>>>>>><<>>><<>>4::>>:<9 XT:A:R NM:i:0 SM:i:0 AM:i:0 X0:i:4 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:37
19:20389:F:275+18M2D19M 147 1 17919 0 18M2D19M = 17644 -314 GTAGTACCAACTGTAAGTCCTTATCTTCATACTTTGT ;44999;499<8<8<<<8<<><<<<><7<;<<<>><< XT:A:R NM:i:2 SM:i:0 AM:i:0 X0:i:4 X1:i:0 XM:i:0 XO:i:1 XG:i:2 MD:Z:18^CA19
9:21597+10M2I25M:R:-209 83 1 21678 0 8M2I27M = 21469 -244 CACCACATCACATATACCAAGCCTGGCTGTGTCTTCT <;9<<5><<<<><<<>><<><>><9>><>>>9>>><> XT:A:R NM:i:2 SM:i:0 AM:i:0 X0:i:5 X1:i:0 XM:i:0 XO:i:1 XG:i:2 MD:Z:35

In this example, the fields are:
{| style="margin: 1em 1em 1em 0; background-color: #f9f9f9; border: 1px #aaa solid; border-collapse: collapse;" border="1"
|-style="background: #f2f2f2; text-align: center;"
! '''Field''' !! '''Alignment 1''' !! '''Alignment 2''' !! '''Alignment 3''' !! '''Alignment 4'''
|-
|QNAME
|1:497:R:-272+13M17D24M
|19:20389:F:275+18M2D19M
|19:20389:F:275+18M2D19M
|9:21597+10M2I25M:R:-209
|-
|FLAG
|113
|99
|147
|83
|-
|RNAME
|1
|1
|1
|1
|-
|POS
|497
|17644
|17919
|21678
|-
|MAPQ
|37
|0
|0
|0
|-
|CIGAR
|37M
|37M
|18M2D19M
|8M2I27M
|-
|MRNM/RNEXT
|15
|=
|=
|=
|-
|MPOS/PNEXT
|100338662
|17919
|17644
|21469
|-
|ISIZE/TLEN
|0
|314
|-314
|-244
|-
|SEQ
|CGGGTCTGACCTGAGGAGAACTGTGCTCCGCCTTCAG
|TATGACTGCTAATAATACCTACACATGTTAGAACCAT
|GTAGTACCAACTGTAAGTCCTTATCTTCATACTTTGT
|CACCACATCACATATACCAAGCCTGGCTGTGTCTTCT
|-
|QUAL
|0;==-==9;>>>>>=>>>>>>>>>>>=>>>>>>>>>>
|>>>>>>>>>>>>>>>>>>>><<>>><<>>4::>>:<9
|;44999;499<8<8<<<8<<><<<<><7<;<<<>><<
|<;9<<5><<<<><<<>><<><>><9>><>>>9>>><>
|-
|TAGs
|XT:A:U NM:i:0 SM:i:37 AM:i:0 X0:i:1 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:37
|XT:A:R NM:i:0 SM:i:0 AM:i:0 X0:i:4 X1:i:0 XM:i:0 XO:i:0 XG:i:0 MD:Z:37
|XT:A:R NM:i:2 SM:i:0 AM:i:0 X0:i:4 X1:i:0 XM:i:0 XO:i:1 XG:i:2 MD:Z:18^CA19
|XT:A:R NM:i:2 SM:i:0 AM:i:0 X0:i:5 X1:i:0 XM:i:0 XO:i:1 XG:i:2 MD:Z:35
|}

BamGenotypeCheck

2010-07-16T17:06:33Z

Pha: fix indenting

'''bamGenotypeCheck''' is a program that verifies whether the reads in particular file match previously known genotypes for an individual (or group of individuals).

== Download bamGenotypeCheck ==

To get a copy go to the [http://www.sph.umich.edu/csg/pha/karma/download/ Karma Download] download page.

== Build bamGenotypeCheck ==

Karma (which includes bamGenotypeCheck) is designed to be reasonably portable.

However, since development occurs only on Ubuntu 9.10 x86 and x64 platforms, and later, there are likely other portability issues.

We support Karma only on Ubuntu 9.10 and later on 64-bit processors.

== Usage ==

A key step in any genetic analysis is to verify whether data being generated matches expectations. This program checks whether reads in a BAM file match previous genotypes for a specific sample.

Using a mathematical model that relates observed sequence reads to an hypothetical true genotype, bamGenotypeCheck tries to decide whether sequence reads match a particular individual or are more likely to be contaminated (including a small proportion of foreign DNA), derived from a closely related individual, or derived from a completely different individual.

== Basic Usage Example ==

Here is a typical command line:

bamGenotypeCheck -r /data/local/ref/karma.ref/human.g1k.v37.fa \
-k BAMfiles.txt -p test.ped -d test.dat -m test.map

== Command Line Options ==

=== Input Files ===

-r ''genome reference in [http://en.wikipedia.org/wiki/Fasta_format simplified FASTA format]''
-a ''allele Frequency file in [[MERLIN format]]''
-p ''pedigree file in [[MERLIN format]]''
-d ''data file in [[MERLIN format]]''
-m ''map file in [[MERLIN format]]''

-k ''a list of BAM files to check''
-c [int] ''stop after reading [int] filtered sequence reads''
-C [int] ''stop after reading [int] reads, filtered or not''

=== Output Options ===

-v ''verbose output''

=== Filtering ===

-b [int] ''exclude bases with quality less than [int]''
-M [int] ''exclude reads with map quality less than [int]''
-f [float] ''drop markers with minor allele frequency smaller than [float]''
-F [int] ''set custom BAM flags filter (not implemented at the moment)''

=== Other Options ===

-e [float] '' set minimum error base error to [float]''

== Principle of Operation ==

Each read group in a BAM file is evaluated independently. This means that in file with multiple read groups, problems will be flagged at the read group level (a plus). However, it also means that it might be hard to discern the correct assignment of read groups with very little data.

For each aligned base that overlaps a known genotype, we calculate the probability the probability that it was derived from a particular known genotype. This comparison considers only bases that overlap previously known genotypes and that meet the base quality and mapping quality thresholds.

Each individual in a pedigree has a different combination of genotypes, and bamGenotypeCheck will systematically search for the individual whose genotypes best match the observed read data.

For more about the technical details, see the page [[Verifying Sample Identities - Implementation]]

== TODO ==

BamGenotypeCheck

2010-07-16T17:05:10Z

Pha: add pointer to fresh karma download

'''bamGenotypeCheck''' is a program that verifies whether the reads in particular file match previously known genotypes for an individual (or group of individuals).

= Download bamGenotypeCheck =

To get a copy go to the [http://www.sph.umich.edu/csg/pha/karma/download/ Karma Download] download page.

= Build bamGenotypeCheck =

Karma (which includes bamGenotypeCheck) is designed to be reasonably portable.

However, since development occurs only on Ubuntu 9.10 x86 and x64 platforms, and later, there are likely other portability issues.

We support Karma only on Ubuntu 9.10 and later on 64-bit processors.

== Usage ==

A key step in any genetic analysis is to verify whether data being generated matches expectations. This program checks whether reads in a BAM file match previous genotypes for a specific sample.

Using a mathematical model that relates observed sequence reads to an hypothetical true genotype, bamGenotypeCheck tries to decide whether sequence reads match a particular individual or are more likely to be contaminated (including a small proportion of foreign DNA), derived from a closely related individual, or derived from a completely different individual.

== Basic Usage Example ==

Here is a typical command line:

bamGenotypeCheck -r /data/local/ref/karma.ref/human.g1k.v37.fa \
-k BAMfiles.txt -p test.ped -d test.dat -m test.map

== Command Line Options ==

=== Input Files ===

-r ''genome reference in [http://en.wikipedia.org/wiki/Fasta_format simplified FASTA format]''
-a ''allele Frequency file in [[MERLIN format]]''
-p ''pedigree file in [[MERLIN format]]''
-d ''data file in [[MERLIN format]]''
-m ''map file in [[MERLIN format]]''

-k ''a list of BAM files to check''
-c [int] ''stop after reading [int] filtered sequence reads''
-C [int] ''stop after reading [int] reads, filtered or not''

=== Output Options ===

-v ''verbose output''

=== Filtering ===

-b [int] ''exclude bases with quality less than [int]''
-M [int] ''exclude reads with map quality less than [int]''
-f [float] ''drop markers with minor allele frequency smaller than [float]''
-F [int] ''set custom BAM flags filter (not implemented at the moment)''

=== Other Options ===

-e [float] '' set minimum error base error to [float]''

== Principle of Operation ==

Each read group in a BAM file is evaluated independently. This means that in file with multiple read groups, problems will be flagged at the read group level (a plus). However, it also means that it might be hard to discern the correct assignment of read groups with very little data.

For each aligned base that overlaps a known genotype, we calculate the probability the probability that it was derived from a particular known genotype. This comparison considers only bases that overlap previously known genotypes and that meet the base quality and mapping quality thresholds.

Each individual in a pedigree has a different combination of genotypes, and bamGenotypeCheck will systematically search for the individual whose genotypes best match the observed read data.

For more about the technical details, see the page [[Verifying Sample Identities - Implementation]]

== TODO ==

BamGenotypeCheck

2010-06-22T14:57:18Z

Pha: /* Input Files */ fix link

'''bamGenotypeCheck''' is a program that verifies whether the reads in particular file match previously known genotypes for an individual (or group of individuals).

== Usage ==

A key step in any genetic analysis is to verify whether data being generated matches expectations. This program checks whether reads in a BAM file match previous genotypes for a specific sample.

Using a mathematical model that relates observed sequence reads to an hypothetical true genotype, bamGenotypeCheck tries to decide whether sequence reads match a particular individual or are more likely to be contaminated (including a small proportion of foreign DNA), derived from a closely related individual, or derived from a completely different individual.

== Basic Usage Example ==

Here is a typical command line:

bamGenotypeCheck -r /data/local/ref/karma.ref/human.g1k.v37.fa \
-k BAMfiles.txt -p test.ped -d test.dat -m test.map

== Command Line Options ==

=== Input Files ===

-r ''genome reference in [http://en.wikipedia.org/wiki/Fasta_format simplified FASTA format]''
-a ''allele Frequency file in [[MERLIN format]]''
-p ''pedigree file in [[MERLIN format]]''
-d ''data file in [[MERLIN format]]''
-m ''map file in [[MERLIN format]]''

-k ''a list of BAM files to check''
-c [int] ''stop after reading [int] filtered sequence reads''
-C [int] ''stop after reading [int] reads, filtered or not''

=== Output Options ===

-v ''verbose output''

=== Filtering ===

-b [int] ''exclude bases with quality less than [int]''
-M [int] ''exclude reads with map quality less than [int]''
-f [float] ''drop markers with minor allele frequency smaller than [float]''
-F [int] ''set custom BAM flags filter (not implemented at the moment)''

=== Other Options ===

-e [float] '' set minimum error base error to [float]''

== Principle of Operation ==

Each read group in a BAM file is evaluated independently. This means that in file with multiple read groups, problems will be flagged at the read group level (a plus). However, it also means that it might be hard to discern the correct assignment of read groups with very little data.

For each aligned base that overlaps a known genotype, we calculate the probability the probability that it was derived from a particular known genotype. This comparison considers only bases that overlap previously known genotypes and that meet the base quality and mapping quality thresholds.

Each individual in a pedigree has a different combination of genotypes, and bamGenotypeCheck will systematically search for the individual whose genotypes best match the observed read data.

For more about the technical details, see the page [[Verifying Sample Identities - Implementation]]

== TODO ==

BamGenotypeCheck

2010-06-22T14:56:05Z

Pha: /* Input Files */ add fasta link

'''bamGenotypeCheck''' is a program that verifies whether the reads in particular file match previously known genotypes for an individual (or group of individuals).

== Usage ==

A key step in any genetic analysis is to verify whether data being generated matches expectations. This program checks whether reads in a BAM file match previous genotypes for a specific sample.

Using a mathematical model that relates observed sequence reads to an hypothetical true genotype, bamGenotypeCheck tries to decide whether sequence reads match a particular individual or are more likely to be contaminated (including a small proportion of foreign DNA), derived from a closely related individual, or derived from a completely different individual.

== Basic Usage Example ==

Here is a typical command line:

bamGenotypeCheck -r /data/local/ref/karma.ref/human.g1k.v37.fa \
-k BAMfiles.txt -p test.ped -d test.dat -m test.map

== Command Line Options ==

=== Input Files ===

-r ''genome reference in [[http://en.wikipedia.org/wiki/Fasta_format|simplified FASTA format]]''
-a ''allele Frequency file in [[MERLIN format]]''
-p ''pedigree file in [[MERLIN format]]''
-d ''data file in [[MERLIN format]]''
-m ''map file in [[MERLIN format]]''

-k ''a list of BAM files to check''
-c [int] ''stop after reading [int] filtered sequence reads''
-C [int] ''stop after reading [int] reads, filtered or not''

=== Output Options ===

-v ''verbose output''

=== Filtering ===

-b [int] ''exclude bases with quality less than [int]''
-M [int] ''exclude reads with map quality less than [int]''
-f [float] ''drop markers with minor allele frequency smaller than [float]''
-F [int] ''set custom BAM flags filter (not implemented at the moment)''

=== Other Options ===

-e [float] '' set minimum error base error to [float]''

== Principle of Operation ==

Each read group in a BAM file is evaluated independently. This means that in file with multiple read groups, problems will be flagged at the read group level (a plus). However, it also means that it might be hard to discern the correct assignment of read groups with very little data.

For each aligned base that overlaps a known genotype, we calculate the probability the probability that it was derived from a particular known genotype. This comparison considers only bases that overlap previously known genotypes and that meet the base quality and mapping quality thresholds.

Each individual in a pedigree has a different combination of genotypes, and bamGenotypeCheck will systematically search for the individual whose genotypes best match the observed read data.

For more about the technical details, see the page [[Verifying Sample Identities - Implementation]]

== TODO ==

BamGenotypeCheck

2010-06-22T14:54:21Z

Pha: /* Input Files */ more detail

'''bamGenotypeCheck''' is a program that verifies whether the reads in particular file match previously known genotypes for an individual (or group of individuals).

== Usage ==

A key step in any genetic analysis is to verify whether data being generated matches expectations. This program checks whether reads in a BAM file match previous genotypes for a specific sample.

Using a mathematical model that relates observed sequence reads to an hypothetical true genotype, bamGenotypeCheck tries to decide whether sequence reads match a particular individual or are more likely to be contaminated (including a small proportion of foreign DNA), derived from a closely related individual, or derived from a completely different individual.

== Basic Usage Example ==

Here is a typical command line:

bamGenotypeCheck -r /data/local/ref/karma.ref/human.g1k.v37.fa \
-k BAMfiles.txt -p test.ped -d test.dat -m test.map

== Command Line Options ==

=== Input Files ===

-r ''FASTA format genome reference''
-a ''allele Frequency file in [[MERLIN format]]''
-p ''pedigree file in [[MERLIN format]]''
-d ''data file in [[MERLIN format]]''
-m ''map file in [[MERLIN format]]''

-k ''a list of BAM files to check''
-c [int] ''stop after reading [int] filtered sequence reads''
-C [int] ''stop after reading [int] reads, filtered or not''

=== Output Options ===

-v ''verbose output''

=== Filtering ===

-b [int] ''exclude bases with quality less than [int]''
-M [int] ''exclude reads with map quality less than [int]''
-f [float] ''drop markers with minor allele frequency smaller than [float]''
-F [int] ''set custom BAM flags filter (not implemented at the moment)''

=== Other Options ===

-e [float] '' set minimum error base error to [float]''

== Principle of Operation ==

Each read group in a BAM file is evaluated independently. This means that in file with multiple read groups, problems will be flagged at the read group level (a plus). However, it also means that it might be hard to discern the correct assignment of read groups with very little data.

For each aligned base that overlaps a known genotype, we calculate the probability the probability that it was derived from a particular known genotype. This comparison considers only bases that overlap previously known genotypes and that meet the base quality and mapping quality thresholds.

Each individual in a pedigree has a different combination of genotypes, and bamGenotypeCheck will systematically search for the individual whose genotypes best match the observed read data.

For more about the technical details, see the page [[Verifying Sample Identities - Implementation]]

== TODO ==

BamGenotypeCheck

2010-06-22T14:53:20Z

Pha: /* Filtering */ add -f option

'''bamGenotypeCheck''' is a program that verifies whether the reads in particular file match previously known genotypes for an individual (or group of individuals).

== Usage ==

A key step in any genetic analysis is to verify whether data being generated matches expectations. This program checks whether reads in a BAM file match previous genotypes for a specific sample.

Using a mathematical model that relates observed sequence reads to an hypothetical true genotype, bamGenotypeCheck tries to decide whether sequence reads match a particular individual or are more likely to be contaminated (including a small proportion of foreign DNA), derived from a closely related individual, or derived from a completely different individual.

== Basic Usage Example ==

Here is a typical command line:

bamGenotypeCheck -r /data/local/ref/karma.ref/human.g1k.v37.fa \
-k BAMfiles.txt -p test.ped -d test.dat -m test.map

== Command Line Options ==

=== Input Files ===

-r ''FASTA format genome reference''
-a ''allele Frequency file''
-p ''pedigree file in [[MERLIN format]]''
-d ''data file in [[MERLIN format]]''
-m ''map file in [[MERLIN format]]''

-k ''a list of BAM files to check''
-c [int] ''stop after reading [int] filtered sequence reads''
-C [int] ''stop after reading [int] reads, filtered or not''

=== Output Options ===

-v ''verbose output''

=== Filtering ===

-b [int] ''exclude bases with quality less than [int]''
-M [int] ''exclude reads with map quality less than [int]''
-f [float] ''drop markers with minor allele frequency smaller than [float]''
-F [int] ''set custom BAM flags filter (not implemented at the moment)''

=== Other Options ===

-e [float] '' set minimum error base error to [float]''

== Principle of Operation ==

Each read group in a BAM file is evaluated independently. This means that in file with multiple read groups, problems will be flagged at the read group level (a plus). However, it also means that it might be hard to discern the correct assignment of read groups with very little data.

For each aligned base that overlaps a known genotype, we calculate the probability the probability that it was derived from a particular known genotype. This comparison considers only bases that overlap previously known genotypes and that meet the base quality and mapping quality thresholds.

Each individual in a pedigree has a different combination of genotypes, and bamGenotypeCheck will systematically search for the individual whose genotypes best match the observed read data.

For more about the technical details, see the page [[Verifying Sample Identities - Implementation]]

== TODO ==

Main Page

2010-06-22T14:48:13Z

Pha: /* Sequence Analysis Tools */ change name

{| style="width:100%; background:#fcfcfc; margin-top:1.2em; border:1px solid #ccc;" |
| style="width:100%; text-align:center; white-space:nowrap; color:#000;" |
<div style="font-size:162%; border:none; margin:0; padding:.1em; color:#000;">Abecasis Group Wiki</div>
|}

[[Image:2009.08_Group_Retreat_Photo.jpg|400px|center|Group Photo]]

== Welcome! ==

Welcome to our brand new wiki!

If you would like to contribute, [[Special:UserLogin|log-in]] or [http://csgwiki.sph.umich.edu/index.php?title=Special:UserLogin&type=signup create an account]. We recommend using your e-mail address or Michigan uniqname as your user id.

For basic instructions, see [http://en.wikipedia.org/wiki/Wikipedia:Tutorial the Wikipedia Tutorial].

== Sequence Analysis Tools ==

We are developing [[Software|software tools]] for the analysis of next generation sequence data.

# Read Mapping with [[Karma]]
# Variant Calling with [[glfSingle]], [[glfTrio]] and [[glfMultiples]]
# Variant Annotations using [[vcfCodingSnps]]
# Quality control using [[C++ Executable: fastQValidator|FastQValidator]], and [[BamGenotypeCheck]]
# File Readers using [[C++ Library: libbam|BamFile]] and [[C++ Library: libfqf|FastQFile]]
# File Conversions using [[bam2FastQ]]

== Projects ==

[[SardiNIA]] - The SardiNIA longitudinal study of aging.

[[1000 Genomes Project Pilot 1 SNP Calling]]

== Learn Genetics ==

Faculty in the group teach in a variety of formal and informal settings. [[Class Notes|Class notes]] and relevant discussion are archived here.

==Abecasis Group Meeting==

2:30-4:00 p.m. on Fridays, room 2610, SPH I.
Please note April 9, please use room 2690, SPH I.

BamGenotypeCheck

2010-06-22T14:47:38Z

Pha: rename again

'''bamGenotypeCheck''' is a program that verifies whether the reads in particular file match previously known genotypes for an individual (or group of individuals).

== Usage ==

A key step in any genetic analysis is to verify whether data being generated matches expectations. This program checks whether reads in a BAM file match previous genotypes for a specific sample.

Using a mathematical model that relates observed sequence reads to an hypothetical true genotype, bamGenotypeCheck tries to decide whether sequence reads match a particular individual or are more likely to be contaminated (including a small proportion of foreign DNA), derived from a closely related individual, or derived from a completely different individual.

== Basic Usage Example ==

Here is a typical command line:

bamGenotypeCheck -r /data/local/ref/karma.ref/human.g1k.v37.fa \
-k BAMfiles.txt -p test.ped -d test.dat -m test.map

== Command Line Options ==

=== Input Files ===

-r ''FASTA format genome reference''
-a ''allele Frequency file''
-p ''pedigree file in [[MERLIN format]]''
-d ''data file in [[MERLIN format]]''
-m ''map file in [[MERLIN format]]''

-k ''a list of BAM files to check''
-c [int] ''stop after reading [int] filtered sequence reads''
-C [int] ''stop after reading [int] reads, filtered or not''

=== Output Options ===

-v ''verbose output''

=== Filtering ===

-b [int] ''exclude bases with quality less than [int]''
-M [int] ''exclude reads with map quality less than [int]''
-F [int] ''set custom BAM flags filter (not implemented at the moment)''

=== Other Options ===

-e [float] '' set minimum error base error to [float]''

== Principle of Operation ==

Each read group in a BAM file is evaluated independently. This means that in file with multiple read groups, problems will be flagged at the read group level (a plus). However, it also means that it might be hard to discern the correct assignment of read groups with very little data.

For each aligned base that overlaps a known genotype, we calculate the probability the probability that it was derived from a particular known genotype. This comparison considers only bases that overlap previously known genotypes and that meet the base quality and mapping quality thresholds.

Each individual in a pedigree has a different combination of genotypes, and bamGenotypeCheck will systematically search for the individual whose genotypes best match the observed read data.

For more about the technical details, see the page [[Verifying Sample Identities - Implementation]]

== TODO ==

GenotypeIDcheck

2010-06-22T14:46:40Z

Pha: moved GenotypeIDcheck to BamGenotypeCheck: a better name was suggested

#REDIRECT [[BamGenotypeCheck]]

BamGenotypeCheck

2010-06-22T14:46:40Z

Pha: moved GenotypeIDcheck to BamGenotypeCheck: a better name was suggested

'''genotypeIdCheck''' is a program that verifies whether the reads in particular file match previously known genotypes for an individual (or group of individuals).

== Usage ==

A key step in any genetic analysis is to verify whether data being generated matches expectations. This program checks whether reads in a BAM file match previous genotypes for a specific sample.

Using a mathematical model that relates observed sequence reads to an hypothetical true genotype, genotypeIdCheck tries to decide whether sequence reads match a particular individual or are more likely to be contaminated (including a small proportion of foreign DNA), derived from a closely related individual, or derived from a completely different individual.

== Basic Usage Example ==

Here is a typical command line:

genotypeIDcheck -r /data/local/ref/karma.ref/human.g1k.v37.fa \
-k BAMfiles.txt -p test.ped -d test.dat -m test.map

== Command Line Options ==

=== Input Files ===

-r ''FASTA format genome reference''
-a ''allele Frequency file''
-p ''pedigree file in [[MERLIN format]]''
-d ''data file in [[MERLIN format]]''
-m ''map file in [[MERLIN format]]''

-k ''a list of BAM files to check''
-c [int] ''stop after reading [int] filtered sequence reads''
-C [int] ''stop after reading [int] reads, filtered or not''

=== Output Options ===

-v ''verbose output''

=== Filtering ===

-b [int] ''exclude bases with quality less than [int]''
-M [int] ''exclude reads with map quality less than [int]''
-F [int] ''set custom BAM flags filter (not implemented at the moment)''

=== Other Options ===

-e [float] '' set minimum error base error to [float]''

== Principle of Operation ==

Each read group in a BAM file is evaluated independently. This means that in file with multiple read groups, problems will be flagged at the read group level (a plus). However, it also means that it might be hard to discern the correct assignment of read groups with very little data.

For each aligned base that overlaps a known genotype, we calculate the probability the probability that it was derived from a particular known genotype. This comparison considers only bases that overlap previously known genotypes and that meet the base quality and mapping quality thresholds.

Each individual in a pedigree has a different combination of genotypes, and genotypeIdCheck will systematically search for the individual whose genotypes best match the observed read data.

For more about the technical details, see the page [[Verifying Sample Identities - Implementation]]

== TODO ==

Software

2010-06-21T15:19:07Z

Pha: change layout a little

= Software Page Overview =

This page points to downloads, documentation, and papers for software that is written here at the [http://genome.sph.umich.edu Center for Statistical Genetics]

= [[Read Mapping]] =

==[[Karma|Karma]]==
Our fast short read aligner, which generates [[Mapping Quality Scores]]

==[[Karma-colorspace|Karma-ColorSpace]]==
QUICKSTART on mapping color space reads

==[[Examples|Examples]]==
Sample command lines with discussion

==[[MapabilityScores]]==
Definitions of various mappability scores adopted at UCSC genome browser.

==Evaluation of Mappers==
[[baseQualityCheck]] is a mature tool to calculate the observed base quality vs. empirical base quality.

= Variant Calling =

==[[glfSingle]]==
Variant calling for a single, deeply sequenced individual

==[[glfTrio]]==
Variant calling for a single, deeply sequenced nuclear family with two parents and one child

==[[glfMultiples]]==
Variant calling for multiple, unrelated individuals

= Variant Annotation =

==[[vcfCodingSnps]]==
Annotate coding variants in a VCF file.

= Quality Control Utilities =

== Validators ==

[[C++ Executable: fastQValidator|FastQValidator]] -- Check that a FASTQ file conforms to specification.

[[GenotypeIDcheck]] -- Check that mapped reads are consistent with known genotypes for each individual.

== File Readers ==

[[C++ Library: libbam|BamFile]] -- Reads a BAM/SAM file.

[[C++ Library: libfqf|FastQFile]] -- Read a FASTQ file sequence by sequence. Validating the sequence as it is read.

== File Conversion ==

[[bam2FastQ]] -- Convert BAM files into FastQ files

= [[Links to Sequence Analysis Tools|Other Useful Links]] =

BamGenotypeCheck

2010-06-21T15:06:37Z

Pha: total rewrite - new program

== Why genotypeIDcheck? ==

When sequencing data arrives from the sequencer machines, several types of errors or contaminations can creep in.

The most serious error that can occur is that the actual sample name for the sample gets swapped with a different sample. In this case, what we think is one person's sequencing data is actually another. Although they are confidential, in some cases, they are intended to be part of a family, and such an error causes problems in later analysis.

Also, since we typically genotype the same samples along with running a full sequence of them, we want to know that the genotyping information matches the sequence, otherwise again, we will find errors in later analysis.

Another set of problems occurs when the actual biological sample is contaminated in some way. It could be the case, for example, that certain common laboratory contaminates will be included in the sequencing data - knowing approximately how much of such contaminants exist (e.g. yeast) is useful.

So, the program genotypeIDcheck was written in order to discover if problems such as these exist, and if so, report them to the user.

To work, genotypeIDcheck needs a) a list of sorted, calibrated BAM files to check, b) a genome reference to check against, c) a pedigree that describes the samples being checked.

On output, genotypeIDcheck shows, for each read group in each BAM file, the individual in the pedigree that comes statistically closest to matching the reads in the read group, and the degree to how close that relationship is. In addition, marker coverage of three different classes of markers is shown for that individual as well.

== Basic Usage Example ==

Here is an example of how genotypeIDcheck works:

genotypeIDcheck -r /data/local/ref/karma.ref/human.g1k.v37.fa \
-k BAMfiles.txt -p test.ped -d test.dat -m test.map

== Command Line Options ==

=== Input Files ===

-r ''KARMA genome reference''
-k ''a filename that contains a list of BAM files to check''
-a ''pedigree Allele Frequency file''
-p ''pedigree .ped file''
-d ''pedigree .dat file''
-m ''pedigree .map file''

=== Output Options ===

-c [int] ''stop after reading [int] filtered sequence reads''
-C [int] ''stop after reading [int] reads, filtered or not''
-v ''verbose output''

=== Filtering ===

-b [int] ''exclude marker positions with base quality less than [int]''
-M [int] ''exclude all reads with map quality less than [int]''
-F [int] ''set custom BAM flags filter (not implemented at the moment)''

=== Other Options ===

-e [float] '' set minimum error estimate to [float]''

== Principle of Operation ==

For computational and output purposes, we consider each read group sample in each BAM file to be distinct from all others.

For each aligned sample, we calculate the probability that the sample is from each individual in the pedigree according to five different probabilities of being identical by descent. So from the base quality (again, assuming calibrated base qualities), the given base in the read, the marker corresponding reference base and genotype chip read data, you can compute the probability of that individual being the one in the sample for the given probability of IBD.

What you actually want is the multiple of theses probabilities taken across all reads in the sample. This becomes impractically small, so we instead sum the logs of the probabilities.

After the sample is read, the sample in the pedigree that contains the highest log-sum of computed Pibd values is assumed to have the strongest relationship to that sample. The ID of the corresponding pedigree sample name is printed, the sample name indicated in the read group is printed, and coverage information is also printed.

For more about the mathematical details, see the page [[Verifying Sample Identities - Implementation]]

== TODO ==

Main Page

2010-06-21T14:23:59Z

Pha: rename BAMidentity check to genotypeIDcheck to match new program name

{| style="width:100%; background:#fcfcfc; margin-top:1.2em; border:1px solid #ccc;" |
| style="width:100%; text-align:center; white-space:nowrap; color:#000;" |
<div style="font-size:162%; border:none; margin:0; padding:.1em; color:#000;">Abecasis Group Wiki</div>
|}

[[Image:2009.08_Group_Retreat_Photo.jpg|400px|center|Group Photo]]

== Welcome! ==

Welcome to our brand new wiki!

If you would like to contribute, [[Special:UserLogin|log-in]] or [http://csgwiki.sph.umich.edu/index.php?title=Special:UserLogin&type=signup create an account]. We recommend using your e-mail address or Michigan uniqname as your user id.

For basic instructions, see [http://en.wikipedia.org/wiki/Wikipedia:Tutorial the Wikipedia Tutorial].

== Sequence Analysis Tools ==

We are developing [[Software|software tools]] for the analysis of next generation sequence data.

# Read Mapping with [[Karma]]
# Variant Calling with [[glfSingle]], [[glfTrio]] and [[glfMultiples]]
# Variant Annotations using [[vcfCodingSnps]]
# Quality control using [[C++ Executable: fastQValidator|FastQValidator]], and [[GenotypeIDcheck]]
# File Readers using [[C++ Library: libbam|BamFile]] and [[C++ Library: libfqf|FastQFile]]
# File Conversions using [[bam2FastQ]]

== Projects ==

[[SardiNIA]] - The SardiNIA longitudinal study of aging.

[[1000 Genomes Project Pilot 1 SNP Calling]]

== Learn Genetics ==

Faculty in the group teach in a variety of formal and informal settings. [[Class Notes|Class notes]] and relevant discussion are archived here.

==Abecasis Group Meeting==

2:30-4:00 p.m. on Fridays, room 2610, SPH I.
Please note April 9, please use room 2690, SPH I.

BamIdentityCheck

2010-06-21T14:22:32Z

Pha: moved BamIdentityCheck to GenotypeIDcheck: now matches the new program name

#REDIRECT [[GenotypeIDcheck]]

BamGenotypeCheck

2010-06-21T14:22:32Z

Pha: moved BamIdentityCheck to GenotypeIDcheck: now matches the new program name

== Basic Usage Example ==

Here is an example of how laneCheck works:

lanecheck --referencegenome NCBI36.fa --dbSNPfile dbSNP.txt
--lanefile lane.lst --pedfile test.ped --datfile test.dat --mapfile test.map
--prefix result

== Command Line Options ==

=== Input Files ===

--referencegenome ''referencegenome file''
--dbSNPfile ''optional'' ''two-column'' ''dbsnp position file, will provide more accurate background mismatch rate if excluding dbSNP positions (e.g. 5 123456)''
--lanefile ''a list of lane file with path''
--pedfile ''genotype information of the samples for checking ''
--datfile ''a companion data file for pedigree file (each row: M snpname, e.g. M rs1234)''
--mapfile ''a companion data file for pedigree file (each row: chr snpname pos, e.g. 5 rs1234 56789)''

=== Basic Output Options ===

--prefix ''specify the prefix name of the output file''

=== Filtering  ===

--minmapquality ''reads with with mapquality falling below this threshold will be excluded''
--genocount ''the maximum number of genotypes compared''
--verbose ''print out detailed information for each hapmap position compared''
--coverage ''print out the proportion of markers in the map file covered by at least one read''
--countbysite ''print out detailed mismatch counts for each base compared''

=== Other Options ===

--memorymap ''use memory map technique for efficient memory sharing of reference genome file''

== Principle of Operation: ==

The overall procedure is that the genotype identity checking program compares internal evidence from the sequence reads themselves to reference genotype information for a panel of candidate individuals. In the case of 1000 Genomes pilot data, these are HapMap genotypes from the same Coriell cell lines that are being sequenced. For each combination of [sequencing run x candidate individual] the program calculates the observed rate of mismatches at both "informative" and "background" locations and reports as "excess mismatch rate"

excess rate = (informative rate - background rate).

"Informative" locations are those where the candidate individual is homozygous, according to the HapMap genotype information, and base calls are compared to the HapMap homozygous allele, rather than to the genome reference sequence. "Background" locations are all sites not known to be polymorphic and not recorded in dbSNP if provided.  A relative high background rate suggests possible problems in sample preparation or read mapping process. 

 

== TODO ==

1. Separate the results by "Read group classifier".

The mapped .bam file may contains sequence data from different instrument runs. The read identifiers often are dot or colon-separated strings of the form 'run_name<sep>read_number'. The 'run_name' may be either an SRR / ERR identifier or the sequencing center's own alpha-numeric internal run identifier. Allow users to input extended regular expression such as '\(^[^.:]+\)[.:].*' hich matches just the part of each read identifier that is common to all reads from one instrument run and which differs between instrument runs.

2. Use model based approach to calculate probability of lane coming from the claimed individual in the index file given a pool of individuals.

BAM Review Action Items

2010-06-08T13:49:07Z

Pha: add new name options, spelling corrections

* Move the examples from the SamFile wiki page to their own page
** include links from the main library page and the SamFile page.
** look into why the one example have two if checks on SamIn status
* Create 1 library for all of our library code rather than having libcsg, libbam, libfqf separated.
** What should this library be called?
*** libdna
*** libdna++
*** libsequence++
*** libDNA
*** libgenotype
* Add an option by class that says whether or not to abort on failure. (or even an option on each method)
** This allows calling code to set that option and then not have to check for failures since the code it calls would abort on a failure.
** Could/should this be achieved using exceptions? User can decide to catch them or let them terminate the program.
*SamFile add a constructor that takes the filename.
** Also takes a flag to indicate open for read/write.
** Also have 2 subclasses one that opens for read, one for write: SamReadFile, SamWriteFile? Or SamFileRead, SamFileWrite?
* Add a function that says: skipInvalidRecords, validateRecords, etc.
** That way, ReadRecord will keep reading records until a valid/parseable one is found.
*SamFileHeader::setTag - instead of having separate ones for PG, RG, etc, have a generic one that takes as a parameter which one it is.
** KeyID, then Value as parameters....(keyID first, then value)
* SamFileHEader::setProgramName, etc...have specific methods for setting fields so users don't need to know the specific tags, etc used for certain values in the header.
** KeyID, then Value as parameters....(keyID first, then value)
* BAM write utility could add a PG field with default settings (user could specify alternate settings) when it writes a file.
* Future methods to add:
** SamFile::setReadSection(const std::string& refName) - take in the reference name by string since that is what most people will know.
*** "" would indicate the ones not associated with a reference.

Verifying Sample Identities - Implementation

2010-04-20T18:18:06Z

Pha: /* Mathematical Details */

== Principle ==

We should be able to verify that the right sample has been sequenced by comparing base calls in a read to known genotypes for a sample. If the sample has been sequenced correctly, the base calls should match previously known genotypes. If the wrong sample has been sequenced, we will see quite a bit more mismatches.

== Mathematical Details ==

For each sample, we would like to calculate the likelihood of a set of reads assuming that we sequenced the correct sample, assuming we sequenced a sample related to the correct sample, or assuming we sequenced an incorrect sample. We would then like to flag samples where it appears likely that the wrong sample has been sequenced.

If we have a list of bases that overlap a known genotype, we can will
describe the probability of a matching of mismatching base using the
following notation:

{| width="100%" cellspacing="1" cellpadding="1" border="1" summary="Summary of Variables Used Below"
|+ Notation
|-
| Variable
| Definition
|-
| ''A/A''
| Previously known genotype; we only consider homozygous sites.
|-
| ''P''''A''
| Frequency of allele A in the population
|-
| ''P''''ibd''
| Probability that the sequenced sample and the target sample share a chromosome. This should be 1.0 when we have sequenced the correct sample and 0.0 if we sequence an unrelated sample. If we sequence a related sample (e.g. a parent or sibling of the target sample), we will see intermediate values.
|-
| <math>\epsilon</math>
| Estimate error rate for the current base in the sequence data.
|}

Then, the probabilities of interest are:

<math>
P(match) = P_{ibd} (1 - \epsilon) + (1 - P_{ibd}) \epsilon
</math>

<math>
P(no match) = P_{ibd} \epsilon + (1 - P_{ibd}) \epsilon = \epsilon ?
</math>

For now, genotypeIDcheck uses:

<math>
P(no match) = 1 - P(match)
</math>

If we also wish to use heterozygous sites, rather than limiting our comparison to reference homozygous sites, we could use:

<math>
P(match) = P(no match) = P_{ibd} 0.5 + (1.0 - P_{ibd}) P_A
</math>

For any given value of <math>P_{ibd}</math> we can calculate this quantity for each read that overlaps a site with a known homozygous genotype. In addition, we can take the product of this quantity across all sites examined -- because this product is likely to be very small, we actually sum the <math>log</math>s of the appropriate quantities rather than multiplying them together.

To decide if we have sequenced the correct sample, we should do the following:

a) Evaluate this log-sum assuming <math>P_{ibd} = 1.0</math>. This assumes that we have sequenced the target sample.

b) Evaluate this log-sum assuming <math>P_{ibd} = 0.0</math>. This assumes we sequenced a different sample, unrelated to the target.

c) Evaluate this log-sum assuming <math>P_{ibd} = 0.5</math>. This assumes we sequenced a sample that shares half the genome with the target sample, perhaps because it is a sibling or parent of the target sample.

d) If desired, evaluate the same log-sum for other intermediate values of <math>P_{ibd}</math>. It may be interesting to set <math>P_{ibd} = 0.95</math> to allow for 5% of reads that are derived from a different sample, for example, due to contamination. It may be interesting to set <math>P_{ibd} = 0.05</math> to consider more distant relatives.

Once the result of evaluating a), b), c) and d) are available, we can decide if the target sample has been sequenced. Sequencing the target sample will mean that the log-sum in a) is the largest. Sequencing a parent or offspring of the target sample will maximize c). Sequencing a completely incorrect sample will maximize b).

If all the log-sums are very similar, then we don't have enough information to make a clear cut decision. Typically, we thousands of genetic markers from a typical SNP chip and whole genome shotgun sequence data, most decisions should be very clear cut.

== Implementation Details ==

After loading genotypes, we generate a genome mask for each position. There are three outcomes of interest:

; Known Genotypes
: These are sites where we have a previously observed a genotype call and where we will be evaluating match / mismatch rates to determine sample identity.

; dbSNP sites
: These are sites that are known to vary among individuals, but for which a known genotype is not available.

; Background sites
: These are all other sites and can be used to estimate the <math>\epsilon</math> error rate parameter.

C++ Class: Pedigree Example

2010-04-13T20:37:06Z

Pha: /* Example */

== Example ==

<source lang="cpp">

#include "Pedigree.h"

int main(int argc, char ** argv);
{
Pedigree ped;

// The data file contains a description of the contents of the
// pedigree file, including for example, a list of marker and
// trait names
ped.Prepare("pedigree.dat");

// The pedigree file contains a list of individuals, stored one
// per row, with specific information about each individual as
// detailed in the data file.
ped.Load("pedigree.ped");

printf("Loaded %d individuals\n", ped.count);
printf("Loaded %d markers\n", ped.markerCount);

// Print out names of the first 10 individuals
for (int i = 0; i < max(ped.count, 10); i++)
printf("Individual #1 is labeled %s\n", (const char *) ped[i].pid);

if (ped.markerCount == 0)
{
printf("No marker data available!\n");
return 0;
}

// Estimate allele frequencies for all markers
ped.EstimateAlleleFrequencies();

// Get some genotype statistics for the first marker
printf("Statistics for marker %s\n", (const char *) ped.markerNames[0]);

// Get information on the first marker
MarkerInfo * info = ped.GetMarkerInfo(0);

// Get information on marker rs12345
MarkerInfo * rs12345 = ped.GetMarkerInfo("rs12345");

// Number of alleles for the first marker
int alleleCount = info.GetAlleleCount();

// Frequency of the first allele
int freq1 = info.freq[1];

// Name of the first allele, typically "A", "C", "G" or "T" for SNPs
String label = info.GetAlleleLabel(1);

// Genotype for a specific individual at a specific marker
ped[0].markers[0].isKnown(); // Test if genotype is known
ped[0].markers[0].isHomozygous(); // Test if genotype is homozygous
ped[0].markers[0][0]; // Index of first allele
ped[0].markers[0][1]; // Index of second allele
}
</source>

C++ Class: Pedigree Example

2010-04-13T20:35:20Z

Pha: /* Example */

C++ Class: Pedigree Example

2010-04-13T20:31:59Z

Pha: /* Example */

C++ Class: Pedigree Example

2010-04-13T20:31:14Z

Pha: Created page with '== Example == <source lang="cpp"> #include "Pedigree.h" int main(int argc, char ** argv); { Pedigree ped; // The data file contains a description of the contents of the…'

Verifying Sample Identities - Implementation

2010-04-13T20:24:18Z

Pha: /* Mathematical Details */

== Principle ==

We should be able to verify that the right sample has been sequenced by comparing base calls in a read to known genotypes for a sample. If the sample has been sequenced correctly, the base calls should match previously known genotypes. If the wrong sample has been sequenced, we will see quite a bit more mismatches.

== Mathematical Details ==

For each sample, we would like to calculate the likelihood of a set of reads assuming that we sequenced the correct sample, assuming we sequenced a sample related to the correct sample, or assuming we sequenced an incorrect sample. We would then like to flag samples where it appears likely that the wrong sample has been sequenced.

If we have a list of bases that overlap a known genotype, we can will
describe the probability of a matching of mismatching base using the
following notation:

{| width="100%" cellspacing="1" cellpadding="1" border="1" summary="Summary of Variables Used Below"
|+ Notation
|-
| Variable
| Definition
|-
| ''A/A''
| Previously known genotype; we only consider homozygous sites.
|-
| ''P''''A''
| Frequency of allele A in the population
|-
| ''P''''ibd''
| Probability that the sequenced sample and the target sample share a chromosome. This should be 1.0 when we have sequenced the correct sample and 0.0 if we sequence an unrelated sample. If we sequence a related sample (e.g. a parent or sibling of the target sample), we will see intermediate values.
|-
| <math>\epsilon</math>
| Estimate error rate for the current base in the sequence data.
|}

Then, the probabilities of interest are:

<math>
P(match) = P_{ibd} (1 - \epsilon) + (1 - P_{ibd}) \epsilon
</math>

<math>
P(no match) = P_{ibd} \epsilon + (1 - P_{ibd}) \epsilon
</math>

If we also wish to use heterozygous sites, rather than limiting our comparison to reference homozygous sites, we could use:

<math>
P(match) = P(no match) = P_{ibd} 0.5 + (1.0 - P_{ibd}) P_A
</math>

For any given value of <math>P_{ibd}</math> we can calculate this quantity for each read that overlaps a site with a known homozygous genotype. In addition, we can take the product of this quantity across all sites examined -- because this product is likely to be very small, we actually sum the <math>log</math>s of the appropriate quantities rather than multiplying them together.

To decide if we have sequenced the correct sample, we should do the following:

a) Evaluate this log-sum assuming <math>P_{ibd} = 1.0</math>. This assumes that we have sequenced the target sample.

b) Evaluate this log-sum assuming <math>P_{ibd} = 0.0</math>. This assumes we sequenced a different sample, unrelated to the target.

c) Evaluate this log-sum assuming <math>P_{ibd} = 0.5</math>. This assumes we sequenced a sample that shares half the genome with the target sample, perhaps because it is a sibling or parent of the target sample.

d) If desired, evaluate the same log-sum for other intermediate values of <math>P_{ibd}</math>. It may be interesting to set <math>P_{ibd} = 0.95</math> to allow for 5% of reads that are derived from a different sample, for example, due to contamination. It may be interesting to set <math>P_{ibd} = 0.05</math> to consider more distant relatives.

Once the result of evaluating a), b), c) and d) are available, we can decide if the target sample has been sequenced. Sequencing the target sample will mean that the log-sum in a) is the largest. Sequencing a parent or offspring of the target sample will maximize c). Sequencing a completely incorrect sample will maximize b).

If all the log-sums are very similar, then we don't have enough information to make a clear cut decision. Typically, we thousands of genetic markers from a typical SNP chip and whole genome shotgun sequence data, most decisions should be very clear cut.

== Implementation Details ==

After loading genotypes, we generate a genome mask for each position. There are three outcomes of interest:

; Known Genotypes
: These are sites where we have a previously observed a genotype call and where we will be evaluating match / mismatch rates to determine sample identity.

; dbSNP sites
: These are sites that are known to vary among individuals, but for which a known genotype is not available.

; Background sites
: These are all other sites and can be used to estimate the <math>\epsilon</math> error rate parameter.

Verifying Sample Identities - Implementation

2010-04-13T20:20:51Z

Pha: /* Mathematical Details */

== Principle ==

We should be able to verify that the right sample has been sequenced by comparing base calls in a read to known genotypes for a sample. If the sample has been sequenced correctly, the base calls should match previously known genotypes. If the wrong sample has been sequenced, we will see quite a bit more mismatches.

== Mathematical Details ==

For each sample, we would like to calculate the likelihood of a set of reads assuming that we sequenced the correct sample, assuming we sequenced a sample related to the correct sample, or assuming we sequenced an incorrect sample. We would then like to flag samples where it appears likely that the wrong sample has been sequenced.

If we have a list of bases that overlap a known genotype, we can will
describe the probability of a matching of mismatching base using the
following notation:

{| width="100%" cellspacing="1" cellpadding="1" border="1" summary="Summary of Variables Used Below"
|+ Notation
|-
| Variable
| Definition
|-
| ''A/A''
| Previously known genotype; we only consider homozygous sites.
|-
| ''P''''A''
| Frequency of allele A in the population
|-
| ''P''''ibd''
| Probability that the sequenced sample and the target sample share a chromosome. This should be 1.0 when we have sequenced the correct sample and 0.0 if we sequence an unrelated sample. If we sequence a related sample (e.g. a parent or sibling of the target sample), we will see intermediate values.
|-
| <math>\epsilon</math>
| Estimate error rate for the current base in the sequence data.
|}

Then, the probabilities of interest are:

<math>
P(match) = P_{ibd} (1 - \epsilon) + (1 - P_{ibd}) \epsilon
</math>

<math>
P(no match) = P_{ibd} \epsilon + (1 - P_{ibd}) \epsilon
</math>

For any given value of <math>P_{ibd}</math> we can calculate this quantity for each read that overlaps a site with a known homozygous genotype. In addition, we can take the product of this quantity across all sites examined -- because this product is likely to be very small, we actually sum the <math>log</math>s of the appropriate quantities rather than multiplying them together.

To decide if we have sequenced the correct sample, we should do the following:

a) Evaluate this log-sum assuming <math>P_{ibd} = 1.0</math>. This assumes that we have sequenced the target sample.

b) Evaluate this log-sum assuming <math>P_{ibd} = 0.0</math>. This assumes we sequenced a different sample, unrelated to the target.

c) Evaluate this log-sum assuming <math>P_{ibd} = 0.5</math>. This assumes we sequenced a sample that shares half the genome with the target sample, perhaps because it is a sibling or parent of the target sample.

d) If desired, evaluate the same log-sum for other intermediate values of <math>P_{ibd}</math>. It may be interesting to set <math>P_{ibd} = 0.95</math> to allow for 5% of reads that are derived from a different sample, for example, due to contamination. It may be interesting to set <math>P_{ibd} = 0.05</math> to consider more distant relatives.

Once the result of evaluating a), b), c) and d) are available, we can decide if the target sample has been sequenced. Sequencing the target sample will mean that the log-sum in a) is the largest. Sequencing a parent or offspring of the target sample will maximize c). Sequencing a completely incorrect sample will maximize b).

If all the log-sums are very similar, then we don't have enough information to make a clear cut decision. Typically, we thousands of genetic markers from a typical SNP chip and whole genome shotgun sequence data, most decisions should be very clear cut.

== Implementation Details ==

After loading genotypes, we generate a genome mask for each position. There are three outcomes of interest:

; Known Genotypes
: These are sites where we have a previously observed a genotype call and where we will be evaluating match / mismatch rates to determine sample identity.

; dbSNP sites
: These are sites that are known to vary among individuals, but for which a known genotype is not available.

; Background sites
: These are all other sites and can be used to estimate the <math>\epsilon</math> error rate parameter.

Verifying Sample Identities - Implementation

2010-04-13T20:09:19Z

Pha: /* Mathematical Details */

== Principle ==

We should be able to verify that the right sample has been sequenced by comparing base calls in a read to known genotypes for a sample. If the sample has been sequenced correctly, the base calls should match previously known genotypes. If the wrong sample has been sequenced, we will see quite a bit more mismatches.

== Mathematical Details ==

For each sample, we would like to calculate the likelihood of a set of reads assuming that we sequenced the correct sample, assuming we sequenced a sample related to the correct sample, or assuming we sequenced an incorrect sample. We would then like to flag samples where it appears likely that the wrong sample has been sequenced.

If we have a list of bases that overlap a known genotype, we can will
describe the probability of a matching of mismatching base using the
following notation:

{| width="100%" cellspacing="1" cellpadding="1" border="1" summary="Summary of Variables Used Below"
|+ Notation
|-
| Variable
| Definition
|-
| ''A/A''
| Previously known genotype; we only consider homozygous sites.
|-
| ''P''''A''
| Frequency of allele A in the population
|-
| ''P''''ibd''
| Probability that the sequenced sample and the target sample share a chromosome. This should be 1.0 when we have sequenced the correct sample and 0.0 if we sequence an unrelated sample. If we sequence a related sample (e.g. a parent or sibling of the target sample), we will see intermediate values.
|-
| <math>\epsilon</math>
| Estimate error rate for the current base in the sequence data.
|}

Then, the probabilities of interest are:

<math>
P(match) = P_{ibd} (1 - \epsilon) + (1 - P_{ibd}) \epsilon
</math>

<math>
P(no match) = P_{ibd} \epsilon + (1 - P_{ibd}) \epsilon
</math>

For any given value of <math>P_{ibd}</math> we can calculate this quantity for each read that overlaps a site with a known homozygous genotype. In addition, we can take the product of this quantity across all sites examined -- because this product is likely to be very small, we actually sum the <math>log</math>s of the appropriate quantities rather than multiplying them together.

To decide if we have sequenced the correct sample, we should do the following:

a) Evaluate this log-sum assuming <math>P_{ibd} = 1.0</math>. This assumes that we have sequenced the target sample.

b) Evaluate this log-sum assuming <math>P_{ibd} = 0.0</math>. This assumes we sequenced a different sample, unrelated to the target.

c) Evaluate this log-sum assuming <math>P_{ibd} = 0.5</math>. This assumes we sequenced a sample that shares half the genome with the target sample, perhaps because it is a sibling or parent of the target sample.

d) If desired, evaluate the same log-sum for other intermediate values of P_{ibd}. It may be interesting to set <math>P_{ibd} = 0.95</math> to allow for 5% of reads that are derived from a different sample, for example, due to contamination. It may be interesting to set <math>P_{ibd} = 0.05</math> to consider more distant relatives.

Once the result of evaluating a), b), c) and d) are available, we can decide if the target sample has been sequenced. Sequencing the target sample will mean that the log-sum in a) is the largest. Sequencing a parent or offspring of the target sample will maximize c). Sequencing a completely incorrect sample will maximize b).

If all the log-sums are very similar, then we don't have enough information to make a clear cut decision. Typically, we thousands of genetic markers from a typical SNP chip and whole genome shotgun sequence data, most decisions should be very clear cut.

Verifying Sample Identities - Implementation

2010-04-13T20:03:43Z

Pha: /* Mathematical Details */

== Principle ==

We should be able to verify that the right sample has been sequenced by comparing base calls in a read to known genotypes for a sample. If the sample has been sequenced correctly, the base calls should match previously known genotypes. If the wrong sample has been sequenced, we will see quite a bit more mismatches.

== Mathematical Details ==

For each sample, we would like to calculate the likelihood of a set of reads assuming that we sequenced the correct sample, assuming we sequenced a sample related to the correct sample, or assuming we sequenced an incorrect sample. We would then like to flag samples where it appears likely that the wrong sample has been sequenced.

If we have a list of bases that overlap a known genotype, we can will
describe the probability of a matching of mismatching base using the
following notation:

{| width="100%" cellspacing="1" cellpadding="1" border="1" summary="Summary of Variables Used Below"
|+ Notation
|-
| Variable
| Definition
|-
| ''A/A''
| Previously known genotype; we only consider homozygous sites.
|-
| ''P''''A''
| Frequency of allele A in the population
|-
| ''P''''ibd''
| Probability that the sequenced sample and the target sample share a chromosome. This should be 1.0 when we have sequenced the correct sample and 0.0 if we sequence an unrelated sample. If we sequence a related sample (e.g. a parent or sibling of the target sample), we will see intermediate values.
|-
| <math>\epsilon</math>
| Estimate error rate for the current base in the sequence data.
|}

Then, the probabilities of interest are:

<math>
P(match) = P_{ibd} (1 - \epsilon) + (1 - P_{ibd}) \epsilon
</math>

<math>
P(no match) = P_{ibd} \epsilon + (1 - P_{ibd}) \epsilon
</math>

For any given value of P_{ibd} we can calculate this quantity for each read that overlaps a site with a known homozygous genotype. In addition, we can take the product of this quantity across all sites examined -- because this product is likely to be very small, we actually sum the <math>log</math>s of the appropriate quantities rather than multiplying them together.

To decide if we have sequenced the correct sample, we should do the following:

a) Evaluate this log-sum assuming P_{ibd} = 1.0. This assumes that we have sequenced the target sample.

b) Evaluate this log-sum assuming P_{ibd} = 0.0. This assumes we sequenced a different sample, unrelated to the target.

c) Evalute this log-sum assuming P_{ibd} = 0.5. This assumes we sequenced a sample that shares half the genome with the target sample, perhaps because it is a sibling or parent of the target sample.

d) If desired, evaluate the same log-sum for other intermediate values of P_{ibd}. It may be interesting to set P_{ibd} = 0.95 to allow for 5% of reads that are derived from a different sample, for example, due to contamination. It may be interesting to set P_{ibd} = 0.05 to consider more distant relatives.

Once the result of evaluating a), b), c) and d) are available, we can decide if the target sample has been sequenced. Sequencing the target sample will mean that the log-sum in a) is the largest. Sequencing a parent or offspring of the target sample will maximize c). Sequencing a completely incorrect sample will maximize b).

If all the log-sums are very similar, then we don't have enough information to make a clear cut decision. Typically, we thousands of genetic markers from a typical SNP chip and whole genome shotgun sequence data, most decisions should be very clear cut.

Verifying Sample Identities - Implementation

2010-04-13T19:52:58Z

Pha: /* Mathematical Details */

Verifying Sample Identities - Implementation

2010-04-13T19:49:48Z

Pha: /* Mathematical Details */

Verifying Sample Identities - Implementation

2010-04-13T19:49:15Z

Pha:

Verifying Sample Identities - Implementation

2010-04-13T19:47:09Z

Pha: Created page with '== Principle == We should be able to verify that the right sample has been sequenced by comparing base calls in a read to known genotypes for a sample. If the sample has been s…'

Karma

2010-04-09T14:30:09Z

Pha: /* Karma CHANGELOG */ fill in 0.9.0 changes

[[Category:Software]]

'''K-tuple Alignment with Rapid Matching Algorithm'''

Karma uses an existing reference to align short reads, such as those generated by Illumina sequencers.

Primary features:

# High performance, high sensitivity
# Large and small gap detection by default
# Multiple gaps per read by default
# Single or paired end reads
# No read length limit
# Quality scores are used to assess quality of maps
# All potential locations are examined exhaustively, none are omitted
# Reasonable memory per CPU ratio on high core count machines

The current version, 0.9.0, is optimized to rapidly map base space reads from Illumina sequencers.

Color space and LS454 sequence alignments are not currently supported. These features will return in Karma 0.9.1.

= Download Karma =

To get a copy go to [http://www.sph.umich.edu/csg/pha/karma/download/ Karma Download]

= Build Karma =

Karma is designed to be reasonably portable.

However, since development occurs only on Ubuntu 9.10 x86 and x64 platforms, there are likely other portability issues.

We support Karma only on Ubuntu 9.10 on 64-bit processors.

== Dependencies ==

Karma requires that the following debian packages be installed on the host Linux machine:

# libssl-dev
# zlib1g-dev

Without these installed, Karma will not build.

== Building ==

Assuming the karma tar file is named karma.tgz, do the following

tar xvzf karma.tgz
cd karma-0.9
make
mkdir ~/bin
cp karma/karma ~/bin

Alternatively, if you want to share the karma binary install it in /usr/local/bin/karma.

== Testing the build ==

To test karma, go to the build tree subdirectory named ''karma'', and type the command:

make test

The test script builds a reference for the small phiX genome, then runs single end as well as paired end alignments. It compares the results of that with known results. Differences are printed to the console, and currently look something like this:

<pre>
diff phiX.sam.good phiX.sam
3c3
< @RG DT:2010-04-08T17:29Z ID:boingboing SM:NA12345
---
> @RG DT:2010-04-08T18:13Z ID:boingboing SM:NA12345
</pre>

Any differences greater than that are an error and need to be fixed by the author.

= Normal Workflow =

Karma works using a set of index and hash files created from an existing reference. Once created, this set of reference index and hash files must always be specified in the command line when aligning reads.

In concept, the simplest workflow is to first create a reference index using ''karma create'', then align reads using ''karma map''. You only have to build the index and hash once.

Because the reference can be large, and because Karma will share the reference among many running instances of Karma, it is useful to put well known references in a common location readily accessible to you and your collaborators.

= Build reference index and hash =

Building a reference index and hash with Karma is straightforward, but because it is time consuming for longer genomes, you typically save the reference index between runs.

The simplest example for creating a reference and index using a wordsize of 11-mer words is:

karma create -i -w 11 phiX.fa

More generally, three primary parameters are necessary for building a Karma reference index:

# a boolean flag indicating base or color space
# the index table word occurrence cutoff value
# the word size

Although the input reference is always expected to be base space and in [http://en.wikipedia.org/wiki/FASTA_format FASTA] format, the binary version of the reference, and the corresponding index and hash files, can be in either color space (ABI SOLiD) or base space (Illumina or LS454). For a given reference [http://en.wikipedia.org/wiki/FASTA_format FASTA] file, you may have either a color or base space binary reference, as well as either color or base space index/hash files, any in varying word sizes or occurrence cutoffs.

Because the index and hash files are dependent on the occurrence cutoff parameter and the word size, the output files created by karma have those values in the file name. This allows you to create a variety of index/hash tables, depending on your expected use (ABI SOLiD, in particular, is sensitive to read length).

== Options for building reference index and hash ==

-r ''reference'' Reference file in [http://en.wikipedia.org/wiki/FASTA_format FASTA] format
-w ''word size'' Word size for index and hash (default 15, typically 10-16)
-O ''occurrence cutoff'' Upper count of number of word positions to store in word positions table (default 5000)
-c Creates a color space reference and index/hash
-i Create the index and hash as well as the binary reference

= Aligning Reads =

Aligning reads to the reference is easy:

karma map -r phiX.fa -w 11 phiX.fastq

or for paired reads:

karma map -r phiX.fa -w 11 phiX-mate1.fastq phiX-mate2.fastq

In both of the above examples, the -r option names the reference originally used to build the index/hash, and the -w 11 specifies that we are using the index/hash built for 11-mer words. Although you can use the default word size of 15 for phiX, the index is 4^15 * 4 = 4GBytes, so a shorter word size is prudent.

Since Karma uses the word size and occurrence cutoff to help construct the actual index and hash filenames, you must specify them the same way you did when you created the reference index and hash.

== Aligning Reads (Illumina) ==

Karma is set up so that the default options work well for mapping Illumina reads to the Human genome.

== Aligning Reads (ABI SOLiD) ==

Karma has been designed to align color space reads. However, in Karma 0.9.0, this functionality is not working.

== Aligning Reads (LS 454) ==

Karma has been designed to align LS 454 reads. However, in Karma 0.9.0, this functionality is not working.

= Karma Performance Tuning =

There are four components to the Karma index and hash. A pure index array, based on an N-mer word index. This is used as a pointer into a word positions table, which is an ordered list of genome positions in which that N-mer word appears. There is a cap called the ''occurrence cutoff'', which once exceeded, causes that index word to be marked as a high repeat pattern. Once marked as high repeat, the N-mer word is instead combined with both the N-mer word preceding it, as well as the N-mer word succeeding it to create a 2 * N-mer word hash key. Two hash tables are populated, a left and a right hash. These are then used when that pattern is found in a read.

== Index Word Size ==

Choosing an appropriate word size for larger genome is critical to performance. The easiest case is for Illumina base space reads with the human genome (3Gbases), where the default 15-mer word size is fine.

For smaller genomes, consider using a smaller word size. Genomes smaller than a few million bases should be perfectly fine with a word size of 11 or 12.

Since the primary index table into the word positions table is 2^(wordsize) * 4 bytes, it can grow large rapidly. All else being equal, a smaller word size leads to longer sets of word positions for each index value. Each increment of word size approximately quadruples storage requirements, and halves runtime. Similarly, each decrement of word size reduces the index table size by 75%, and doubles runtime. These approximations are old, but serve a useful rule of thumb.

For ABI SOLiD reads, the word size is critical, due to the shorter length of reads as compared to Illumina or LS 454.

The optimal minimum word size is chosen such that it is 1/4 the minimum expected average read length. It also must be chosen to be 1/2 the minimum expected read length, since at least 2 full words must exist in the read.

So for 48-mer reads, a reasonable value of word size is 12. Although the base space default of 15 is fine, too, Karma is able to take advantage of a higher number of index words per read, yielding substantial speedups even with the shorter read. Similarly, 52-mer reads would map better with a 13-mer word size, and 56-mer reads would map best with a 14-mer word size.

== Occurrence Cutoff ==

The occurrence cutoff value determines how quickly an N-mer pattern is declared to be ''high repeat'' and left out of the index in favor of a hash. The default value of 5000 seems adequate for Illumina reads with the human genome. If ultimate performance is necessary, some experimentation is called for with this value.

== Shared Memory ==

Karma uses memory mapped files to share the potentially large reference index and hash data structures.

Karma uses this to great effect on our 8 processors with hyperthreading enabled. 16 copies of karma can share one reference index and hash, yielding a very acceptable memory per CPU ratio of around 1GB/CPU.

A problem with large reference index and hash data structures is that they are more prone to being paged out. On a shared machine that is being used extensively even just simple disk I/O, memory pages are being reclaimed such that Karma will become swapped out.

While Karma can recover on its own, it is best to either run in a production manner on dedicated machines, or to run a program such as the utility ''mapfile'' found in the utilities sub-folder. This program continually touches each page of the data structures in sequential order, forcing them to the head of the disk buffer pool, so they don't get aged out of the queue.

= Modifying the Reference Header =

''NB: This feature is not yet complete''

To facilitate SAM RG values being set automatically in a production environment, we keep a header in the binary version of the reference. The header can be viewed and edited using the header subcommands here.

To view the header:

karma header -r phiX.fa

To view and edit the header:

karma header -r phiX.fa -e

= Other test and check capabilities =

Due to the size and complexity of Karma input, output and index files, various checks and tests are useful, so we include some diagnostics capabilities:

Tests for external files:

karma check [options...] file.bam file.fastq file.sam file.fa file.umfa

Tests internal to Karma:
karma test [options...]
-d -> debug
-s [int] -> set random number seed [12345]

= Karma File structure =

Upon successfully building references, you will obtain a list of reference files like below:

{| style="width: 571px; height: 288px" border="1" cellspacing="1" cellpadding="1" width="571"
|-
|
|
Base Space

| Color Space
|-
|
Reference genome

|
NCBI37-bs.umfa

| NCBI37-cs.umfa
|-
|
Word Index

|
NCBI37-bs.15.5000.umwiwp

NCBI37-bs.15.5000.umwihi

|
NCBI37-cs.15.5000.umwiwp

NCBI37-cs.15.5000.umwihi

|-
|
Word Hash (Left)

|
NCBI37-bs.15.5000.umwhl

| NCBI37-cs.15.5000.umwhl
|-
|
Word Hash (Right)

|
NCBI37-bs.15.5000.umwhr

| NCBI37-cs.15.5000.umwhr
|}

 

 

= Karma TODO List =

# command line help is muddled up - UserOptions.h needs work
# color space read handling needs to be re-integrated and tested
# LS 454 needs to be tested, and the code adapted
# pre-process some number of records to establish an appropriate max insert size
# finish reference header view/edit code
# investigate and document maximum memory use during ''create'' sub-command
# finish and improve check and test commands

= Karma CHANGELOG =

* Karma 0.9.0
* reference may now contain an arbitrary number of chromosomes
* local re-alignment is drastically improved (handle small and large gaps better)
* command line is re-vamped - now easier to use
* create/naming/using reference index and hashes is easier
* Karma 0.8.8S
* add first version of local re-alignment
* bump max number of chromosomes to 200
* we no longer do Smith-Waterman on each candidate location
* Karma 0.8.8
* Karma 0.8.6

= Other useful links =

[http://lh3lh3.users.sourceforge.net/bioinfo.shtml Heng Li's thoughts about aligners]

[http://lh3lh3.users.sourceforge.net/udb.shtml Benchmark of Dictionary Structures]

Karma

2010-04-09T14:01:41Z

Pha: add TODO items and begin listing change history

[[Category:Software]]

'''K-tuple Alignment with Rapid Matching Algorithm'''

Karma uses an existing reference to align short reads, such as those generated by Illumina sequencers.

Primary features:

# High performance, high sensitivity
# Large and small gap detection by default
# Multiple gaps per read by default
# Single or paired end reads
# No read length limit
# Quality scores are used to assess quality of maps
# All potential locations are examined exhaustively, none are omitted
# Reasonable memory per CPU ratio on high core count machines

The current version, 0.9.0, is optimized to rapidly map base space reads from Illumina sequencers.

Color space and LS454 sequence alignments are not currently supported. These features will return in Karma 0.9.1.

= Download Karma =

To get a copy go to [http://www.sph.umich.edu/csg/pha/karma/download/ Karma Download]

= Build Karma =

Karma is designed to be reasonably portable.

However, since development occurs only on Ubuntu 9.10 x86 and x64 platforms, there are likely other portability issues.

We support Karma only on Ubuntu 9.10 on 64-bit processors.

== Dependencies ==

Karma requires that the following debian packages be installed on the host Linux machine:

# libssl-dev
# zlib1g-dev

Without these installed, Karma will not build.

== Building ==

Assuming the karma tar file is named karma.tgz, do the following

tar xvzf karma.tgz
cd karma-0.9
make
mkdir ~/bin
cp karma/karma ~/bin

Alternatively, if you want to share the karma binary install it in /usr/local/bin/karma.

== Testing the build ==

To test karma, go to the build tree subdirectory named ''karma'', and type the command:

make test

The test script builds a reference for the small phiX genome, then runs single end as well as paired end alignments. It compares the results of that with known results. Differences are printed to the console, and currently look something like this:

<pre>
diff phiX.sam.good phiX.sam
3c3
< @RG DT:2010-04-08T17:29Z ID:boingboing SM:NA12345
---
> @RG DT:2010-04-08T18:13Z ID:boingboing SM:NA12345
</pre>

Any differences greater than that are an error and need to be fixed by the author.

= Normal Workflow =

Karma works using a set of index and hash files created from an existing reference. Once created, this set of reference index and hash files must always be specified in the command line when aligning reads.

In concept, the simplest workflow is to first create a reference index using ''karma create'', then align reads using ''karma map''. You only have to build the index and hash once.

Because the reference can be large, and because Karma will share the reference among many running instances of Karma, it is useful to put well known references in a common location readily accessible to you and your collaborators.

= Build reference index and hash =

Building a reference index and hash with Karma is straightforward, but because it is time consuming for longer genomes, you typically save the reference index between runs.

The simplest example for creating a reference and index using a wordsize of 11-mer words is:

karma create -i -w 11 phiX.fa

More generally, three primary parameters are necessary for building a Karma reference index:

# a boolean flag indicating base or color space
# the index table word occurrence cutoff value
# the word size

Although the input reference is always expected to be base space and in [http://en.wikipedia.org/wiki/FASTA_format FASTA] format, the binary version of the reference, and the corresponding index and hash files, can be in either color space (ABI SOLiD) or base space (Illumina or LS454). For a given reference [http://en.wikipedia.org/wiki/FASTA_format FASTA] file, you may have either a color or base space binary reference, as well as either color or base space index/hash files, any in varying word sizes or occurrence cutoffs.

Because the index and hash files are dependent on the occurrence cutoff parameter and the word size, the output files created by karma have those values in the file name. This allows you to create a variety of index/hash tables, depending on your expected use (ABI SOLiD, in particular, is sensitive to read length).

== Options for building reference index and hash ==

-r ''reference'' Reference file in [http://en.wikipedia.org/wiki/FASTA_format FASTA] format
-w ''word size'' Word size for index and hash (default 15, typically 10-16)
-O ''occurrence cutoff'' Upper count of number of word positions to store in word positions table (default 5000)
-c Creates a color space reference and index/hash
-i Create the index and hash as well as the binary reference

= Aligning Reads =

Aligning reads to the reference is easy:

karma map -r phiX.fa -w 11 phiX.fastq

or for paired reads:

karma map -r phiX.fa -w 11 phiX-mate1.fastq phiX-mate2.fastq

In both of the above examples, the -r option names the reference originally used to build the index/hash, and the -w 11 specifies that we are using the index/hash built for 11-mer words. Although you can use the default word size of 15 for phiX, the index is 4^15 * 4 = 4GBytes, so a shorter word size is prudent.

Since Karma uses the word size and occurrence cutoff to help construct the actual index and hash filenames, you must specify them the same way you did when you created the reference index and hash.

== Aligning Reads (Illumina) ==

Karma is set up so that the default options work well for mapping Illumina reads to the Human genome.

== Aligning Reads (ABI SOLiD) ==

Karma has been designed to align color space reads. However, in Karma 0.9.0, this functionality is not working.

== Aligning Reads (LS 454) ==

Karma has been designed to align LS 454 reads. However, in Karma 0.9.0, this functionality is not working.

= Karma Performance Tuning =

There are four components to the Karma index and hash. A pure index array, based on an N-mer word index. This is used as a pointer into a word positions table, which is an ordered list of genome positions in which that N-mer word appears. There is a cap called the ''occurrence cutoff'', which once exceeded, causes that index word to be marked as a high repeat pattern. Once marked as high repeat, the N-mer word is instead combined with both the N-mer word preceding it, as well as the N-mer word succeeding it to create a 2 * N-mer word hash key. Two hash tables are populated, a left and a right hash. These are then used when that pattern is found in a read.

== Index Word Size ==

Choosing an appropriate word size for larger genome is critical to performance. The easiest case is for Illumina base space reads with the human genome (3Gbases), where the default 15-mer word size is fine.

For smaller genomes, consider using a smaller word size. Genomes smaller than a few million bases should be perfectly fine with a word size of 11 or 12.

Since the primary index table into the word positions table is 2^(wordsize) * 4 bytes, it can grow large rapidly. All else being equal, a smaller word size leads to longer sets of word positions for each index value. Each increment of word size approximately quadruples storage requirements, and halves runtime. Similarly, each decrement of word size reduces the index table size by 75%, and doubles runtime. These approximations are old, but serve a useful rule of thumb.

For ABI SOLiD reads, the word size is critical, due to the shorter length of reads as compared to Illumina or LS 454.

The optimal minimum word size is chosen such that it is 1/4 the minimum expected average read length. It also must be chosen to be 1/2 the minimum expected read length, since at least 2 full words must exist in the read.

So for 48-mer reads, a reasonable value of word size is 12. Although the base space default of 15 is fine, too, Karma is able to take advantage of a higher number of index words per read, yielding substantial speedups even with the shorter read. Similarly, 52-mer reads would map better with a 13-mer word size, and 56-mer reads would map best with a 14-mer word size.

== Occurrence Cutoff ==

The occurrence cutoff value determines how quickly an N-mer pattern is declared to be ''high repeat'' and left out of the index in favor of a hash. The default value of 5000 seems adequate for Illumina reads with the human genome. If ultimate performance is necessary, some experimentation is called for with this value.

== Shared Memory ==

Karma uses memory mapped files to share the potentially large reference index and hash data structures.

Karma uses this to great effect on our 8 processors with hyperthreading enabled. 16 copies of karma can share one reference index and hash, yielding a very acceptable memory per CPU ratio of around 1GB/CPU.

A problem with large reference index and hash data structures is that they are more prone to being paged out. On a shared machine that is being used extensively even just simple disk I/O, memory pages are being reclaimed such that Karma will become swapped out.

While Karma can recover on its own, it is best to either run in a production manner on dedicated machines, or to run a program such as the utility ''mapfile'' found in the utilities sub-folder. This program continually touches each page of the data structures in sequential order, forcing them to the head of the disk buffer pool, so they don't get aged out of the queue.

= Modifying the Reference Header =

''NB: This feature is not yet complete''

To facilitate SAM RG values being set automatically in a production environment, we keep a header in the binary version of the reference. The header can be viewed and edited using the header subcommands here.

To view the header:

karma header -r phiX.fa

To view and edit the header:

karma header -r phiX.fa -e

= Other test and check capabilities =

Due to the size and complexity of Karma input, output and index files, various checks and tests are useful, so we include some diagnostics capabilities:

Tests for external files:

karma check [options...] file.bam file.fastq file.sam file.fa file.umfa

Tests internal to Karma:
karma test [options...]
-d -> debug
-s [int] -> set random number seed [12345]

= Karma File structure =

Upon successfully building references, you will obtain a list of reference files like below:

{| style="width: 571px; height: 288px" border="1" cellspacing="1" cellpadding="1" width="571"
|-
|
|
Base Space

| Color Space
|-
|
Reference genome

|
NCBI37-bs.umfa

| NCBI37-cs.umfa
|-
|
Word Index

|
NCBI37-bs.15.5000.umwiwp

NCBI37-bs.15.5000.umwihi

|
NCBI37-cs.15.5000.umwiwp

NCBI37-cs.15.5000.umwihi

|-
|
Word Hash (Left)

|
NCBI37-bs.15.5000.umwhl

| NCBI37-cs.15.5000.umwhl
|-
|
Word Hash (Right)

|
NCBI37-bs.15.5000.umwhr

| NCBI37-cs.15.5000.umwhr
|}

 

 

= Karma TODO List =

# command line help is muddled up - UserOptions.h needs work
# color space read handling needs to be re-integrated and tested
# LS 454 needs to be tested, and the code adapted
# pre-process some number of records to establish an appropriate max insert size
# finish reference header view/edit code
# investigate and document maximum memory use during ''create'' sub-command
# finish and improve check and test commands

= Karma CHANGELOG =

* Karma 0.9.0
* Karma 0.8.8S
* Karma 0.8.8
* Karma 0.8.6

= Other useful links =

[http://lh3lh3.users.sourceforge.net/bioinfo.shtml Heng Li's thoughts about aligners]

[http://lh3lh3.users.sourceforge.net/udb.shtml Benchmark of Dictionary Structures]

Karma

2010-04-08T19:22:27Z

Pha: add features/benefits

[[Category:Software]]

'''K-tuple Alignment with Rapid Matching Algorithm'''

Karma uses an existing reference to align short reads, such as those generated by Illumina sequencers.

Primary features:

# High performance, high sensitivity
# Large and small gap detection by default
# Multiple gaps per read by default
# Single or paired end reads
# No read length limit
# Quality scores are used to assess quality of maps
# All potential locations are examined exhaustively, none are omitted
# Reasonable memory per CPU ratio on high core count machines

The current version, 0.9.0, is optimized to rapidly map base space reads from Illumina sequencers.

Color space and LS454 sequence alignments are not currently supported. These features will return in Karma 0.9.1.

= Download Karma =

To get a copy go to [http://www.sph.umich.edu/csg/pha/karma/download/ Karma Download]

= Build Karma =

Karma is designed to be reasonably portable.

However, since development occurs only on Ubuntu 9.10 x86 and x64 platforms, there are likely other portability issues.

We support Karma only on Ubuntu 9.10 on 64-bit processors.

== Dependencies ==

Karma requires that the following debian packages be installed on the host Linux machine:

# libssl-dev
# zlib1g-dev

Without these installed, Karma will not build.

== Building ==

Assuming the karma tar file is named karma.tgz, do the following

tar xvzf karma.tgz
cd karma-0.9
make
mkdir ~/bin
cp karma/karma ~/bin

Alternatively, if you want to share the karma binary install it in /usr/local/bin/karma.

== Testing the build ==

To test karma, go to the build tree subdirectory named ''karma'', and type the command:

make test

The test script builds a reference for the small phiX genome, then runs single end as well as paired end alignments. It compares the results of that with known results. Differences are printed to the console, and currently look something like this:

<pre>
diff phiX.sam.good phiX.sam
3c3
< @RG DT:2010-04-08T17:29Z ID:boingboing SM:NA12345
---
> @RG DT:2010-04-08T18:13Z ID:boingboing SM:NA12345
</pre>

Any differences greater than that are an error and need to be fixed by the author.

= Normal Workflow =

Karma works using a set of index and hash files created from an existing reference. Once created, this set of reference index and hash files must always be specified in the command line when aligning reads.

In concept, the simplest workflow is to first create a reference index using ''karma create'', then align reads using ''karma map''. You only have to build the index and hash once.

Because the reference can be large, and because Karma will share the reference among many running instances of Karma, it is useful to put well known references in a common location readily accessible to you and your collaborators.

= Build reference index and hash =

Building a reference index and hash with Karma is straightforward, but because it is time consuming for longer genomes, you typically save the reference index between runs.

The simplest example for creating a reference and index using a wordsize of 11-mer words is:

karma create -i -w 11 phiX.fa

More generally, three primary parameters are necessary for building a Karma reference index:

# a boolean flag indicating base or color space
# the index table word occurrence cutoff value
# the word size

Although the input reference is always expected to be base space and in [http://en.wikipedia.org/wiki/FASTA_format FASTA] format, the binary version of the reference, and the corresponding index and hash files, can be in either color space (ABI SOLiD) or base space (Illumina or LS454). For a given reference [http://en.wikipedia.org/wiki/FASTA_format FASTA] file, you may have either a color or base space binary reference, as well as either color or base space index/hash files, any in varying word sizes or occurrence cutoffs.

Because the index and hash files are dependent on the occurrence cutoff parameter and the word size, the output files created by karma have those values in the file name. This allows you to create a variety of index/hash tables, depending on your expected use (ABI SOLiD, in particular, is sensitive to read length).

== Options for building reference index and hash ==

-r ''reference'' Reference file in [http://en.wikipedia.org/wiki/FASTA_format FASTA] format
-w ''word size'' Word size for index and hash (default 15, typically 10-16)
-O ''occurrence cutoff'' Upper count of number of word positions to store in word positions table (default 5000)
-c Creates a color space reference and index/hash
-i Create the index and hash as well as the binary reference

= Aligning Reads =

Aligning reads to the reference is easy:

karma map -r phiX.fa -w 11 phiX.fastq

or for paired reads:

karma map -r phiX.fa -w 11 phiX-mate1.fastq phiX-mate2.fastq

In both of the above examples, the -r option names the reference originally used to build the index/hash, and the -w 11 specifies that we are using the index/hash built for 11-mer words. Although you can use the default word size of 15 for phiX, the index is 4^15 * 4 = 4GBytes, so a shorter word size is prudent.

Since Karma uses the word size and occurrence cutoff to help construct the actual index and hash filenames, you must specify them the same way you did when you created the reference index and hash.

== Aligning Reads (Illumina) ==

Karma is set up so that the default options work well for mapping Illumina reads to the Human genome.

== Aligning Reads (ABI SOLiD) ==

Karma has been designed to align color space reads. However, in Karma 0.9.0, this functionality is not working.

== Aligning Reads (LS 454) ==

Karma has been designed to align LS 454 reads. However, in Karma 0.9.0, this functionality is not working.

= Karma Performance Tuning =

There are four components to the Karma index and hash. A pure index array, based on an N-mer word index. This is used as a pointer into a word positions table, which is an ordered list of genome positions in which that N-mer word appears. There is a cap called the ''occurrence cutoff'', which once exceeded, causes that index word to be marked as a high repeat pattern. Once marked as high repeat, the N-mer word is instead combined with both the N-mer word preceding it, as well as the N-mer word succeeding it to create a 2 * N-mer word hash key. Two hash tables are populated, a left and a right hash. These are then used when that pattern is found in a read.

== Index Word Size ==

Choosing an appropriate word size for larger genome is critical to performance. The easiest case is for Illumina base space reads with the human genome (3Gbases), where the default 15-mer word size is fine.

For smaller genomes, consider using a smaller word size. Genomes smaller than a few million bases should be perfectly fine with a word size of 11 or 12.

Since the primary index table into the word positions table is 2^(wordsize) * 4 bytes, it can grow large rapidly. All else being equal, a smaller word size leads to longer sets of word positions for each index value. Each increment of word size approximately quadruples storage requirements, and halves runtime. Similarly, each decrement of word size reduces the index table size by 75%, and doubles runtime. These approximations are old, but serve a useful rule of thumb.

For ABI SOLiD reads, the word size is critical, due to the shorter length of reads as compared to Illumina or LS 454.

The optimal minimum word size is chosen such that it is 1/4 the minimum expected average read length. It also must be chosen to be 1/2 the minimum expected read length, since at least 2 full words must exist in the read.

So for 48-mer reads, a reasonable value of word size is 12. Although the base space default of 15 is fine, too, Karma is able to take advantage of a higher number of index words per read, yielding substantial speedups even with the shorter read. Similarly, 52-mer reads would map better with a 13-mer word size, and 56-mer reads would map best with a 14-mer word size.

== Occurrence Cutoff ==

The occurrence cutoff value determines how quickly an N-mer pattern is declared to be ''high repeat'' and left out of the index in favor of a hash. The default value of 5000 seems adequate for Illumina reads with the human genome. If ultimate performance is necessary, some experimentation is called for with this value.

== Shared Memory ==

Karma uses memory mapped files to share the potentially large reference index and hash data structures.

Karma uses this to great effect on our 8 processors with hyperthreading enabled. 16 copies of karma can share one reference index and hash, yielding a very acceptable memory per CPU ratio of around 1GB/CPU.

A problem with large reference index and hash data structures is that they are more prone to being paged out. On a shared machine that is being used extensively even just simple disk I/O, memory pages are being reclaimed such that Karma will become swapped out.

While Karma can recover on its own, it is best to either run in a production manner on dedicated machines, or to run a program such as the utility ''mapfile'' found in the utilities sub-folder. This program continually touches each page of the data structures in sequential order, forcing them to the head of the disk buffer pool, so they don't get aged out of the queue.

= Modifying the Reference Header =

''NB: This feature is not yet complete''

To facilitate SAM RG values being set automatically in a production environment, we keep a header in the binary version of the reference. The header can be viewed and edited using the header subcommands here.

To view the header:

karma header -r phiX.fa

To view and edit the header:

karma header -r phiX.fa -e

= Other test and check capabilities =

Due to the size and complexity of Karma input, output and index files, various checks and tests are useful, so we include some diagnostics capabilities:

Tests for external files:

karma check [options...] file.bam file.fastq file.sam file.fa file.umfa

Tests internal to Karma:
karma test [options...]
-d -> debug
-s [int] -> set random number seed [12345]

= Karma File structure =

Upon successfully building references, you will obtain a list of reference files like below:

{| style="width: 571px; height: 288px" border="1" cellspacing="1" cellpadding="1" width="571"
|-
|
|
Base Space

| Color Space
|-
|
Reference genome

|
NCBI37-bs.umfa

| NCBI37-cs.umfa
|-
|
Word Index

|
NCBI37-bs.15.5000.umwiwp

NCBI37-bs.15.5000.umwihi

|
NCBI37-cs.15.5000.umwiwp

NCBI37-cs.15.5000.umwihi

|-
|
Word Hash (Left)

|
NCBI37-bs.15.5000.umwhl

| NCBI37-cs.15.5000.umwhl
|-
|
Word Hash (Right)

|
NCBI37-bs.15.5000.umwhr

| NCBI37-cs.15.5000.umwhr
|}

 

 

= Karma TODO List =

= Karma CHANGELOG =

= Other useful links =

[http://lh3lh3.users.sourceforge.net/bioinfo.shtml Heng Li's thoughts about aligners]

[http://lh3lh3.users.sourceforge.net/udb.shtml Benchmark of Dictionary Structures]

Karma

2010-04-08T19:13:57Z

Pha: /* Other useful links */

[[Category:Software]]

'''K-tuple Alignment with Rapid Matching Algorithm'''

Karma uses an existing reference to align short reads, such as those generated by Illumina sequencers.

The current version, 0.9.0, is optimized to rapidly map base space reads from Illumina sequencers.

Color space and LS454 sequence alignments are not currently supported. These features will return in Karma 0.9.1.

= Download Karma =

To get a copy go to [http://www.sph.umich.edu/csg/pha/karma/download/ Karma Download]

= Build Karma =

Karma is designed to be reasonably portable.

However, since development occurs only on Ubuntu 9.10 x86 and x64 platforms, there are likely other portability issues.

We support Karma only on Ubuntu 9.10 on 64-bit processors.

== Dependencies ==

Karma requires that the following debian packages be installed on the host Linux machine:

# libssl-dev
# zlib1g-dev

Without these installed, Karma will not build.

== Building ==

Assuming the karma tar file is named karma.tgz, do the following

tar xvzf karma.tgz
cd karma-0.9
make
mkdir ~/bin
cp karma/karma ~/bin

Alternatively, if you want to share the karma binary install it in /usr/local/bin/karma.

== Testing the build ==

To test karma, go to the build tree subdirectory named ''karma'', and type the command:

make test

The test script builds a reference for the small phiX genome, then runs single end as well as paired end alignments. It compares the results of that with known results. Differences are printed to the console, and currently look something like this:

<pre>
diff phiX.sam.good phiX.sam
3c3
< @RG DT:2010-04-08T17:29Z ID:boingboing SM:NA12345
---
> @RG DT:2010-04-08T18:13Z ID:boingboing SM:NA12345
</pre>

Any differences greater than that are an error and need to be fixed by the author.

= Normal Workflow =

Karma works using a set of index and hash files created from an existing reference. Once created, this set of reference index and hash files must always be specified in the command line when aligning reads.

In concept, the simplest workflow is to first create a reference index using ''karma create'', then align reads using ''karma map''. You only have to build the index and hash once.

Because the reference can be large, and because Karma will share the reference among many running instances of Karma, it is useful to put well known references in a common location readily accessible to you and your collaborators.

= Build reference index and hash =

Building a reference index and hash with Karma is straightforward, but because it is time consuming for longer genomes, you typically save the reference index between runs.

The simplest example for creating a reference and index using a wordsize of 11-mer words is:

karma create -i -w 11 phiX.fa

More generally, three primary parameters are necessary for building a Karma reference index:

# a boolean flag indicating base or color space
# the index table word occurrence cutoff value
# the word size

Although the input reference is always expected to be base space and in [http://en.wikipedia.org/wiki/FASTA_format FASTA] format, the binary version of the reference, and the corresponding index and hash files, can be in either color space (ABI SOLiD) or base space (Illumina or LS454). For a given reference [http://en.wikipedia.org/wiki/FASTA_format FASTA] file, you may have either a color or base space binary reference, as well as either color or base space index/hash files, any in varying word sizes or occurrence cutoffs.

Because the index and hash files are dependent on the occurrence cutoff parameter and the word size, the output files created by karma have those values in the file name. This allows you to create a variety of index/hash tables, depending on your expected use (ABI SOLiD, in particular, is sensitive to read length).

== Options for building reference index and hash ==

-r ''reference'' Reference file in [http://en.wikipedia.org/wiki/FASTA_format FASTA] format
-w ''word size'' Word size for index and hash (default 15, typically 10-16)
-O ''occurrence cutoff'' Upper count of number of word positions to store in word positions table (default 5000)
-c Creates a color space reference and index/hash
-i Create the index and hash as well as the binary reference

= Aligning Reads =

Aligning reads to the reference is easy:

karma map -r phiX.fa -w 11 phiX.fastq

or for paired reads:

karma map -r phiX.fa -w 11 phiX-mate1.fastq phiX-mate2.fastq

In both of the above examples, the -r option names the reference originally used to build the index/hash, and the -w 11 specifies that we are using the index/hash built for 11-mer words. Although you can use the default word size of 15 for phiX, the index is 4^15 * 4 = 4GBytes, so a shorter word size is prudent.

Since Karma uses the word size and occurrence cutoff to help construct the actual index and hash filenames, you must specify them the same way you did when you created the reference index and hash.

== Aligning Reads (Illumina) ==

Karma is set up so that the default options work well for mapping Illumina reads to the Human genome.

== Aligning Reads (ABI SOLiD) ==

Karma has been designed to align color space reads. However, in Karma 0.9.0, this functionality is not working.

== Aligning Reads (LS 454) ==

Karma has been designed to align LS 454 reads. However, in Karma 0.9.0, this functionality is not working.

= Karma Performance Tuning =

There are four components to the Karma index and hash. A pure index array, based on an N-mer word index. This is used as a pointer into a word positions table, which is an ordered list of genome positions in which that N-mer word appears. There is a cap called the ''occurrence cutoff'', which once exceeded, causes that index word to be marked as a high repeat pattern. Once marked as high repeat, the N-mer word is instead combined with both the N-mer word preceding it, as well as the N-mer word succeeding it to create a 2 * N-mer word hash key. Two hash tables are populated, a left and a right hash. These are then used when that pattern is found in a read.

== Index Word Size ==

Choosing an appropriate word size for larger genome is critical to performance. The easiest case is for Illumina base space reads with the human genome (3Gbases), where the default 15-mer word size is fine.

For smaller genomes, consider using a smaller word size. Genomes smaller than a few million bases should be perfectly fine with a word size of 11 or 12.

Since the primary index table into the word positions table is 2^(wordsize) * 4 bytes, it can grow large rapidly. All else being equal, a smaller word size leads to longer sets of word positions for each index value. Each increment of word size approximately quadruples storage requirements, and halves runtime. Similarly, each decrement of word size reduces the index table size by 75%, and doubles runtime. These approximations are old, but serve a useful rule of thumb.

For ABI SOLiD reads, the word size is critical, due to the shorter length of reads as compared to Illumina or LS 454.

The optimal minimum word size is chosen such that it is 1/4 the minimum expected average read length. It also must be chosen to be 1/2 the minimum expected read length, since at least 2 full words must exist in the read.

So for 48-mer reads, a reasonable value of word size is 12. Although the base space default of 15 is fine, too, Karma is able to take advantage of a higher number of index words per read, yielding substantial speedups even with the shorter read. Similarly, 52-mer reads would map better with a 13-mer word size, and 56-mer reads would map best with a 14-mer word size.

== Occurrence Cutoff ==

The occurrence cutoff value determines how quickly an N-mer pattern is declared to be ''high repeat'' and left out of the index in favor of a hash. The default value of 5000 seems adequate for Illumina reads with the human genome. If ultimate performance is necessary, some experimentation is called for with this value.

== Shared Memory ==

Karma uses memory mapped files to share the potentially large reference index and hash data structures.

Karma uses this to great effect on our 8 processors with hyperthreading enabled. 16 copies of karma can share one reference index and hash, yielding a very acceptable memory per CPU ratio of around 1GB/CPU.

A problem with large reference index and hash data structures is that they are more prone to being paged out. On a shared machine that is being used extensively even just simple disk I/O, memory pages are being reclaimed such that Karma will become swapped out.

While Karma can recover on its own, it is best to either run in a production manner on dedicated machines, or to run a program such as the utility ''mapfile'' found in the utilities sub-folder. This program continually touches each page of the data structures in sequential order, forcing them to the head of the disk buffer pool, so they don't get aged out of the queue.

= Modifying the Reference Header =

''NB: This feature is not yet complete''

To facilitate SAM RG values being set automatically in a production environment, we keep a header in the binary version of the reference. The header can be viewed and edited using the header subcommands here.

To view the header:

karma header -r phiX.fa

To view and edit the header:

karma header -r phiX.fa -e

= Other test and check capabilities =

Due to the size and complexity of Karma input, output and index files, various checks and tests are useful, so we include some diagnostics capabilities:

Tests for external files:

karma check [options...] file.bam file.fastq file.sam file.fa file.umfa

Tests internal to Karma:
karma test [options...]
-d -> debug
-s [int] -> set random number seed [12345]

= Karma File structure =

Upon successfully building references, you will obtain a list of reference files like below:

{| style="width: 571px; height: 288px" border="1" cellspacing="1" cellpadding="1" width="571"
|-
|
|
Base Space

| Color Space
|-
|
Reference genome

|
NCBI37-bs.umfa

| NCBI37-cs.umfa
|-
|
Word Index

|
NCBI37-bs.15.5000.umwiwp

NCBI37-bs.15.5000.umwihi

|
NCBI37-cs.15.5000.umwiwp

NCBI37-cs.15.5000.umwihi

|-
|
Word Hash (Left)

|
NCBI37-bs.15.5000.umwhl

| NCBI37-cs.15.5000.umwhl
|-
|
Word Hash (Right)

|
NCBI37-bs.15.5000.umwhr

| NCBI37-cs.15.5000.umwhr
|}

 

 

= Karma TODO List =

= Karma CHANGELOG =

= Other useful links =

[http://lh3lh3.users.sourceforge.net/bioinfo.shtml Heng Li's thoughts about aligners]

[http://lh3lh3.users.sourceforge.net/udb.shtml Benchmark of Dictionary Structures]

Karma

2010-04-08T19:13:01Z

Pha: /* Karma Performance Tuning */

[[Category:Software]]

'''K-tuple Alignment with Rapid Matching Algorithm'''

Karma uses an existing reference to align short reads, such as those generated by Illumina sequencers.

The current version, 0.9.0, is optimized to rapidly map base space reads from Illumina sequencers.

Color space and LS454 sequence alignments are not currently supported. These features will return in Karma 0.9.1.

= Download Karma =

To get a copy go to [http://www.sph.umich.edu/csg/pha/karma/download/ Karma Download]

= Build Karma =

Karma is designed to be reasonably portable.

However, since development occurs only on Ubuntu 9.10 x86 and x64 platforms, there are likely other portability issues.

We support Karma only on Ubuntu 9.10 on 64-bit processors.

== Dependencies ==

Karma requires that the following debian packages be installed on the host Linux machine:

# libssl-dev
# zlib1g-dev

Without these installed, Karma will not build.

== Building ==

Assuming the karma tar file is named karma.tgz, do the following

tar xvzf karma.tgz
cd karma-0.9
make
mkdir ~/bin
cp karma/karma ~/bin

Alternatively, if you want to share the karma binary install it in /usr/local/bin/karma.

== Testing the build ==

To test karma, go to the build tree subdirectory named ''karma'', and type the command:

make test

The test script builds a reference for the small phiX genome, then runs single end as well as paired end alignments. It compares the results of that with known results. Differences are printed to the console, and currently look something like this:

<pre>
diff phiX.sam.good phiX.sam
3c3
< @RG DT:2010-04-08T17:29Z ID:boingboing SM:NA12345
---
> @RG DT:2010-04-08T18:13Z ID:boingboing SM:NA12345
</pre>

Any differences greater than that are an error and need to be fixed by the author.

= Normal Workflow =

Karma works using a set of index and hash files created from an existing reference. Once created, this set of reference index and hash files must always be specified in the command line when aligning reads.

In concept, the simplest workflow is to first create a reference index using ''karma create'', then align reads using ''karma map''. You only have to build the index and hash once.

Because the reference can be large, and because Karma will share the reference among many running instances of Karma, it is useful to put well known references in a common location readily accessible to you and your collaborators.

= Build reference index and hash =

Building a reference index and hash with Karma is straightforward, but because it is time consuming for longer genomes, you typically save the reference index between runs.

The simplest example for creating a reference and index using a wordsize of 11-mer words is:

karma create -i -w 11 phiX.fa

More generally, three primary parameters are necessary for building a Karma reference index:

# a boolean flag indicating base or color space
# the index table word occurrence cutoff value
# the word size

Although the input reference is always expected to be base space and in [http://en.wikipedia.org/wiki/FASTA_format FASTA] format, the binary version of the reference, and the corresponding index and hash files, can be in either color space (ABI SOLiD) or base space (Illumina or LS454). For a given reference [http://en.wikipedia.org/wiki/FASTA_format FASTA] file, you may have either a color or base space binary reference, as well as either color or base space index/hash files, any in varying word sizes or occurrence cutoffs.

Because the index and hash files are dependent on the occurrence cutoff parameter and the word size, the output files created by karma have those values in the file name. This allows you to create a variety of index/hash tables, depending on your expected use (ABI SOLiD, in particular, is sensitive to read length).

== Options for building reference index and hash ==

-r ''reference'' Reference file in [http://en.wikipedia.org/wiki/FASTA_format FASTA] format
-w ''word size'' Word size for index and hash (default 15, typically 10-16)
-O ''occurrence cutoff'' Upper count of number of word positions to store in word positions table (default 5000)
-c Creates a color space reference and index/hash
-i Create the index and hash as well as the binary reference

= Aligning Reads =

Aligning reads to the reference is easy:

karma map -r phiX.fa -w 11 phiX.fastq

or for paired reads:

karma map -r phiX.fa -w 11 phiX-mate1.fastq phiX-mate2.fastq

In both of the above examples, the -r option names the reference originally used to build the index/hash, and the -w 11 specifies that we are using the index/hash built for 11-mer words. Although you can use the default word size of 15 for phiX, the index is 4^15 * 4 = 4GBytes, so a shorter word size is prudent.

Since Karma uses the word size and occurrence cutoff to help construct the actual index and hash filenames, you must specify them the same way you did when you created the reference index and hash.

== Aligning Reads (Illumina) ==

Karma is set up so that the default options work well for mapping Illumina reads to the Human genome.

== Aligning Reads (ABI SOLiD) ==

Karma has been designed to align color space reads. However, in Karma 0.9.0, this functionality is not working.

== Aligning Reads (LS 454) ==

Karma has been designed to align LS 454 reads. However, in Karma 0.9.0, this functionality is not working.

= Karma Performance Tuning =

There are four components to the Karma index and hash. A pure index array, based on an N-mer word index. This is used as a pointer into a word positions table, which is an ordered list of genome positions in which that N-mer word appears. There is a cap called the ''occurrence cutoff'', which once exceeded, causes that index word to be marked as a high repeat pattern. Once marked as high repeat, the N-mer word is instead combined with both the N-mer word preceding it, as well as the N-mer word succeeding it to create a 2 * N-mer word hash key. Two hash tables are populated, a left and a right hash. These are then used when that pattern is found in a read.

== Index Word Size ==

Choosing an appropriate word size for larger genome is critical to performance. The easiest case is for Illumina base space reads with the human genome (3Gbases), where the default 15-mer word size is fine.

For smaller genomes, consider using a smaller word size. Genomes smaller than a few million bases should be perfectly fine with a word size of 11 or 12.

Since the primary index table into the word positions table is 2^(wordsize) * 4 bytes, it can grow large rapidly. All else being equal, a smaller word size leads to longer sets of word positions for each index value. Each increment of word size approximately quadruples storage requirements, and halves runtime. Similarly, each decrement of word size reduces the index table size by 75%, and doubles runtime. These approximations are old, but serve a useful rule of thumb.

For ABI SOLiD reads, the word size is critical, due to the shorter length of reads as compared to Illumina or LS 454.

The optimal minimum word size is chosen such that it is 1/4 the minimum expected average read length. It also must be chosen to be 1/2 the minimum expected read length, since at least 2 full words must exist in the read.

So for 48-mer reads, a reasonable value of word size is 12. Although the base space default of 15 is fine, too, Karma is able to take advantage of a higher number of index words per read, yielding substantial speedups even with the shorter read. Similarly, 52-mer reads would map better with a 13-mer word size, and 56-mer reads would map best with a 14-mer word size.

== Occurrence Cutoff ==

The occurrence cutoff value determines how quickly an N-mer pattern is declared to be ''high repeat'' and left out of the index in favor of a hash. The default value of 5000 seems adequate for Illumina reads with the human genome. If ultimate performance is necessary, some experimentation is called for with this value.

== Shared Memory ==

Karma uses memory mapped files to share the potentially large reference index and hash data structures.

Karma uses this to great effect on our 8 processors with hyperthreading enabled. 16 copies of karma can share one reference index and hash, yielding a very acceptable memory per CPU ratio of around 1GB/CPU.

A problem with large reference index and hash data structures is that they are more prone to being paged out. On a shared machine that is being used extensively even just simple disk I/O, memory pages are being reclaimed such that Karma will become swapped out.

While Karma can recover on its own, it is best to either run in a production manner on dedicated machines, or to run a program such as the utility ''mapfile'' found in the utilities sub-folder. This program continually touches each page of the data structures in sequential order, forcing them to the head of the disk buffer pool, so they don't get aged out of the queue.

= Modifying the Reference Header =

''NB: This feature is not yet complete''

To facilitate SAM RG values being set automatically in a production environment, we keep a header in the binary version of the reference. The header can be viewed and edited using the header subcommands here.

To view the header:

karma header -r phiX.fa

To view and edit the header:

karma header -r phiX.fa -e

= Other test and check capabilities =

Due to the size and complexity of Karma input, output and index files, various checks and tests are useful, so we include some diagnostics capabilities:

Tests for external files:

karma check [options...] file.bam file.fastq file.sam file.fa file.umfa

Tests internal to Karma:
karma test [options...]
-d -> debug
-s [int] -> set random number seed [12345]

= Karma File structure =

Upon successfully building references, you will obtain a list of reference files like below:

{| style="width: 571px; height: 288px" border="1" cellspacing="1" cellpadding="1" width="571"
|-
|
|
Base Space

| Color Space
|-
|
Reference genome

|
NCBI37-bs.umfa

| NCBI37-cs.umfa
|-
|
Word Index

|
NCBI37-bs.15.5000.umwiwp

NCBI37-bs.15.5000.umwihi

|
NCBI37-cs.15.5000.umwiwp

NCBI37-cs.15.5000.umwihi

|-
|
Word Hash (Left)

|
NCBI37-bs.15.5000.umwhl

| NCBI37-cs.15.5000.umwhl
|-
|
Word Hash (Right)

|
NCBI37-bs.15.5000.umwhr

| NCBI37-cs.15.5000.umwhr
|}

 

 

= Karma TODO List =

= Karma CHANGELOG =

= Other useful links =

[http://lh3lh3.users.sourceforge.net/bioinfo.shtml Heng Li's thoughts about aligner]

[http://lh3lh3.users.sourceforge.net/udb.shtml Benchmark of Dictionary Structures]

Karma

2010-04-08T19:11:14Z

Pha: fill in build info

[[Category:Software]]

'''K-tuple Alignment with Rapid Matching Algorithm'''

Karma uses an existing reference to align short reads, such as those generated by Illumina sequencers.

The current version, 0.9.0, is optimized to rapidly map base space reads from Illumina sequencers.

Color space and LS454 sequence alignments are not currently supported. These features will return in Karma 0.9.1.

= Download Karma =

To get a copy go to [http://www.sph.umich.edu/csg/pha/karma/download/ Karma Download]

= Build Karma =

Karma is designed to be reasonably portable.

However, since development occurs only on Ubuntu 9.10 x86 and x64 platforms, there are likely other portability issues.

We support Karma only on Ubuntu 9.10 on 64-bit processors.

== Dependencies ==

Karma requires that the following debian packages be installed on the host Linux machine:

# libssl-dev
# zlib1g-dev

Without these installed, Karma will not build.

== Building ==

Assuming the karma tar file is named karma.tgz, do the following

tar xvzf karma.tgz
cd karma-0.9
make
mkdir ~/bin
cp karma/karma ~/bin

Alternatively, if you want to share the karma binary install it in /usr/local/bin/karma.

== Testing the build ==

To test karma, go to the build tree subdirectory named ''karma'', and type the command:

make test

The test script builds a reference for the small phiX genome, then runs single end as well as paired end alignments. It compares the results of that with known results. Differences are printed to the console, and currently look something like this:

<pre>
diff phiX.sam.good phiX.sam
3c3
< @RG DT:2010-04-08T17:29Z ID:boingboing SM:NA12345
---
> @RG DT:2010-04-08T18:13Z ID:boingboing SM:NA12345
</pre>

Any differences greater than that are an error and need to be fixed by the author.

= Normal Workflow =

Karma works using a set of index and hash files created from an existing reference. Once created, this set of reference index and hash files must always be specified in the command line when aligning reads.

In concept, the simplest workflow is to first create a reference index using ''karma create'', then align reads using ''karma map''. You only have to build the index and hash once.

Because the reference can be large, and because Karma will share the reference among many running instances of Karma, it is useful to put well known references in a common location readily accessible to you and your collaborators.

= Build reference index and hash =

Building a reference index and hash with Karma is straightforward, but because it is time consuming for longer genomes, you typically save the reference index between runs.

The simplest example for creating a reference and index using a wordsize of 11-mer words is:

karma create -i -w 11 phiX.fa

More generally, three primary parameters are necessary for building a Karma reference index:

# a boolean flag indicating base or color space
# the index table word occurrence cutoff value
# the word size

Although the input reference is always expected to be base space and in [http://en.wikipedia.org/wiki/FASTA_format FASTA] format, the binary version of the reference, and the corresponding index and hash files, can be in either color space (ABI SOLiD) or base space (Illumina or LS454). For a given reference [http://en.wikipedia.org/wiki/FASTA_format FASTA] file, you may have either a color or base space binary reference, as well as either color or base space index/hash files, any in varying word sizes or occurrence cutoffs.

Because the index and hash files are dependent on the occurrence cutoff parameter and the word size, the output files created by karma have those values in the file name. This allows you to create a variety of index/hash tables, depending on your expected use (ABI SOLiD, in particular, is sensitive to read length).

== Options for building reference index and hash ==

-r ''reference'' Reference file in [http://en.wikipedia.org/wiki/FASTA_format FASTA] format
-w ''word size'' Word size for index and hash (default 15, typically 10-16)
-O ''occurrence cutoff'' Upper count of number of word positions to store in word positions table (default 5000)
-c Creates a color space reference and index/hash
-i Create the index and hash as well as the binary reference

= Aligning Reads =

Aligning reads to the reference is easy:

karma map -r phiX.fa -w 11 phiX.fastq

or for paired reads:

karma map -r phiX.fa -w 11 phiX-mate1.fastq phiX-mate2.fastq

In both of the above examples, the -r option names the reference originally used to build the index/hash, and the -w 11 specifies that we are using the index/hash built for 11-mer words. Although you can use the default word size of 15 for phiX, the index is 4^15 * 4 = 4GBytes, so a shorter word size is prudent.

Since Karma uses the word size and occurrence cutoff to help construct the actual index and hash filenames, you must specify them the same way you did when you created the reference index and hash.

== Aligning Reads (Illumina) ==

Karma is set up so that the default options work well for mapping Illumina reads to the Human genome.

== Aligning Reads (ABI SOLiD) ==

Karma has been designed to align color space reads. However, in Karma 0.9.0, this functionality is not working.

== Aligning Reads (LS 454) ==

Karma has been designed to align LS 454 reads. However, in Karma 0.9.0, this functionality is not working.

= Karma Performance Tuning =

There are four components to the Karma index and hash. A pure index, based on an N-mer word index. This is used as a pointer into a word positions table, which is an ordered list of genome positions in which that N-mer word appears. There is a cap called the ''occurrence cutoff'', which once crossed, causes that index word to be marked as a high repeat pattern. Once marked as high repeat, the N-mer word is now combined with both the N-mer word preceding it, as well as the N-mer word succeeding it to create a 2 * N-mer word hash key. Two hash tables are populated, a left and a right hash.

== Index Word Size ==

Choosing an appropriate word size for larger genome is critical to performance. The easiest case is for Illumina base space reads with the human genome (3Gbases), where the default 15-mer word size is fine.

For smaller genomes, consider using a smaller word size. Genomes smaller than a few million bases should be perfectly fine with a word size of 11 or 12.

Since the primary index table into the word positions table is 2^(wordsize) * 4 bytes, it can grow large rapidly. All else being equal, a smaller word size leads to longer sets of word positions for each index value. Each increment of word size approximately quadruples storage requirements, and halves runtime. Similarly, each decrement of word size reduces the index table size by 75%, and doubles runtime. These approximations are old, but serve a useful rule of thumb.

For ABI SOLiD reads, the word size is critical, due to the shorter length of reads as compared to Illumina or LS 454.

The optimal minimum word size is chosen such that it is 1/4 the minimum expected average read length. It also must be chosen to be 1/2 the minimum expected read length, since at least 2 full words must exist in the read.

So for 48-mer reads, a reasonable value of word size is 12. Although the base space default of 15 is fine, too, Karma is able to take advantage of a higher number of index words per read, yielding substantial speedups even with the shorter read. Similarly, 52-mer reads would map better with a 13-mer word size, and 56-mer reads would map best with a 14-mer word size.

== Occurrence Cutoff ==

The occurrence cutoff value determines how quickly an N-mer pattern is declared to be ''high repeat'' and left out of the index in favor of a hash. The default value of 5000 seems adequate for Illumina reads with the human genome. If ultimate performance is necessary, some experimentation is called for with this value.

== Shared Memory ==

Karma uses memory mapped files to share the potentially large reference index and hash data structures.

Karma uses this to great effect on our 8 processors with hyperthreading enabled. 16 copies of karma can share one reference index and hash, yielding a very acceptable memory per CPU ratio of around 1GB/CPU.

A problem with large reference index and hash data structures is that they are more prone to being paged out. On a shared machine that is being used extensively even just simple disk I/O, memory pages are being reclaimed such that Karma will become swapped out.

While Karma can recover on its own, it is best to either run in a production manner on dedicated machines, or to run a program such as the utility ''mapfile'' found in the utilities sub-folder. This program continually touches each page of the data structures in sequential order, forcing them to the head of the disk buffer pool, so they don't get aged out of the queue.

= Modifying the Reference Header =

''NB: This feature is not yet complete''

To facilitate SAM RG values being set automatically in a production environment, we keep a header in the binary version of the reference. The header can be viewed and edited using the header subcommands here.

To view the header:

karma header -r phiX.fa

To view and edit the header:

karma header -r phiX.fa -e

= Other test and check capabilities =

Due to the size and complexity of Karma input, output and index files, various checks and tests are useful, so we include some diagnostics capabilities:

Tests for external files:

karma check [options...] file.bam file.fastq file.sam file.fa file.umfa

Tests internal to Karma:
karma test [options...]
-d -> debug
-s [int] -> set random number seed [12345]

= Karma File structure =

Upon successfully building references, you will obtain a list of reference files like below:

{| style="width: 571px; height: 288px" border="1" cellspacing="1" cellpadding="1" width="571"
|-
|
|
Base Space

| Color Space
|-
|
Reference genome

|
NCBI37-bs.umfa

| NCBI37-cs.umfa
|-
|
Word Index

|
NCBI37-bs.15.5000.umwiwp

NCBI37-bs.15.5000.umwihi

|
NCBI37-cs.15.5000.umwiwp

NCBI37-cs.15.5000.umwihi

|-
|
Word Hash (Left)

|
NCBI37-bs.15.5000.umwhl

| NCBI37-cs.15.5000.umwhl
|-
|
Word Hash (Right)

|
NCBI37-bs.15.5000.umwhr

| NCBI37-cs.15.5000.umwhr
|}

 

 

= Karma TODO List =

= Karma CHANGELOG =

= Other useful links =

[http://lh3lh3.users.sourceforge.net/bioinfo.shtml Heng Li's thoughts about aligner]

[http://lh3lh3.users.sourceforge.net/udb.shtml Benchmark of Dictionary Structures]

Karma

2010-04-08T19:05:38Z

Pha: /* Dependencies */ fill em in

[[Category:Software]]

'''K-tuple Alignment with Rapid Matching Algorithm'''

Karma uses an existing reference to align short reads, such as generated by Illumina sequencers.

The current version, 0.9.0, is optimized to rapidly map base space reads from Illumina sequencers.

Color space and LS454 sequence alignments are not working. These features will return in Karma 0.9.1.

= Download Karma =

To get a copy go to [http://www.sph.umich.edu/csg/pha/karma/download/ Karma Download]

= Build Karma =

== Dependencies ==

Karma requires that the following debian packages be installed on the host Linux machine:

# libssl-dev
# zlib1g-dev

Without these installed, Karma will not build.

== Building ==

== Testing the build ==

To test karma, go to the build tree subdirectory named ''karma'', and type the command:

make test

The test script builds a reference for the small phiX genome, then runs single end as well as paired end alignments. It compares the results of that with known results. Differences are printed to the console, and currently look something like this:

<pre>
diff phiX.sam.good phiX.sam
3c3
< @RG DT:2010-04-08T17:29Z ID:boingboing SM:NA12345
---
> @RG DT:2010-04-08T18:13Z ID:boingboing SM:NA12345
</pre>

Any differences greater than that are an error and need to be fixed by the author.

= Normal Workflow =

Karma works using a set of index and hash files created from an existing reference. Once created, this set of reference index and hash files must always be specified in the command line when aligning reads.

In concept, the simplest workflow is to first create a reference index using ''karma create'', then align reads using ''karma map''. You only have to build the index and hash once.

Because the reference can be large, and because Karma will share the reference among many running instances of Karma, it is useful to put well known references in a common location readily accessible to you and your collaborators.

= Build reference index and hash =

Building a reference index and hash with Karma is straightforward, but because it is time consuming for longer genomes, you typically save the reference index between runs.

The simplest example for creating a reference and index using a wordsize of 11-mer words is:

karma create -i -w 11 phiX.fa

More generally, three primary parameters are necessary for building a Karma reference index:

# a boolean flag indicating base or color space
# the index table word occurrence cutoff value
# the word size

Although the input reference is always expected to be base space and in [http://en.wikipedia.org/wiki/FASTA_format FASTA] format, the binary version of the reference, and the corresponding index and hash files, can be in either color space (ABI SOLiD) or base space (Illumina or LS454). For a given reference [http://en.wikipedia.org/wiki/FASTA_format FASTA] file, you may have either a color or base space binary reference, as well as either color or base space index/hash files, any in varying word sizes or occurrence cutoffs.

Because the index and hash files are dependent on the occurrence cutoff parameter and the word size, the output files created by karma have those values in the file name. This allows you to create a variety of index/hash tables, depending on your expected use (ABI SOLiD, in particular, is sensitive to read length).

== Options for building reference index and hash ==

-r ''reference'' Reference file in [http://en.wikipedia.org/wiki/FASTA_format FASTA] format
-w ''word size'' Word size for index and hash (default 15, typically 10-16)
-O ''occurrence cutoff'' Upper count of number of word positions to store in word positions table (default 5000)
-c Creates a color space reference and index/hash
-i Create the index and hash as well as the binary reference

= Aligning Reads =

Aligning reads to the reference is easy:

karma map -r phiX.fa -w 11 phiX.fastq

or for paired reads:

karma map -r phiX.fa -w 11 phiX-mate1.fastq phiX-mate2.fastq

In both of the above examples, the -r option names the reference originally used to build the index/hash, and the -w 11 specifies that we are using the index/hash built for 11-mer words. Although you can use the default word size of 15 for phiX, the index is 4^15 * 4 = 4GBytes, so a shorter word size is prudent.

Since Karma uses the word size and occurrence cutoff to help construct the actual index and hash filenames, you must specify them the same way you did when you created the reference index and hash.

== Aligning Reads (Illumina) ==

Karma is set up so that the default options work well for mapping Illumina reads to the Human genome.

== Aligning Reads (ABI SOLiD) ==

Karma has been designed to align color space reads. However, in Karma 0.9.0, this functionality is not working.

== Aligning Reads (LS 454) ==

Karma has been designed to align LS 454 reads. However, in Karma 0.9.0, this functionality is not working.

= Karma Performance Tuning =

There are four components to the Karma index and hash. A pure index, based on an N-mer word index. This is used as a pointer into a word positions table, which is an ordered list of genome positions in which that N-mer word appears. There is a cap called the ''occurrence cutoff'', which once crossed, causes that index word to be marked as a high repeat pattern. Once marked as high repeat, the N-mer word is now combined with both the N-mer word preceding it, as well as the N-mer word succeeding it to create a 2 * N-mer word hash key. Two hash tables are populated, a left and a right hash.

== Index Word Size ==

Choosing an appropriate word size for larger genome is critical to performance. The easiest case is for Illumina base space reads with the human genome (3Gbases), where the default 15-mer word size is fine.

For smaller genomes, consider using a smaller word size. Genomes smaller than a few million bases should be perfectly fine with a word size of 11 or 12.

Since the primary index table into the word positions table is 2^(wordsize) * 4 bytes, it can grow large rapidly. All else being equal, a smaller word size leads to longer sets of word positions for each index value. Each increment of word size approximately quadruples storage requirements, and halves runtime. Similarly, each decrement of word size reduces the index table size by 75%, and doubles runtime. These approximations are old, but serve a useful rule of thumb.

For ABI SOLiD reads, the word size is critical, due to the shorter length of reads as compared to Illumina or LS 454.

The optimal minimum word size is chosen such that it is 1/4 the minimum expected average read length. It also must be chosen to be 1/2 the minimum expected read length, since at least 2 full words must exist in the read.

So for 48-mer reads, a reasonable value of word size is 12. Although the base space default of 15 is fine, too, Karma is able to take advantage of a higher number of index words per read, yielding substantial speedups even with the shorter read. Similarly, 52-mer reads would map better with a 13-mer word size, and 56-mer reads would map best with a 14-mer word size.

== Occurrence Cutoff ==

The occurrence cutoff value determines how quickly an N-mer pattern is declared to be ''high repeat'' and left out of the index in favor of a hash. The default value of 5000 seems adequate for Illumina reads with the human genome. If ultimate performance is necessary, some experimentation is called for with this value.

== Shared Memory ==

Karma uses memory mapped files to share the potentially large reference index and hash data structures.

Karma uses this to great effect on our 8 processors with hyperthreading enabled. 16 copies of karma can share one reference index and hash, yielding a very acceptable memory per CPU ratio of around 1GB/CPU.

A problem with large reference index and hash data structures is that they are more prone to being paged out. On a shared machine that is being used extensively even just simple disk I/O, memory pages are being reclaimed such that Karma will become swapped out.

While Karma can recover on its own, it is best to either run in a production manner on dedicated machines, or to run a program such as the utility ''mapfile'' found in the utilities sub-folder. This program continually touches each page of the data structures in sequential order, forcing them to the head of the disk buffer pool, so they don't get aged out of the queue.

= Modifying the Reference Header =

''NB: This feature is not yet complete''

To facilitate SAM RG values being set automatically in a production environment, we keep a header in the binary version of the reference. The header can be viewed and edited using the header subcommands here.

To view the header:

karma header -r phiX.fa

To view and edit the header:

karma header -r phiX.fa -e

= Other test and check capabilities =

Due to the size and complexity of Karma input, output and index files, various checks and tests are useful, so we include some diagnostics capabilities:

Tests for external files:

karma check [options...] file.bam file.fastq file.sam file.fa file.umfa

Tests internal to Karma:
karma test [options...]
-d -> debug
-s [int] -> set random number seed [12345]

= Karma File structure =

Upon successfully building references, you will obtain a list of reference files like below:

{| style="width: 571px; height: 288px" border="1" cellspacing="1" cellpadding="1" width="571"
|-
|
|
Base Space

| Color Space
|-
|
Reference genome

|
NCBI37-bs.umfa

| NCBI37-cs.umfa
|-
|
Word Index

|
NCBI37-bs.15.5000.umwiwp

NCBI37-bs.15.5000.umwihi

|
NCBI37-cs.15.5000.umwiwp

NCBI37-cs.15.5000.umwihi

|-
|
Word Hash (Left)

|
NCBI37-bs.15.5000.umwhl

| NCBI37-cs.15.5000.umwhl
|-
|
Word Hash (Right)

|
NCBI37-bs.15.5000.umwhr

| NCBI37-cs.15.5000.umwhr
|}

 

 

= Karma TODO List =

= Karma CHANGELOG =

= Other useful links =

[http://lh3lh3.users.sourceforge.net/bioinfo.shtml Heng Li's thoughts about aligner]

[http://lh3lh3.users.sourceforge.net/udb.shtml Benchmark of Dictionary Structures]

Karma

2010-04-08T19:03:00Z

Pha: fill in performance details

[[Category:Software]]

'''K-tuple Alignment with Rapid Matching Algorithm'''

Karma uses an existing reference to align short reads, such as generated by Illumina sequencers.

The current version, 0.9.0, is optimized to rapidly map base space reads from Illumina sequencers.

Color space and LS454 sequence alignments are not working. These features will return in Karma 0.9.1.

= Download Karma =

To get a copy go to [http://www.sph.umich.edu/csg/pha/karma/download/ Karma Download]

= Build Karma =

== Dependencies ==

== Building ==

== Testing the build ==

To test karma, go to the build tree subdirectory named ''karma'', and type the command:

make test

The test script builds a reference for the small phiX genome, then runs single end as well as paired end alignments. It compares the results of that with known results. Differences are printed to the console, and currently look something like this:

<pre>
diff phiX.sam.good phiX.sam
3c3
< @RG DT:2010-04-08T17:29Z ID:boingboing SM:NA12345
---
> @RG DT:2010-04-08T18:13Z ID:boingboing SM:NA12345
</pre>

Any differences greater than that are an error and need to be fixed by the author.

= Normal Workflow =

Karma works using a set of index and hash files created from an existing reference. Once created, this set of reference index and hash files must always be specified in the command line when aligning reads.

In concept, the simplest workflow is to first create a reference index using ''karma create'', then align reads using ''karma map''. You only have to build the index and hash once.

Because the reference can be large, and because Karma will share the reference among many running instances of Karma, it is useful to put well known references in a common location readily accessible to you and your collaborators.

= Build reference index and hash =

Building a reference index and hash with Karma is straightforward, but because it is time consuming for longer genomes, you typically save the reference index between runs.

The simplest example for creating a reference and index using a wordsize of 11-mer words is:

karma create -i -w 11 phiX.fa

More generally, three primary parameters are necessary for building a Karma reference index:

# a boolean flag indicating base or color space
# the index table word occurrence cutoff value
# the word size

Although the input reference is always expected to be base space and in [http://en.wikipedia.org/wiki/FASTA_format FASTA] format, the binary version of the reference, and the corresponding index and hash files, can be in either color space (ABI SOLiD) or base space (Illumina or LS454). For a given reference [http://en.wikipedia.org/wiki/FASTA_format FASTA] file, you may have either a color or base space binary reference, as well as either color or base space index/hash files, any in varying word sizes or occurrence cutoffs.

Because the index and hash files are dependent on the occurrence cutoff parameter and the word size, the output files created by karma have those values in the file name. This allows you to create a variety of index/hash tables, depending on your expected use (ABI SOLiD, in particular, is sensitive to read length).

== Options for building reference index and hash ==

-r ''reference'' Reference file in [http://en.wikipedia.org/wiki/FASTA_format FASTA] format
-w ''word size'' Word size for index and hash (default 15, typically 10-16)
-O ''occurrence cutoff'' Upper count of number of word positions to store in word positions table (default 5000)
-c Creates a color space reference and index/hash
-i Create the index and hash as well as the binary reference

= Aligning Reads =

Aligning reads to the reference is easy:

karma map -r phiX.fa -w 11 phiX.fastq

or for paired reads:

karma map -r phiX.fa -w 11 phiX-mate1.fastq phiX-mate2.fastq

In both of the above examples, the -r option names the reference originally used to build the index/hash, and the -w 11 specifies that we are using the index/hash built for 11-mer words. Although you can use the default word size of 15 for phiX, the index is 4^15 * 4 = 4GBytes, so a shorter word size is prudent.

Since Karma uses the word size and occurrence cutoff to help construct the actual index and hash filenames, you must specify them the same way you did when you created the reference index and hash.

== Aligning Reads (Illumina) ==

Karma is set up so that the default options work well for mapping Illumina reads to the Human genome.

== Aligning Reads (ABI SOLiD) ==

Karma has been designed to align color space reads. However, in Karma 0.9.0, this functionality is not working.

== Aligning Reads (LS 454) ==

Karma has been designed to align LS 454 reads. However, in Karma 0.9.0, this functionality is not working.

= Karma Performance Tuning =

There are four components to the Karma index and hash. A pure index, based on an N-mer word index. This is used as a pointer into a word positions table, which is an ordered list of genome positions in which that N-mer word appears. There is a cap called the ''occurrence cutoff'', which once crossed, causes that index word to be marked as a high repeat pattern. Once marked as high repeat, the N-mer word is now combined with both the N-mer word preceding it, as well as the N-mer word succeeding it to create a 2 * N-mer word hash key. Two hash tables are populated, a left and a right hash.

== Index Word Size ==

Choosing an appropriate word size for larger genome is critical to performance. The easiest case is for Illumina base space reads with the human genome (3Gbases), where the default 15-mer word size is fine.

For smaller genomes, consider using a smaller word size. Genomes smaller than a few million bases should be perfectly fine with a word size of 11 or 12.

Since the primary index table into the word positions table is 2^(wordsize) * 4 bytes, it can grow large rapidly. All else being equal, a smaller word size leads to longer sets of word positions for each index value. Each increment of word size approximately quadruples storage requirements, and halves runtime. Similarly, each decrement of word size reduces the index table size by 75%, and doubles runtime. These approximations are old, but serve a useful rule of thumb.

For ABI SOLiD reads, the word size is critical, due to the shorter length of reads as compared to Illumina or LS 454.

The optimal minimum word size is chosen such that it is 1/4 the minimum expected average read length. It also must be chosen to be 1/2 the minimum expected read length, since at least 2 full words must exist in the read.

So for 48-mer reads, a reasonable value of word size is 12. Although the base space default of 15 is fine, too, Karma is able to take advantage of a higher number of index words per read, yielding substantial speedups even with the shorter read. Similarly, 52-mer reads would map better with a 13-mer word size, and 56-mer reads would map best with a 14-mer word size.

== Occurrence Cutoff ==

The occurrence cutoff value determines how quickly an N-mer pattern is declared to be ''high repeat'' and left out of the index in favor of a hash. The default value of 5000 seems adequate for Illumina reads with the human genome. If ultimate performance is necessary, some experimentation is called for with this value.

== Shared Memory ==

Karma uses memory mapped files to share the potentially large reference index and hash data structures.

Karma uses this to great effect on our 8 processors with hyperthreading enabled. 16 copies of karma can share one reference index and hash, yielding a very acceptable memory per CPU ratio of around 1GB/CPU.

A problem with large reference index and hash data structures is that they are more prone to being paged out. On a shared machine that is being used extensively even just simple disk I/O, memory pages are being reclaimed such that Karma will become swapped out.

While Karma can recover on its own, it is best to either run in a production manner on dedicated machines, or to run a program such as the utility ''mapfile'' found in the utilities sub-folder. This program continually touches each page of the data structures in sequential order, forcing them to the head of the disk buffer pool, so they don't get aged out of the queue.

= Modifying the Reference Header =

''NB: This feature is not yet complete''

To facilitate SAM RG values being set automatically in a production environment, we keep a header in the binary version of the reference. The header can be viewed and edited using the header subcommands here.

To view the header:

karma header -r phiX.fa

To view and edit the header:

karma header -r phiX.fa -e

= Other test and check capabilities =

Due to the size and complexity of Karma input, output and index files, various checks and tests are useful, so we include some diagnostics capabilities:

Tests for external files:

karma check [options...] file.bam file.fastq file.sam file.fa file.umfa

Tests internal to Karma:
karma test [options...]
-d -> debug
-s [int] -> set random number seed [12345]

= Karma File structure =

Upon successfully building references, you will obtain a list of reference files like below:

{| style="width: 571px; height: 288px" border="1" cellspacing="1" cellpadding="1" width="571"
|-
|
|
Base Space

| Color Space
|-
|
Reference genome

|
NCBI37-bs.umfa

| NCBI37-cs.umfa
|-
|
Word Index

|
NCBI37-bs.15.5000.umwiwp

NCBI37-bs.15.5000.umwihi

|
NCBI37-cs.15.5000.umwiwp

NCBI37-cs.15.5000.umwihi

|-
|
Word Hash (Left)

|
NCBI37-bs.15.5000.umwhl

| NCBI37-cs.15.5000.umwhl
|-
|
Word Hash (Right)

|
NCBI37-bs.15.5000.umwhr

| NCBI37-cs.15.5000.umwhr
|}

 

 

= Karma TODO List =

= Karma CHANGELOG =

= Other useful links =

[http://lh3lh3.users.sourceforge.net/bioinfo.shtml Heng Li's thoughts about aligner]

[http://lh3lh3.users.sourceforge.net/udb.shtml Benchmark of Dictionary Structures]

Karma

2010-04-08T18:36:26Z

Pha: add mapping examples

Karma

2010-04-08T18:25:15Z

Pha: populate alignment section

Karma

2010-04-08T18:07:48Z

Pha: /* Testing the build */ more detail

Karma

2010-04-08T18:07:05Z

Pha: /* Testing the build */ flesh it out