NAME
Maq - Mapping and Assembly with Qualities
SYNOPSIS
maq command [options] arguments
maq.pl command [options] arguments
DESCRIPTION
Maq is a software that builds mapping assemblies from short reads
generated by the next-generation sequencing machines. It is
particularly designed for Illumina-Solexa 1G Genetic Analyzer, and has
a preliminary functionality to handle AB SOLiD data.
With Maq you can:
· Fast align Illumina/SOLiD reads to the reference genome. With the
default options, one million pairs of reads can be mapped to the
human genome in about 10 CPU hours with less than 1G memory.
· Accurately measure the error probability of the alignment of each
individual read.
· Call the consensus genotypes, including homozygous and heterozygous
polymorphisms, with a Phred probabilistic quality assigned to each
base.
· Find short indels with paired end reads.
· Accurately find large scale genomic deletions and translocations with
paired end reads.
· Discover potential CNVs by checking read depth.
· Evaluate the accuracy of raw base qualities from sequencers and help
to check the systematic errors.
However, Maq can NOT:
· Do de novo assembly. (Maq can only call the consensus by mapping
reads to a known reference.)
· Map shorts reads against themselves. (Maq can only find complete
overlap between reads.)
· Align capillary reads or 454 reads to the reference. (Maq cannot
align reads longer than 63bp.)
MAQ COMMANDS
Key Commands
fasta2bfa maq fasta2bfa in.ref.fasta out.ref.bfa
Convert sequences in FASTA format to Maq’s BFA (binary
FASTA) format.
fastq2bfq maq fastq2bfq [-n nreads] in.read.fastq
out.read.bfq│out.prefix
Convert reads in FASTQ format to Maq’s BFQ (binary FASTQ)
format.
OPTIONS:
-n INT number of reads per file [not specified]
map maq map [-n nmis] [-a maxins] [-c] [-1 len1] [-2 len2] [-d
adap3] [-m mutrate] [-u unmapped] [-e maxerr] [-M c│g] [-N]
[-H allhits] [-C maxhits] out.aln.map in.ref.bfa
in.read1.bfq [in.read2.bfq] 2> out.map.log
Map reads to the reference sequences.
OPTIONS:
-n INT Number of maximum mismatches that can always be
found [2]
-a INT Maximum outer distance for a correct read pair
[250]
-A INT Maximum outer distance of two RF paied read (0 for
disable) [0]
-c Map reads in the colour space (for SOLiD only)
-1 INT Read length for the first read, 0 for auto [0]
-2 INT Read length for the second read, 0 for auto [0]
-m FLOAT Mutation rate between the reference sequences and
the reads [0.001]
-d FILE Specify a file containing a single line of the
3’-adapter sequence [null]
-u FILE Dump unmapped reads and reads containing more than
nmis mismatches to a separate file [null]
-e INT Threshold on the sum of mismatching base qualities
[70]
-H FILE Dump multiple/all 01-mismatch hits to FILE [null]
-C INT Maximum number of hits to output. Unlimited if
larger than 512. [250]
-M c│g methylation alignment mode. All C (or G) on the
forward strand will be changed to T (or A). This
option is for testing only.
-N store the mismatch position in the output file
out.aln.map. When this option is in use, the
maximum allowed read length is 55bp.
NOTE:
* Paired end reads should be prepared in two files, one for
each end, with reads are sorted in the same order. This
means the k-th read in the first file is mated with the
k-th read in the second file. The corresponding read names
must be identical up to the tailing ‘/1’ or ‘/2’. For
example, such a pair of read names are allowed:
‘EAS1_1_5_100_200/1’ and ‘EAS1_1_5_100_200/2’. The tailing
‘/[12]’ is usually generated by the GAPipeline to
distinguish the two ends in a pair.
* The output is a compressed binary file. It is affected by
the endianness.
* The best way to run this command is to provide about 1 to
3 million reads as input. More reads consume more memory.
* Option -n controls the sensitivity of the alignment. By
default, a hit with up to 2 mismatches can be always
found. Higher -n finds more hits and also improves the
accuracy of mapping qualities. However, this is done at
the cost of speed.
* Alignments with many high-quality mismatches should be
discarded as false alignments or possible contaminations.
This behaviour is controlled by option -e. The -e
threshold is only calculated approximately because base
qualities are divided by 10 at a certain stage of the
alignment. The -Q option in the assemble command precisely
set the threshold.
* A pair of reads are said to be correctly paired if and
only if the orientation is FR and the outer distance of
the pair is no larger than maxins. There is no limit on
the minimum insert size. This setting is determined by the
paired end alignment algorithm used in Maq. Requiring a
minimum insert size will lead to some wrong alignments
with highly overestimated mapping qualities.
* Currently, read pairs from Illumina/Solexa long-insert
library have RF read orientation. The maximum insert size
is set by option -A. However, long-insert library is also
mixed with a small fraction of short-insert read pairs. -a
should also be set correctly.
* Sometimes 5’-end or even the entire 3’-adapter sequence
may be sequenced. Providing -d renders Maq to eliminate
the adapter contaminations.
* Given 2 million reads as input, maq usually takes 800MB
memory.
mapmerge maq mapmerge out.aln.map in.aln1.map in.aln2.map [...]
Merge a batch of read alignments together.
NOTE:
* In theory, this command can merge unlimited number of
alignments. However, as mapmerge will be reading all the
inputs at the same time, it may hit the limit of the
maximum number of opening files set by the OS. At present,
this has to be manually solved by endusers.
* Command mapmerge can be used to merge alignment files with
different read lengths. All the subsequent analyses do not
assume fixed length any more.
rmdup maq rmdup out.rmdup.map in.ori.map
Remove pairs with identical outer coordinates. In principle,
pairs with identical outer coordinates should happen rarely.
However, due to the amplification in sample preparation,
this occurs much more frequently than by chance. Practical
analyses show that removing duplicates helps to improve the
overall accuracy of SNP calling.
assemble maq assemble [-sp] [-m maxmis] [-Q maxerr] [-r hetrate] [-t
coef] [-q minQ] [-N nHap] out.cns in.ref.bfa in.aln.map 2>
out.cns.log
Call the consensus sequences from read mapping.
OPTIONS:
-t FLOAT Error dependency coefficient [0.93]
-r FLOAT Fraction of heterozygotes among all sites [0.001]
-s Take single end mapping quality as the final
mapping quality; otherwise paired end mapping
quality will be used
-p Discard paired end reads that are not mapped in
correct pairs
-m INT Maximum number of mismatches allowed for a read to
be used in consensus calling [7]
-Q INT Maximum allowed sum of quality values of mismatched
bases [60]
-q INT Minimum mapping quality allowed for a read to be
used in consensus calling [0]
-N INT Number of haplotypes in the pool (>=2) [2]
NOTE:
* Option -Q sets a limit on the maximum sum of mismatching
base qualities. Reads containing many high-quality
mismatches should be discarded.
* Option -N sets the number of haplotypes in a pool. It is
designed for resequencing of samples by pooling multiple
strains/individuals together. For diploid genome
resequencing, this option equals 2.
glfgen maq glfgen [-sp] [-m maxmis] [-Q maxerr] [-r hetrate] [-t
coef] [-q minQ] [-N nHap] out.cns in.ref.bfa in.aln.map 2>
out.cns.log
Calculate log-likelihood for all genotypes and store the
results in GLF format (Genotyping Likelihood Format). Please
check MAQ website for detailed descriptions of the file
format and the related utilities.
indelpe maq indelpe in.ref.bfa in.aln.map > out.indelpe
Call consistent indels from paired end reads. The output is
TAB delimited with each line consisting of chromosome, start
position, type of the indel, number of reads across the
indel, size of the indel and inserted/deleted nucleotides
(separated by colon), number of indels on the reverse
strand, number of indels on the forward strand, 5’ sequence
ahead of the indel, 3’ sequence following the indel, number
of reads aligned without indels and three additional columns
for filters.
At the 3rd column, type of the indel, a star indicates the
indel is confirmed by reads from both strands, a plus means
the indel is hit by at least two reads but from the same
strand, a minus shows the indel is only found on one read,
and a dot means the indel is too close to another indel and
is filtered out.
Users are recommended to run through ‘maq.pl indelpe’ to
correct the number of reads mapped without indels. For more
details, see the ‘maq.pl indelpe’ section.
indelsoa maq indelsoa in.ref.bfa in.aln.map > out.indelsoa
Call potential homozygous indels and break points by
detecting the abnormal alignment pattern around indels and
break points. The output is also TAB delimited with each
line consisting of chromosome, approximate coordinate,
length of the abnormal region, number of reads mapped across
the position, number of reads on the left-hand side of the
position and number of reads on the right-hand side. The
last column can be ignored.
The output contains many false positives. A recommended
filter could be:
awk ’$5+$6-$4 >= 3 && $4 <= 1’ in.indelsoa
Note that this command does not aim to be an accurate indel
detector, but mainly helps to avoid some false positives in
substitution calling. In addition, it only works well given
deep depth (~40X for example); otherwise the false negative
rate would be very high.
Format Converting
sol2sanger maq sol2sanger in.sol.fastq out.sanger.fastq
Convert Solexa FASTQ to standard/Sanger FASTQ format.
bfq2fastq maq bfq2fastq in.read.bfq out.read.fastq
Convert Maq’s BFQ format to standard FASTQ format.
mapass2maq maq mapass2maq in.mapass2.map out.maq.map
Convert obsolete mapass2’s map format to Maq’s map format.
The old format does not contain read names.
Information Extracting
mapview maq mapview [-bN] in.aln.map > out.aln.txt
Display the read alignment in plain text. For reads aligned
before the Smith-Waterman alignment, each line consists of
read name, chromosome, position, strand, insert size from
the outer coorniates of a pair, paired flag, mapping
quality, single-end mapping quality, alternative mapping
quality, number of mismatches of the best hit, sum of
qualities of mismatched bases of the best hit, number of
0-mismatch hits of the first 24bp, number of 1-mismatch hits
of the first 24bp on the reference, length of the read, read
sequence and its quality. Alternative mapping quality
always equals to mapping quality if the reads are not
paired. If reads are paired, it equals to the smaller
mapping quality of the two ends. This alternative mapping
quality is actually the mapping quality of an abnormal pair.
The fifth column, paired flag, is a bitwise flag. Its lower
4 bits give the orientation: 1 stands for FF, 2 for FR, 4
for RF, and 8 for RR, where FR means that the read with
smaller coordinate is on the forward strand, and its mate is
on the reverse strand. Only FR is allowed for a correct
pair. The higher bits of this flag give further information.
If the pair meets the paired end requirement, 16 will be
set. If the two reads are mapped to different chromosomes,
32 will be set. If one of the two reads cannot be mapped at
all, 64 will be set. The flag for a correct pair always
equals to 18.
For reads aligned by the Smith-Waterman alignment
afterwards, the flag is always 130. A line consists of read
name, chromosome, position, strand, insert size, flag
(always 130), position of the indel on the read (0 if no
indel), length of the indels (positive for insertions and
negative for deletions), mapping quality of its mate, number
of mismatches of the best hit, sum of qualities of
mismatched bases of the best hit, two zeros, length of the
read, read sequence and its quality. The mate of a
130-flagged read always gets a flag 18.
Flag 192 indicates that the read is not mapped but its mate
is mapped. For such a read pair, one read has flag 64 and
the other has 192.
OPTIONS:
-b do not display the read sequence and the quality
-N display the positions where mismatches occur. This
flag only works with a .map file generated by ‘maq
map -N’.
mapcheck maq mapcheck [-s] [-m maxmis] [-q minQ] in.ref.bfa
in.aln.map > out.mapcheck
Read quality check. The mapcheck first reports the
composition and the depth of the reference. After that there
is a form. The first column indicates the position on a
read. Following four columns which show the nucleotide
composition, substitution rates between the reference and
reads will be given. These rates and the numbers in the
following columns are scaled to 999 and rounded to nearest
integer. The next group of columns show the distribution of
base qualities along the reads at a quality interval of 10.
A decay in quality can usually be observed, which means
bases at the end of read are less accurate. The last group
of columns present the fraction of substitutions for read
bases at a quality interval. This measures the accuracy of
base quality estimation. Idealy, we expect to see 1 in the
3? column, 10 in the 2? column and 100 in the 1? column.
OPTIONS:
-s Take single end mapping quality as the final
mapping quality
-m INT Maximum number of mismatahces allowed for a read to
be counted [4]
-q INT Minimum mapping quality allowed for a read to be
counted [30]
pileup maq pileup [-spvP] [-m maxmis] [-Q maxerr] [-q minQ] [-l
sitefile] in.ref.bfa in.aln.map > out.pileup
Display the alignment in a ‘pileup’ text format. Each line
consists of chromosome, position, reference base, depth and
the bases on reads that cover this position. If -v is added
on the command line, base qualities and mapping qualities
will be presented in the sixth and seventh columns in order.
The fifth column always starts with ‘@’. In this column,
read bases identical to the reference are showed in comma
‘,’ or dot ‘.’, and read bases different from the reference
in letters. A comma or a upper case indicates that the base
comes from a read aligned on the forward strand, while a dot
or a lower case on the reverse strand.
This command is for users who want to develop their own SNP
callers.
OPTIONS:
-s Take single end mapping quality as the final
mapping quality
-p Discard paired end reads that are not mapped as
correct pairs
-v Output verbose information including base qualities
and mapping qualities
-m INT Maximum number of mismatches allowed for a read to
be used [7]
-Q INT Maximum allowed number of quality values of
mismatches [60]
-q INT Minimum mapping quality allowed for a read to be
used [0]
-l FILE File containing the sites at which pileup will be
printed out. In this file the first column gives
the names of the reference and the second the
coordinates. Additional columns will be ignored.
[null]
-P also output the base position on the read
cns2fq maq cns2fq [-Q minMapQ] [-n minNeiQ] [-d minDepth] [-D
maxDepth] in.cns > out.cns.fastq
Extract the consensus sequences in FASTQ format. In the
sequence lines, bases in lower case are essentially repeats
or do not have sufficient coverage; bases in upper case
indicate regions where SNPs can be reliably called. In the
quality lines, ASCII of a character minus 33 gives the PHRED
quality.
OPTIONS:
-Q INT Minimum mapping quality [40]
-d INT Minimum read depth [3]
-n INT Minimum neighbouring quality [20]
-D INT Maximum read dpeth. >=255 for unlimited. [255]
cns2snp maq cns2snp in.cns > out.snp
Extract SNP sites. Each line consists of chromosome,
position, reference base, consensus base, Phred-like
consensus quality, read depth, the average number of hits of
reads covering this position, the highest mapping quality of
the reads covering the position, the minimum consensus
quality in the 3bp flanking regions at each side of the site
(6bp in total), the second best call, log likelihood ratio
of the second best and the third best call, and the third
best call.
The 5th column is the key criterion when you judge the
reliability of a SNP. However, as this quality is only
calculated assuming site independency, you should also
consider other columns to get more accurate SNP calls.
Script command ‘maq.pl SNPfilter’ is designed for this (see
below).
The 7th column implies whether the site falls in a
repetitive region. If no read covering the site can be
mapped with high mapping quality, the flanking region is
possibly repetitive or in the lack of good reads. A SNP at
such site is usually not reliable.
The 8th column roughly gives the copy number of the flanking
region in the reference genome. In most cases, this number
approaches 1.00, which means the region is about unique.
Sometimes you may see non-zero read depth but 0.00 at the
7th column. This indicates that all the reads covering the
position have at least two mismatches. Maq only counts the
number of 0- and 1-mismatch hits to the reference. This is
due to a complex technical issue.
The 9th column gives the neighbouring quality. Filtering on
this column is also required to get reliable SNPs. This idea
is inspired by NQS, although NQS is initially designed for a
single read instead of a consensus.
cns2view maq cns2view in.cns > out.view
Show detailed information at all sites. The output format is
identical to cns2snp report.
cns2ref maq cns2ref in.cns > out.ref.fasta
Extract the reference sequence.
cns2win maq cns2win [-w winsize] [-c chr] [-b begin] [-e end] [-q
minQ] in.cns > out.win
Extract information averaged in a tilling window. The output
is TAB delimited, which consists of reference name,
coordinate divided by 1,000,000, SNP rate, het rate, raw
read depth, read depth in approximately unique regions, the
average number of hits of reads in the window and percent
GC.
OPTIONS:
-w INT Size of a window [1000]
-c STR Destinated reference sequence; otherwise all
references will be used [null]
-b INT Start position, 0 for no constraint [0]
-e INT End position, 0 for no constraint [0]
-q INT Minimum consensus quality of the sites to be used
[0]
Simulation Related
fakemut maq fakemut [-r mutrate] [-R indelfrac] in.ref.fasta >
out.fakeref.fasta 2> out.fake.snp
Randomly introduce substitutions and indels to the
reference. Substitutions and sinlge base-pair indels can be
added.
OPTIONS:
-r FLOAT Mutation rate [0.001]
-R FLOAT Fraction of mutations to be indels [0.1]
simutrain maq simutrain out.simupars.dat in.read.fastq
Estimate/train parameters for read simulation.
simulate maq simulate [-d insize] [-s stdev] [-N nReads] [-1
readLen1] [-2 readLen2] [-r mutRate] [-R indelFrac] [-h]
out.read1.fastq out.read2.fastq in.ref.fasta in.simupars.dat
Simulate paired end reads. File in.simupars.dat determines
the read lengths and quality distribution. It is generated
from simutrain, or can be downloaded from Maq website. In
the output read files, a read name consists of the reference
sequence name and the outer coordinates of the pair of
simulated reads. By default, simulate assumes reads come
from a diploid sequence which is generated by adding two
different sets of mutations, including one base-pair indels,
to in.ref.fasta.
OPTIONS:
-d INT mean of the outer distance of insert sizes [170]
-s INT standard deviation of insert sizes [20]
-N INT number of pairs of reads to be generated [1000000]
-1 INT length of the first read [set by in.simupars.dat]
-2 INT length of the second read [set by in.simupars.dat]
-r FLOAT mutation rate [0.001]
-R FLOAT fraction of 1bp indels [0.1]
-h add all mutations to in.ref.fasta and generate
reads from the single mutated sequence (haploid
mode)
NOTE:
* Reads generated from this command are independent, which
deviates from the truth. Whereas alignment evaluation is
less affected by this, evaluation on SNP calling should be
performed with caution. Error dependency may be one of the
major causes of wrong SNP calls.
simustat maq simustat in.simu-aln.map > out.simustat
Evaluate mapping qualities from simulated reads.
SOLiD Related
fasta2csfa maq fasta2csfa in.nucl-ref.fasta > out.colour-ref.fasta
Convert nucleotide FASTA to colour-coded FASTA. Flag -c
should be then applied to map command. In the output, letter
‘A’ stands for color 0, ‘C’ for 1, ‘G’ for 2 and ‘T’ for 3.
Each sequence in the output is 1bp shorter than the input.
csmap2nt maq csmap2nt out.nt.map in.ref.nt.bfa in.cs.map
Convert color alignment to nucleotide alignment. The input
in.ref.nt.bfa is the nucleotide binary FASTA reference file.
It must correspond to the original file from which the color
reference is converted. Nucleotide consensus can be called
from the resultant alignment.
Miscellaneous/Advanced Commands
submap maq submap [-q minMapQ] [-Q maxSumErr] [-m maxMM] [-p]
out.map in.map
Filter bad alignments in in.map. Command-line options are
described in the ‘assemble’ command.
eland2maq maq eland2maq [-q defqual] out.map in.list in.eland
Convert eland alignment to maq’s .map format. File in.list
consists of the sequence names that appear at the seventh
column of the eland alignment file in.eland and the name you
expect to see in maq alignment. The following is an example:
cX.fa chrX
c1.fa chr1
c2.fa chr2
If you are aligning reads in several batches using eland, it
is important to use the same in.list for the conversion. In
addition, maq will load all the alignments and sort them in
the memory. If you have concatenate several eland outputs
into one huge file, you should separate it into smaller
files to prevent maq from eating all your machine memory.
This command actually aims to show Eland alignment in
Maqview. As no quality information is available, the
resultant maq alignment file should not be used to call
consensus genotypes.
export2maq maq export2maq [-1 read1len] [-2 read2len] [-a maxdist] [-n]
out.map in.list in.export
Convert Illumina’s Export format to Maq’s .map format.
Export format is a new alignment format since
SolexaPipeline-0.3.0 which also calculates mapping qualities
like maq. The resultant file can be used to call consensus
genotypes as most of necessary information is available for
maq to do this accurately.
OPTIONS:
-1 INT Length of the first read [0]
-2 INT Length of the second read [0]
-a INT Maximum outer distance for a correct read pair
[250]
-n Retain filtered reads
MAQ-PERL COMMANDS
demo maq.pl demo [-h] [-s] [-N nPairs] [-d outDir] in.fasta
in.simudat
Demonstrate the use of maq and its companion scripts. This
command will simulate reads from a FASTA file in.fasta. The
sequence length and qualities are determined by in.simudat
which is generated from maq simutrain or can be downloaded
from Maq website. The simulated reads will then be mapped
with maq.pl easyrun. The alignment accuracy is evaluated by
maq simustat, the consensus accuracy by maq simucns, and the
SNP accuracy by maq_eval.pl.
By default, paired end reads will be simulated and a diploid
sequence will be generated from the input by adding
mutations to either haploid type. The insert size and
mutation rate are controlled by maq simulate.
OPTIONS:
-h simulate a haploid sequence instead of a diploid
sequence
-s use single-end mode to align reads instead of
paired-end mode
-N INT number of pairs of reads to be simulated [1000000]
-d DIR output directory [maqdemo]
NOTE:
* The output files from maq_eval.pl have not been
documented, but you may make a good guess at some of these
files.
* This command just demonstrates the use of the maq suite.
The accuracy on real data is almost always lower than what
you see from pure simulation.
easyrun maq.pl easyrun [-1 read1Len] [-d out.dir] [-n nReads] [-A
3adapter] [-e minDep] [-q minCnsQ] [-p] [-2 read2Len] [-a
maxIns] [-S] [-N] in.ref.fasta in1.fastq [in2.fastq]
Analyses pipeline for small genomes. Easyrun command will
run most of analyses implemented in maq. By default, easyrun
assumes all the input read sequences files are single-end
and independent; when -p is specified, two read sequence
files are required, one for each end.
Several files will be generated in out.dir, among which the
following files are the key output:
cns.final.snp final SNP calls with low quality ones
filtered out
cns.fq consensus sequences and qualities in the
FASTQ format
OPTIONS:
-d DIR output directory [easyrun]
-n INT number of reads/pairs in one batch of alignment
[2000000]
-S apply split-read analysis of short indels (maybe
very slow)
-N INT number of haplotypes/strains in the pool (>=2) [2]
-A FILE file for 3’-adapter. The file should contain a
single line of sequence [null]
-1 INT length of the first read, 0 for auto [0]
-e INT minimum read depth required to call a SNP (for
SNPfilter) [3]
-q INT minimum consensus quality for SNPs in cns.final.snp
[30]
-p switch to paired end alignment mode
-2 INT length of the second read when -p is applied [0]
-a INT maximum insert size when -p is applied [250]
NOTES:
* For SNP calling on pooled samples, users should set
correct ‘-N’ as well as ‘-E 0’.
* The input file can be maq’s binary format. maq.pl will
automatically detect the file format.
SNPfilter maq.pl SNPfilter [-d minDep] [-D maxDep] [-Q maxMapQ] [-q
minCnsQ] [-w indelWinSize] [-n minNeiQ] [-F in.indelpe] [-f
in.indelsoa] [-s minScore] [-m maxAcross] [-a] [-N
maxWinSNP] [-W densWinSize] in.cns2snp.snp >
out.filtered.snp
Rule out SNPs that are covered by few reads (specified by
-d), by too many reads (specified by -D), near (specified by
-w) to a potential indel, falling in a possible repetitve
region (characterized by -Q), or having low-quality
neighbouring bases (specified by -n). If maxWinSNP or more
SNPs appear in any densWinSize window, they will also be
filtered out together.
OPTIONS:
-d INT Minimum read depth required to call a SNP [3]
-D INT Maximum read depth required to call a SNP (<255,
otherwise ignored) [256]
-Q INT Required maximum mapping quality of reads covering
the SNP [40]
-q INT Minimum consensus quality [20]
-n INT Minimum adjacent consensus quality [20]
-w INT Size of the window around the potential indels.
SNPs that are close to indels will be suppressed
[3]
-F FILE The indelpe output [null]
-f FILE The indelsoa output [null]
-s INT Minimum score for a soa-indel to be considered [3]
-m INT Maximum number of reads that can be mapped across
a soa-indel [1]
-a Alternative filter for single end alignment
indelpe maq.pl indelpe in.indelpe > out.indelpe
Correct the number of reads mapped without indels for
homopolymer tracts. This command modify the 4th, 10th and
the last three columns of in.indelpe and output the result
in out.indelpe. After the correction, the following awk
command gives putative homozygous indels:
awk ’($3=="*"││$3=="+") && $6+$7>=3 && ($6+$7)/$4>=0.75’
and the following gives heterozygotes:
awk ’($3=="*"││$3=="+") && $6+$7>=3 && ($6+$7)/$4<0.75’
Please note that this indelpe command just implements
several heuristic rules. It does not correct for impure
homopolymer runs or di-nucleotide/triplet repeats.
Consequently, the two awk commands only give approximate
hom/het indels.
EXAMPLES
· Easyrun script:
maq.pl easyrun -d easyrun ref.fasta part1.fastq part2.fastq
· Key commands behind easyrun:
maq fasta2bfa ref.fasta ref.bfa;
maq fastq2bfq part1.fastq part1.bfq;
maq fastq2bfq part2.fastq part2.bfq;
maq map part1.map ref.bfa part1.bfq;
maq map part2.map ref.bfa part2.bfq;
maq mapmerge aln.map part1.map part2.map;
maq assemble cns.cns ref.bfa aln.map;
LICENSE
GNU General Public License, version 3 (GPLv3)
AVAILABILITY
<http://maq.sourceforge.net>
AUTHOR
Heng Li <lh3@sanger.ac.uk>