Divide-and-conquer in a #Makefile : recursivity and #parallelism.

This post is my notebook about implementing a divide-and-conquer strategy in GNU make.
Say you have a list of 'N' VCFs files. You want to create a list of:

• common SNPs in vcf1 and vcf2
• common SNPs in vcf3 and the previous list
• common SNPs in vcf4 and the previous list
• (...)
• common SNPs in vcfN and the previous list

Yes, I know I can do this using:grep -v '^#' f.vcf|cut -f 1,2,4,5 | sort | uniq

Using a linear Makefile it could look like:

list2: vcf1 vcf2
    grep -v '^#' $^ |cut -f 1,2,4,5 | sort | uniq > $@
list3: vcf3 list2
    grep -v '^#' $^ |cut -f 1,2,4,5 | sort | uniq > $@
list4: vcf4 list3
    grep -v '^#' $^ |cut -f 1,2,4,5 | sort | uniq > $@
list5: vcf5 list4
    grep -v '^#' $^ |cut -f 1,2,4,5 | sort | uniq > $@
(...)

We can speed-up the workflow using the parallel option of make -j (number-of-parallel-jobs) and using a divide-and-conquer strategy. Here, the targets 'list1_2' and 'list3_4' can be processed independently in parallel.

list1_2: vcf1 vcf2
    grep -v '^#' $^ |cut -f 1,2,4,5 | sort | uniq > $@

list3_4: vcf3 vcf4
    grep -v '^#' $^ |cut -f 1,2,4,5 | sort | uniq > $@

list1_4: list1_2 list3_4
    grep -v '^#' $^ |cut -f 1,2,4,5 | sort | uniq > $@

By using the internal 'Make' functions $(eval), $(call), $(shell) we can define a recursive method "recursive". This method takes two arguments which are the 0-based indexes of a VCF in the list of VCFs. Here is the Makefile:
Running the makefile:

$ make 
gunzip -c Sample18.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target0_1
gunzip -c Sample13.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target1_2
LC_ALL=C comm -12 target0_1 target1_2 > target0_2
gunzip -c Sample1.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target2_3
gunzip -c Sample19.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target3_4
gunzip -c Sample12.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target4_5
LC_ALL=C comm -12 target3_4 target4_5 > target3_5
LC_ALL=C comm -12 target2_3 target3_5 > target2_5
LC_ALL=C comm -12 target0_2 target2_5 > target0_5
gunzip -c Sample17.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target5_6
gunzip -c Sample16.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target6_7
LC_ALL=C comm -12 target5_6 target6_7 > target5_7
gunzip -c Sample9.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target7_8
gunzip -c Sample15.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target8_9
gunzip -c Sample5.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target9_10
LC_ALL=C comm -12 target8_9 target9_10 > target8_10
LC_ALL=C comm -12 target7_8 target8_10 > target7_10
LC_ALL=C comm -12 target5_7 target7_10 > target5_10
LC_ALL=C comm -12 target0_5 target5_10 > target0_10
gunzip -c Sample14.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target10_11
gunzip -c Sample3.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target11_12
LC_ALL=C comm -12 target10_11 target11_12 > target10_12
gunzip -c Sample11.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target12_13
gunzip -c Sample2.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target13_14
gunzip -c Sample6.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target14_15
LC_ALL=C comm -12 target13_14 target14_15 > target13_15
LC_ALL=C comm -12 target12_13 target13_15 > target12_15
LC_ALL=C comm -12 target10_12 target12_15 > target10_15
gunzip -c Sample20.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target15_16
gunzip -c Sample10.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target16_17
LC_ALL=C comm -12 target15_16 target16_17 > target15_17
gunzip -c Sample4.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target17_18
gunzip -c Sample8.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target18_19
gunzip -c Sample7.vcf.gz | grep -v '^#' | cut -f 1,2,4,5 | LC_ALL=C sort | uniq > target19_20
LC_ALL=C comm -12 target18_19 target19_20 > target18_20
LC_ALL=C comm -12 target17_18 target18_20 > target17_20
LC_ALL=C comm -12 target15_17 target17_20 > target15_20
LC_ALL=C comm -12 target10_15 target15_20 > target10_20
LC_ALL=C comm -12 target0_10 target10_20 > target0_20

and here is the generated workflow (drawn with make2graph ).
:
That's it
Pierre.

Notes about bwa 0.6.2 and the multiple hits.

Here are my notes about the way bwa 0.6.2 handles the multiple hits.


I've created a sub-sequence (~6Mb ) of the human chr22 (see the Makefile below). The sequence doesn't contain any base 'N' or any lowercase (=repeatMasked) letter.
1E6 reads (forward and reverse) have been generated from this ~chr22 with samtools wgsim. No mutation and no sequencing error have been allowed.

This ~chr22 sequence have been duplicated and renamed as 'DUP'. Both sequences have been merged in one fasta file named 'twoChrom.fa'

The reads have been aligned with BWA (v. 0.6.2-r126) on this genome.

There's only one SAM alignment per read

$ samtools view align.sorted.bam  | wc -l
2000000

1% of the reads were not properly paired/mapped on the genome

$ samtools view  -f 2 align.sorted.bam | wc -l
1999926

The properly mapped reads contain some informations about the alternative hits in 'XA'

$samtools view  -f 2 align.sorted.bam | grep 22_10_402_0 | verticalize -n

>>> 1
$1   22_10_402_0:0:0_0:0:0_c1608
$2   99
$3   22
$4   10
$5   0
$6   70M
$7   =
$8   333
$9   393
$10  GGGACGGTCATGCAATCTGGACAACATTCACCTTTAAAAGTTTATTGATCTTTTGTGACATGCACGTGGG
$11  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
$12  XT:A:R
$13  NM:i:0
$14  SM:i:0
$15  AM:i:0
$16  X0:i:2
$17  X1:i:0
$18  XM:i:0
$19  XO:i:0
$20  XG:i:0
$21  MD:Z:70
$22  XA:Z:DUP,+10,70M,0;
<<< 1

>>> 2
$1   22_10_402_0:0:0_0:0:0_c1608
$2   147
$3   22
$4   333
$5   0
$6   70M
$7   =
$8   10
$9   -393
$10  ATACCTGCCAGATGAGTCACTGGCAAAAGGTGCTGCTCCCTGGTGAGGGAGAAACACCAGGGGCTGGGAG
$11  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
$12  XT:A:R
$13  NM:i:0
$14  SM:i:0
$15  AM:i:0
$16  X0:i:2
$17  X1:i:0
$18  XM:i:0
$19  XO:i:0
$20  XG:i:0
$21  MD:Z:70
$22  XA:Z:DUP,-333,70M,0;
<<< 2

50% mapped on chr22, 50% mapped on DUP

$ samtools view  -f 2 align.sorted.bam | cut -d '      ' -f 3 | sort | uniq -c
 999860 22
1000066 DUP

Some pairs have been unmapped because the mate have been mapped on the other chromosome !

$ samtools-0.1.18/samtools view  -F 2 align.sorted.bam | grep 22_1210321_1210924 | verticalize -n

>>> 1
$1   22_1210321_1210924_0:0:0_0:0:0_e28f5
$2   81
$3   22
$4   1210855
$5   0
$6   70M
$7   DUP
$8   1210321
$9   0
$10  AGAGACTGAGTCCGGCTAGAGAACAGGGTGGAGCCCCTTTGGACCTTAGAGCTGGGCCTTTGGGCCTTGG
$11  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
$12  XT:A:R
$13  NM:i:0
$14  SM:i:0
$15  AM:i:0
$16  X0:i:6
$17  X1:i:0
$18  XM:i:0
$19  XO:i:0
$20  XG:i:0
$21  MD:Z:70
$22  XA:Z:22,-1996342,70M,0;DUP,-1996342,70M,0;22,+2551778,70M,0;DUP,+2551778,70M,0;DUP,-1210855,70M,0;
<<< 1

>>> 2
$1   22_1210321_1210924_0:0:0_0:0:0_e28f5
$2   161
$3   DUP
$4   1210321
$5   0
$6   70M
$7   22
$8   1210855
$9   0
$10  CCCGGGCCGGATGGCTCGCCTGCGCGGCCAGCTCCGGGCCGAAGCGGCTTCGCGGTCCGAGGTGCCGCGG
$11  IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII
$12  XT:A:R
$13  NM:i:0
$14  SM:i:0
$15  AM:i:0
$16  X0:i:2
$17  X1:i:0
$18  XM:i:0
$19  XO:i:0
$20  XG:i:0
$21  MD:Z:70
$22  XA:Z:22,+1210321,70M,0;
<<< 2

The Makefile

SAMDIR=/usr/local/package/samtools-0.1.18
SAMTOOLS=$(SAMDIR)/samtools
BCFTOOLS=$(SAMDIR)/bcftools/bcftools
BWA=/usr/local/package/bwa-0.6.2/bwa
REF1=chr22.fa
REF2=twoChrom.fa

.INTERMEDIATE : align.sam random_1.sai random_2.sai align.bam variations.bcf

%.bam : %.sam
        $(SAMTOOLS) view -o $@ -b -S -T $(REF2) $<
%.bam.bai : %.bam
        $(SAMTOOLS) index $<



align.sorted.bam : align.bam
        $(SAMTOOLS) sort $< align.sorted


align.sam : random_1.sai random_2.sai  
        $(BWA) sampe -a 600 $(REF2) $^ random_1.fq.gz random_2.fq.gz > $@

$(REF1):
        curl -s "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/chromosomes/$(REF1).gz" |\
        gunzip -c | sed 's/^>chr/>/' | grep -v '[Nnatgc]' | head -n 100000 > $@


random_1.sai :  random_1.fq.gz $(REF2).bwt
        $(BWA) aln -f $@ $(REF2) $<

random_2.sai :  random_2.fq.gz $(REF2).bwt
        $(BWA) aln -f $@ $(REF2) $<

random_1.fq.gz random_2.fq.gz : $(REF1)
        $(SAMDIR)/misc/wgsim -N 1000000 -e 0.0 -r 0.0 -d 400 $< random_1.fq random_2.fq > wgsim.output
        gzip -f --best random_1.fq random_2.fq

$(REF2): $(REF1)
        cp $< $@
        sed 's/^>.*/>DUP/' $< >> $@

$(REF2).bwt : $(REF2)
        $(BWA) index -a bwtsw $<

$(REF2).fai :  $(REF2)
        $(SAMTOOLS) faidx $< 



clean:
        rm -f chr22.* *.bam *.vcf *.bcf *.sai *.gz *.fq *.bai  wgsim.output *.sam

That's it,

Pierre

The LZW compression algorithm as a measure of the short-reads complexity

The LZW algorithm, is a dictionary-based universal lossless data compression algorithm. The algorithm is easy to implement, here is a pseudocode (copied from there):

string s;
char ch;
...

s = empty string;
while (there is still data to be read)
{
    ch = read a character;
    if (dictionary contains s+ch)
    {
 s = s+ch;
    }
    else
    {
 encode s to output file;
 add s+ch to dictionary;
 s = ch;
    }
}
encode s to output file;

And here is my C++ implementation: the size of the dictionary reflects the complexity of the sequence: http://code.google.com/p/variationtoolkit/source/browse/trunk/src/lzw.h.

I've used this complexity to plot the number-of-reads=f(size-LZW);

Exome data

#complexity mapped unmapped sample
9 23274 21 TTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
10 1676 31379 CTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
11 2365 455 CCCTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
12 1523 5118 AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACGGAAAAAA
13 1941 2827 GTTCCTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTTT
14 2253 2495 CCCCCCCCCCCCCCCCCCCCCCCCCCACCCCCCCCCCCACCCCCCACCCACACC
15 2774 2908 ATAAAATAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGGAG
16 3965 3149 AAAAAAAACCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCTCCCCCCCCCTCCT
17 6944 4020 CTTTTTTTTTTTCTTTTCTTTTTTTTTTCCCTCTTTTTTTTTTTTTTTTTTTTC
18 11607 5143 TTGGTTTTTTTTTTTTTTTTTTTTTTTGGTTTGTTTTTTTTTTTTTTTTACCCT
19 19659 6724 GGGGGGGGGGGGGGGGGGGGGAGGAGGAAGGGGAGGAAGGGAGGAGGAAAGAGA
20 32504 9412 ACCCTAACCCTACCCCTAACCCTAACCCTAACCCAACCCTAACCCTAACCCTAA
21 50824 13984 TAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCC
22 77399 19651 CTAACCCTAACCCTAACCCTAACCCTAACCCTAACTCTAACCCTAACCCTAACC
23 114774 28966 GATCTCCCTAACCCTAACCCTACCCTAACCCTAACCCTAACCCTAACCCTAACC
24 176229 43729 TCCGATCTACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCC
25 316878 67402 TTCCGATCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAAC
26 721852 104378 TTCCGATCTGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTTAGGGTT
27 2028152 164968 CCGATCTAGGGTTAGGGTTAGGGTTGGGGTTAGGGTTAGGGTTAGGGTTAGGGT
28 6108817 284769 GCTGTGGTCTTCATCTGCAGGTGTCTGACTTCCAGCAACTGCTGGCCTGTGCCA
29 16907095 553236 GGGCACTGCAGGGCCCTCTTGCTTACTGTATAGTGGTGGCACGCCGCCCGCTGG
30 37103260 1130111 GTGATTTGGGCTGGGGCCTGGCCATGTGTATTTTTTTAAATTTCCACTGATGAT
31 55419720 1911772 GACCTGAGGAGAACTGTGCTCCGCCTTCAGAGTACCACCGAAATCTGTGCAGAG
32 47163550 2041580 GCCATGTGTATTTTTTTAAATTTCCACTGATGATTTTGCTGCATGGCCGGTGTT
33 17867328 1073476 CTGTATCCCACCAGCAATGTCTAGGAATGCCTGTTTCTCCACAAAGTGTTTACT
34 2014482 212089 TTTGCTGTCTCTTAGCCCAGACTTCCCGTGTCCTTTNNACCNGGCCTTTGAGAG
35 72637 30914 ACATCAANCTCAGGCACNTGGCCCAGGTCTGGCACTTAGAAGTAGTTCTCTGGG
36 8496 8247 AGGATATCTGGGNTGCNNCCGGAGTCGCAGTGTCTTGGGCCGCCTGAAGGTGAG
37 905 1506 AAGCATTACTGGAAACATCCTCATTGTGTTNTCTGNGACCANTNACCCTCACTN
38 58 102 TCGAGCNNCGTTGACTTCAGGNGGTCTNCTACCAGCAGCTCGNAATAGTTGCAC
39 1 0 AANTTCNAACGACTGTANNTCATNNGGCNNTGCNGGNCCNANAAACTGGCTGAG

Whole Genome data

#complexity mapped unmapped sample
13 1 2728 GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGG
14 0 608 GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGTGGGGGGGGGGG
15 0 2181 AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAANAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
16 7 2095 GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGTGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGTGGGGGGG
17 27 1924 AAAAAAAAAAAAAAAACTAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAT
18 41 2558 AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGGAAAAAATAAAAAAGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
19 66 2961 GGGGGGGGGGTGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGTGTGTGGGCG
20 127 3391 AGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGCGGGGGGGGGGGGCGGGGGGGGGAGGAGGGGGGGG
21 181 4244 NNNNNNNNNNNNNNNNNNNNNNNATTANNNNNNNNNNNNNNNNNNNTAANNNNNNNNNNNNNNNNNNNANNNNNNNNNNNANNNNNNNNNNNNNNNNNNN
22 371 5242 AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGGGGGGGAGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAGAGAAGAAAAAAAAAAA
23 627 6308 GGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGTGGGGGCGCGCGGCCGGGGGCGCGGGT
24 1398 8204 GGTATAATGCTAGGTATAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA
25 3990 10923 GAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGAAAGA
26 8469 15085 CATCAGAATACAGCTAACAAGAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACAAAAAAAAAAAAAAAAAAAAAA
27 14494 20162 GCGGTGGCGGGGGCCCGCGGGCCCCCCGCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGTGGGGGGGGGGGCGG
28 24273 26976 TTCTTTCTTTCTTTCTTTCTTTCTTTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTTCTTCCTCCTTTTCTTTCCTTTTCTTTCT
29 37918 35975 ACCAACACCACACCCCCACCCCCCCCCCCCCCCCACCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCACCCCCAACCCCTAACCCTAACCCTAACC
30 58261 48523 CCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCT
31 81164 62111 CTAACCCTAACCCTAACCCTAACCCTAACCCTCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCT
32 112886 79802 TAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCCAACCCCAACCCTAACCCCAAC
33 154666 101551 CCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCGACCCCTAACCCGA
34 204241 130861 CCCTAGCCCCTCCCTATCCCTAACCCTAACCCTAACCCTAAAACCCTAACCCTAAAACCCTAACCCTAAAACCCTAACCCTAACCCTAACCCAACCCTAA
35 267951 166738 TCACCCTCACCCTCACCCTCACCCTAACCCTCACCCTCACCCTCACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCTAACCCTAACC
36 347051 210144 TAACCCTAACCCTAACCCTAACCCCTAACCCTAACCCCAACCCTAACCCTAACCCTAACCCTAACCCCATCACTAACCTGTAACCCTCACCCTAACCCTA
37 447808 259191 CCCCACCCCCATCCCTAACCCGACCCTCAACCCAACCCCGAACCCAAACCCCAACCCCAACCCAAACCCAAACCCTAACCCTAACCCAAACCCTAACCCA
38 581941 320139 CTAACCCTAACCCTAACCCTAACCCTAACCCTTACCCTTACCCTTAACCCTCAACCCAACCCTAACACTAACCCTAACCCTAACCCCAAACCCAAGCCCA
39 770244 392189 CCCTACCCCTANCCCTACCCCTACCCCTAACCCTAACCCTAACCCTAACNCTAACCCAACCCCTCACACTACCCATCACCCCCACACCCTACCCCTACCC
40 1040306 485156 CCTTCACTCTACGCTTATCTCCCTACCTACCCCTAACCCTATCCCTAACCCTAACCCTATCCCTAACCCTAACCCTACCCCTAACCCTTACCCTAACCCA
41 1472994 592142 CAACCCGAGTACAATGGAAACGAATGGAATGGAATGAAATGGAATGGAATGGAATGGAATAGAATGGAATGGAATGGAATGGAATCAACCCGAGTGCAAT
42 2167116 711100 CCATCACCCCACCCCTACCCCTAACGCCACCCCTACCCCCAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCAGATCGGA
43 3307945 822818 CTTCCCCTACCCCCAACCCCGATCCCGAACCCAACCCCTAGCCCTACCCTTAACCCATCCCCATCCCTACCCCTAACCCTAACCCTAACCCTAAGCCAAC
44 5436093 924140 AGGGTTAGGGTAAGGGTTAGGGTTAGGGTTAGGGTTAGGGGTAGGGTTCGGGATTGGAAAGAGCGGCGGGTTGGGGGAGGGTTATGGGATCTGATGAAAT
45 10232300 1030747 CTAACCCTAACCCTAACCCTAACCTAACCCATCCCCCAGCCAACCTTTACCCTCACCCCTCCTCTGACCCTAACCCTCAACCTTCCCCTGCCTCGGAATC
46 21775389 1204827 TAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTCGCGGTACCCTCAGCCGGCCCGCCCGCCCGGGTCTGACCTGAGGAGAACTGTGCTCC
47 47888153 1605181 CTAACCCTAACCCTAACCCTAACCCTAACCTAAACCTTATCGCTATGCTTACCAGTAGCCTGAACCTGACCAATACACTAACCCTCACCCGGAAAATAAA
48 98581306 2536484 TAACCCTAACCCTAACCCTAACCCTAACCCTCGCGGTACCCTCAGCCGGCCCGCCCGCCCGGGTCTGACCTGAGGAGAACTGTGCTCCGCCTTCAGAGTA
49 175000036 4206340 GTGGTTTTTGTCTGCCAGTTCATGGTAATCACAGTGATTTCAAGGGGGGGTAAAAAAAGGAGGTGTGAGAGGGGCCCCCGGTTTCCACACAGACACCACA
50 246244603 6207853 CCTAACCCTAACCCTACCCATATACCTAACCCTAAAATTAAAAGTAATCATAACCCTAACCTTAGTTCTGCAACTACGGCTACACACACGTGCAGACCTA
51 247156482 7048827 CCTCTGGTGGCCCTGTCCGGGCATGACAGAAGGCGCGCACCCTTGACTTCTGTTCACTTCTCACTATGTCCCCTCAGCCCCTATCTCTGAATGGCCTGGC
52 155610209 5305604 GCGGTACCCTCAGCCGGCCCGCCCGCCCGGGTCTGACCTGAGGAGAACTGTGCTCCGCCTTCAGAGTACCAACGAAATCTGTGCAGAGGACAACGCAGCT
53 53933711 2288949 AGCGTCGCAACTCAAATGCAGCATTCCTAATGCACACATGACACCCAAAATATAACAGACATATTACTCATGGAGGGGGAGGGTGAGTGTGAGGGTGAGG
54 9490136 496289 TTTCACCAGAAGTAGGCCTCTTCCTGACAGGCAGCTGCACCACTGCCGGGCGCTGTGCCCTACCTTTGCTCTGCCCGCTGGAGACGGGGTTTGTCATGGG
55 991089 67105 AATTTCTGGAATGGATTATTAAACAGAGAGTCTGTAAGCACTTAGAAAAGGCCGCGGTGAGCCCCAGGGGCCAGCACTGCTCGAAATGTACAGCATTTCT
56 105830 11519 GGAAAATTTCTGGAATGGATTATTACAGAGTCTGTAAGCACTTAGAAAAGGCCGCGGTGAGTCCCAGGGGCCAGCACTGCTCGAAATGTACAGCATTTCT
57 21081 3128 TNACNGANGNNTNNNGTNTATTGNTCCAANAATCGNAGANNGAGAGGTTAAANTNNNNNNCNNNGATTNTGGGTTGTCTATTGATGTTTTTGGTCTATTC
58 6263 964 ATCNAGAGGCCAAGCCCAGCCTGTCNGCTTTNGTGTATAAAGNTCTCATGGAACAGAGCTGTGAGCCTGCCGNNTGTNGTCNNNNNNTNCGCCTGGNNAN
59 1360 246 ATTNGCCGGATGTGGTGGTGGGCGCCTGTAGTCCCAACTACTCAGGAGGCTGAAGCAGGAGAATGGCNAGAACNCGNNAGATGGNNGNTGNNGTNAGCCG
60 198 26 TCGGTCAACAAAATGGGTGACAGAGACCTACGCACGGATTATAATNNANCNGGCNCCANCCCGAGTGNTNNNCGGGGATTGGATGGNNCANNNTCCATAG
61 8 3 TGGAAAATNACTAGCNNGGAAGCAGACTNCGGGCCANANANATANNCAGTCACTTTANGCCCNGNANGGTGGNTCACANCTGTNATCCTANGNCNTTGGN
62 1 0 AGTTACGTGCTTACAGAATACTTTNTTTTGAGGTCAATANNANNANTAAGTNANGNATNCNNGATATCCTAGNGGGAATTCTCCGNCCTTCTGGAAGCTG

I'm sure there's must be something to say about this, but I just don't have time :-)

That's it,

Pierre

Memory-Mapping the Human Genome with 'mmap': my notebook

In this post, I've explored how to use a memory-mapped file to handle the 3Go of the Human Genome sequence as if it was entirely loaded in memory.
Via wikipedia: "A memory-mapped file is a segment of virtual memory which has been assigned a direct byte-for-byte correlation with some portion of a file or file-like resource. This resource is typically a file that is physically present on-disk... In computing, mmap is a POSIX-compliant Unix system call that maps files or devices into memory. It is a method of memory-mapped file I/O. It naturally implements demand paging, because initially file contents are not entirely read from disk and do not use physical RAM at all. The actual reads from disk are performed in "lazy" manner, after a specific location is accessed."
Using a C++ program, I'm going to memory-map the fasta sequence of the human genome (indexed with samtools faidx) and search the position of some short-reads using a BoyerMoore algorithm:

Required members:

/* maps a chromosome to its samtools faidx index */
map<string,faidx1_t> name2index;
/* used to get the size of the file */
struct stat buf;
/* genome fasta file file descriptor */
int fd;
/* the mmap (memory mapped) pointer */
char *mapptr;

Opening the mmap

string faidx(fastaFile);
string line;
faidx.append(".fai");
/* open *.fai file */
ifstream in(faidx.c_str(),ios::in);
/* read indexes in the .fai file that was created with samtools */
while(getline(in,line,'\n'))
 {
 faidx1_t index;
 //parse the faidx line...
 //...
 name2index.insert(make_pair(chrom,index));
 }
/* close index file */
in.close();

/* get the whole size of the fasta file */
stat(fasta, &buf);
/* open the fasta file */
fd = open(fastaFile, O_RDONLY);
/* open a memory mapped file associated to this fasta file descriptor */
mapptr = (char*)mmap(0, buf.st_size, PROT_READ, MAP_SHARED, fd, 0);

It's a kind of MAGIC: Getting the base at index 'position-th' of chromosome 'chrom'

std::map<std::string,faidx1_t>::iterator r=name2index.find(chrom);
faidx1_t& index=r->second;
char base=at(&index,position);
(...)
/* returns the base at position 'index' for the chromosome indexed by faidx */
char at(const FaidxPtr faidx,int64_t index)
 {
 long pos= faidx->offset +
  index / faidx->line_blen * faidx->line_len +
  index % faidx->line_blen
  ;
 /* here is the magic: no need to fseek/fread/ftell the file */
 return toupper(mapptr[pos]);
 }

Mapping the short reads

I've hacked a simple Boyer-Moore-Horspool algorithm from ttp://en.wikipedia.org/wiki/Boyer-Moore-Horspool_algorithm. Of course, you wouldn't use this algorithm to map your short reads for real :-) .

Disposing the mmap

/* close memory mapped map */
if(mapptr!=NULL) munmap(mapptr,buf.st_size);
 /* dispose fasta file descriptor */
if(fd!=-1) close(fd);

Compile & run

$  g++ -Wall testmmap.cpp -lz

$ ./a.out -g /path/tp/hg19.fa /path/to/my.fastq.gz

ATCATTTTCCTCCTAACAGATTAAAAATCAAGAAATATAAACCAGATGTAGCAG chr11 93065362
(...)

Source code

That's it,

Pierre

Visualizing my twitter network with Zoom.it

I wrote a small Java tool to download my twitter network as a GEXF file. This tool is available on github at:

https://github.com/lindenb/(...)/TwitterGraph.java

java -jar twittergraph.jar -o twittergraph.gexf 7431072 #my twitter ID

This tool doesn't use the OAuth API, so it have to wait for a few minutes, and retry to connect, every times it reaches the twitter API quotas (150 requests per hour). In the end it took one night to download the data from my network (~390 friends).

<gexf
xmlns="http://www.gexf.net/1.1draft"
xmlns:viz="http://www.gexf.net/1.1draft/viz"
xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
version="1.1"
xsi:schemaLocation="http://www.gexf.net/1.1draft http://www.gexf.net/1.1draft/gexf.xsd">
<meta lastmodifieddate="2011-02-04">
<creator>Gephi 0.7</creator>
<description/>
</meta>
<graph defaultedgetype="directed" timeformat="double" mode="dynamic">
<attributes class="node" mode="static">
<attribute id="name" title="name" type="string"/>
<attribute id="screenName" title="screenName" type="string"/>
<attribute id="imageUrl" title="imageUrl" type="string">
<default>http://a3.twimg.com/sticky/default_profile_images/default_profile_1_reasonably_small.png</default>
</attribute>
<attribute id="location" title="location" type="string"/>
<attribute id="description" title="description" type="string"/>
<attribute id="protectedProfile" title="protectedProfile" type="boolean"/>
<attribute id="friends" title="friends" type="integer"/>
<attribute id="followers" title="followers" type="integer"/>
<attribute id="listed" title="listed" type="integer"/>
<attribute id="utc_offset" title="utc offset" type="integer"/>
<attribute id="statuses_count" title="statuses count" type="integer"/>
</attributes>
<nodes>
<node id="6612402" label="sciencebase">
<attvalues>
<attvalue for="name" value="David Bradley"/>
<attvalue for="screenName" value="sciencebase"/>
<attvalue for="imageUrl" value="http://a3.twimg.com/profile_images/1142396198/twitter-blue-bradley_normal.jpg"/>
<attvalue for="location" value="Cambridge, UK"/>
<attvalue for="description" value="Science Writer David Bradley based in Cambridge, UK. Physical and life sciences news and views + technology, internet, web commentary."/>
<attvalue for="protectedProfile" value="false"/>
<attvalue for="friends" value="2022"/>
<attvalue for="followers" value="9197"/>
<attvalue for="listed" value="1065"/>
<attvalue for="utc_offset" value="0"/>
<attvalue for="statuses_count" value="7526"/>
</attvalues>
</node>
<node id="19344270" label="EMBOcomm">
<attvalues>
<attvalue for="name" value="Suzanne Beveridge"/>
<attvalue for="screenName" value="EMBOcomm"/>
<attvalue for="imageUrl" value="http://a0.twimg.com/profile_images/1189685782/S_Beveridge5100_normal.JPG"/>
<attvalue for="location" value="Heidelberg"/>
<attvalue for="description" value="Follow me for the latest from EMBO, the European Molecular Biology Organization"/>
<attvalue for="protectedProfile" value="false"/>
<attvalue for="friends" value="396"/>
<attvalue for="followers" value="697"/>
<attvalue for="listed" value="59"/>
<attvalue for="utc_offset" value="3600"/>
<attvalue for="statuses_count" value="632"/>
</attvalues>
</node>
<node id="20153702" label="walshtp">
<attvalues>
<attvalue for="name" value="Tom Walsh"/>
<attvalue for="screenName" value="walshtp"/>
<attvalue for="imageUrl" value="http://a3.twimg.com/profile_images/644287976/IMG_0815_normal.JPG"/>
<attvalue for="location" value="Dundee, Scotland"/>
<attvalue for="description" value="Scientific programmer and sysadmin. "/>
<attvalue for="protectedProfile" value="false"/>
<attvalue for="friends" value="129"/>
<attvalue for="followers" value="99"/>
<attvalue for="listed" value="8"/>
<attvalue for="utc_offset" value="0"/>
<attvalue for="statuses_count" value="783"/>
</attvalues>
</node>
<node id="15150655" label="konradfoerstner">
<attvalues>
<attvalue for="name" value="Konrad Förstner"/>
<attvalue for="screenName" value="konradfoerstner"/>
<attvalue for="imageUrl" value="http://a3.twimg.com/profile_images/643611092/konrad_avantar2_normal.jpeg"/>
<attvalue for="location" value="here and there"/>
<attvalue for="description" value="Idealist, Scientist, Includist, Data analyst, Open Source|Data|Access, Coder, Command line friend, CouchSurfer, Konrad"/>
<attvalue for="protectedProfile" value="false"/>
<attvalue for="friends" value="266"/>
<attvalue for="followers" value="167"/>
<attvalue for="listed" value="17"/>
<attvalue for="utc_offset" value="3600"/>
<attvalue for="statuses_count" value="1948"/>
</attvalues>
</node>

(...)

<edge id="E3811" source="14899756" target="14295341">
<attvalues>
<attvalue for="weight" value="1.0"/>
</attvalues>
</edge>
<edge id="E4816" source="14899756" target="19542750">
<attvalues>
<attvalue for="weight" value="1.0"/>
</attvalues>
</edge>
<edge id="E4830" source="14899756" target="60065276">
<attvalues>
<attvalue for="weight" value="1.0"/>
</attvalues>
</edge>
<edge id="E339" source="14899756" target="617133">
<attvalues>
<attvalue for="weight" value="1.0"/>
</attvalues>
</edge>
<edge id="E4807" source="14899756" target="15276911">
<attvalues>
<attvalue for="weight" value="1.0"/>
</attvalues>
</edge>
<edge id="E4822" source="14899756" target="26506721">
<attvalues>
<attvalue for="weight" value="1.0"/>
</attvalues>
</edge>
<edge id="E4824" source="14899756" target="27023131">
<attvalues>
<attvalue for="weight" value="1.0"/>
</attvalues>
</edge>
<edge id="E4819" source="14899756" target="22406785">
<attvalues>
<attvalue for="weight" value="1.0"/>
</attvalues>
</edge>
<edge id="E4808" source="14899756" target="16170580">
<attvalues>
<attvalue for="weight" value="1.0"/>
</attvalues>
</edge>
<edge id="E1237" source="14899756" target="4339911">
<attvalues>
<attvalue for="weight" value="1.0"/>
</attvalues>
</edge>
<edge id="E4828" source="14899756" target="56564230">
<attvalues>
<attvalue for="weight" value="1.0"/>
</attvalues>
</edge>
<edge id="E4826" source="14899756" target="33838201">
<attvalues>
<attvalue for="weight" value="1.0"/>
</attvalues>
</edge>
<edge id="E4815" source="14899756" target="19002481">
<attvalues>
<attvalue for="weight" value="1.0"/>
</attvalues>
</edge>
</edges>
</graph>
</gexf>

The GEXF file was then opened with Gephi, processed with the ForceAtlas algorithm and exported as a PDF file.

The PDF file was uploaded on scribd: http://www.scribd.com/doc/48415306/My-Twitter-Network

My Twitter Network

I then, downloaded the PDF from scribd.com, quickly copied the URL of the generated PDF and pasted it into http://zoom.it/.

Here is the result ! :-)



That's it !

Pierre

Playing with the Wordle algorithm: a tag cloud of Mesh Terms

The paper describing Wordle has been recently published (http://www.research.ibm.com/visual/papers/wordle_final2.pdf ). The algorithm was briefly described: “The most distinctive geometric aspect of a Wordle is the layout algorithm, which packs words to make efficient use of space. While many space-filling visualizations exist, they typically work by recursiveLayout proceeds according to this pseudocode:

sort words by weight, decreasing
for each word w:
     w.position := makeInitialPosition(w);
     while w intersects other words:
        updatePosition(w);

The two key procedures here are "makeInitialPosition" and "updatePosition". The makeInitialPosition routine picks a point at random according to a distribution that takes into account the desired overall shape, and, if desired, alphabetical order. The updatePosition routine moves the word on a spiral of increasing radius, radiating from the word's starting position. The updatePosition routine is aware of constraints on the overall shape of the Wordle. Constraining the layout to a rectangular shape causes updatePosition to prefer positions inside of the strict boundaries of the playing field; a blobby overall shape accepts boundary violations. The rectangular constraint is relaxed when the spiral radius exceeds either playing field dimension. ”

Your browser does not support the <CANVAS> element !

Just for fun , I've implemented my own version of the Wordle Alogrithm. The Java code is available on github at http://github.com/lindenb/jsandbox/blob/master/src/sandbox/MyWordle.java. I won't describe the program here, but I'll just say that the code invokes a java.awt.font.TextLayout class to get the shape of the text:

Graphics2D g=(...)
FontRenderContext frc = g.getFontRenderContext();
Font font=new Font("Dialog",Font.BOLD,fontSize);
TextLayout textLayout=new TextLayout(w.getText(), font, frc);
Shape shape=textLayout.getOutline(null);

and a java.awt.geom.Area to test if two shapes intersects.

Ok, let's test this code. FIrst I'm going to dump a pubmed query as XML with another simple tool named PubmedDump. This XML file is then parsed with a javascript program called from the SAX parser saxstream.jar (previously described here):
mesh.js

importPackage(Packages.sandbox);
importPackage(Packages.java.io);
importPackage(Packages.java.awt);
var content=null;
var mesh2count={};

function startElement(uri,localName,name,atts)
{
if(name=="DescriptorName")
{
content="";
}
}
function characters(s)
{
if(content!=null) content+=s;
}
function endElement(uri,localName,name)
{
if(content!=null)
{
var count=mesh2count[content];
if(count===undefined) count=0;
mesh2count[content]=count+1;
}
content=null;
}
function endDocument()
{
var w= new MyWordle();
for(var s in mesh2count)
{
var word= new MyWordle.Word(s,mesh2count[s]);
w.add(word);
}
w.setUseArea(true);/* use shape area instead of bounding boxes */
w.setAllowRotate(true);
w.setSortType(1);/* sort by weight */
w.doLayout();

var f=new File("result.svg");
w.saveAsSVG(f);
}

This javascript code counts the occurence of each MESH term (<DescriptorName>), and when the document has been parsed, a new instance of MyWordle class is created, filled and the result is saved to a SVG file.

Invocation

Here is an example for the query "Rotavirus NSP3 NSP1":

java -jar pubmeddump.jar "Rotavirus NSP3 NSP1" |\
java -cp mywordle.jar:saxscript.jar org.lindenb.tinytools.SAXScript -n -f mesh.js

Result

MyWordle.java

View more presentations from Pierre Lindenbaum.

That's it,
Pierre

Pairwise Alignments in C++ (for fun). My notebook.

With Next Generation Sequencing technologies, I have more and more opportunities to program in C++. C/C++ was my first programming language and I now enjoy my past experience with java for designing some classes, using the design patterns, etc... . Thus, I've recently played with the pairwise alignments in C++ and, in this post I'll describe the classes I've used (hey, I wrote this for fun, I didn't intend to write a formal and optimized code about this subject (I'm not a specialist of this anyway). Furthermore, I will not describe what is Dynamic Programming).

The Classes

Sequence

Sequence is an abstract sequence of characters (= string). It has a 'size' and should return the index-th character.

class Sequence
{
public:
Sequence() { }
virtual ~Sequence() { }
virtual size_type size() const=0;
virtual char at(size_type index) const=0;
}*SequencePtr;

CCSequence

CCSequence is a concrete implementation of Sequence using a std::string as the string holder.

class CCSequence:public Sequence
{
private:
std::string str;
public:
CCSequence(const std::string& str):Sequence(),str(str) { }
virtual ~CCSequence() { }
virtual size_type size() const
{
return (size_type)str.size();
}
virtual char at(size_type index) const
{
return str[index];
}
};

PtrSequence

PtrSequence is another implementation of Sequence but here, it uses a C pointer char* as the string holder.

class PtrSequence:public Sequence
{
private:
const char* str;
std::size_t length;
public:
PtrSequence(const char* str,size_t length):Sequence(),str(str),length(length) { }
PtrSequence(const char* str):str(str),length(0UL) { length=std::strlen(str); }
virtual ~PtrSequence() {}
virtual size_type size() const
{
return (size_type)this->length;
}
virtual char at(size_type index) const
{
return str[index];
}
};

PathItem

This class is a component of the pairwise alignment. It stores the positions on both sequences. A position equals to '-1' would be a gap.

typedef struct PathItem
{
Sequence::size_type x;
Sequence::size_type y;
PathItem(Sequence::size_type x,Sequence::size_type y):x(x),y(y)
{
}
Sequence::size_type get(int side) const { return (side==0?x:y);}
}*PathItemPtr;

Path

the class Path is a solution of the matrix traversal by the dynamic algorithm. It is a vector of PathItem and it also implements Sequence as the consensus of the alignment:

typedef class Path: public Sequence
{
private:
const SequencePtr seqX;
const SequencePtr seqY;
std::vector<PathItem> path;

std::auto_ptr<Sequence> _asSeq(int side) const
{
const SequencePtr seq=(side==0?seqX:seqY);
std::vector<PathItem>::const_iterator r=path.begin();
std::vector<PathItem>::const_iterator r_end=path.end();
std::ostringstream os;
while(r!=r_end)
{
Sequence::size_type n=r->get(side);
os << (n==-1?'-':seq->at(n));
++r;
}
return std::auto_ptr<Sequence>(new CCSequence(os.str()));
}
public:
Path(const SequencePtr seqX,const SequencePtr seqY,std::vector<PathItem> path):Sequence(),seqX(seqX),seqY(seqY),path(path)
{
}

virtual ~Path()
{
}

const SequencePtr X() const
{
return seqX;
}
const SequencePtr Y() const
{
return seqY;
}
std::auto_ptr<Sequence> consensusX() const
{
return _asSeq(0);
}
std::auto_ptr<Sequence> consensusY() const
{
return _asSeq(1);
}
std::auto_ptr<std::string> mid() const
{
std::vector<PathItem>::const_iterator r=path.begin();
std::vector<PathItem>::const_iterator r_end=path.end();
std::ostringstream os;
while(r!=r_end)
{
if(r->x==-1 || r->y==-1)
{
os << ' ';
}
else if(std::toupper(X()->at(r->x))==std::toupper(Y()->at(r->y)))
{
os << '|';
}
else
{
os << ' ';
}
++r;
}
return std::auto_ptr<std::string>(new std::string(os.str()));
}
virtual size_type size() const
{
return path.size();
}
virtual char at(size_type index) const
{
const PathItem& i=path.at(index);
if(i.x==-1 && i.y==-1) return '-';
if(i.x==-1) return Y()->at(i.y);
if(i.y==-1) return X()->at(i.x);
char cx= X()->at(i.x);
char cy= Y()->at(i.y);
return cx==cy?cx:'X';
}

std::ostream& print(std::ostream& out) const
{
out << *(consensusX()) << "\n"
<< *(mid()) << "\n"
<< *(consensusY()) << "\n"
;
return out;
}
}*PathPtr;

SubstitutionMatrix

SubstitutionMatrix is an abstract class. Its function compare returns the cost of the substitution of the symbol 'c1' by the symbol 'c2':

template<typename SCORE>
class SubstitutionMatrix
{
public:
virtual SCORE compare(char c1,char c2) const=0;
};

Identity

Identity is a simple implementation of SubstitutionMatrix. It returns 1 if the symbols are identical, else it returns -1.

template<typename SCORE>
class Identity:public SubstitutionMatrix<SCORE>
{
public:
virtual SCORE compare(char c1,char c2) const
{
return (SCORE)(std::toupper(c1)==std::toupper(c2)?1:-1);
}
};

Blosum

Blosum is an abstract implementation of SubstitutionMatrix for a BLOSUM matrix:

/** cf. ftp://ftp.ncbi.nih.gov/blast/matrices */
template<typename SCORE>
class Blosum:public SubstitutionMatrix<SCORE>
{
protected:
virtual int aa2index(char c) const
{
switch(std::toupper(c))
{
case 'A': return 0;
(...)
}
}
public:
virtual SCORE compare(char c1,char c2) const=0;
};

Blosum62

It is a concrete implementation of Blosum

static const signed char _blosum62[]={
4,-1,-2,-2,0,-1,-1,0,-2,-1,-1,-1,-1,-2,-1,1,0,-3,-2,0,-2,-1,0,-4,
(...)
};
/** ftp://ftp.ncbi.nih.gov/blast/matrices/BLOSUM62 */
template<typename SCORE>
class Blosum62:public Blosum<SCORE>
{
public:
virtual SCORE compare(char c1,char c2) const
{
return (SCORE)::_blosum62[this->aa2index(c1)*24+ this->aa2index(c2)];
}
};

Array2D

This abstract class holds the matrix of type 'T' for the dynamic algorithm. For a given matrix we need to know its width, its height , getthe value at position(x,y) and we want to set the value at position(x,y).

emplate<typename T>
class Array2D
{
public:
Array2D() {}
virtual ~Array2D() {}
virtual size_type width() const=0;
virtual size_type height() const=0;
virtual T get(size_type x,size_type y) const=0;
virtual void set(size_type x,size_type y,T val)=0;
protected:
std::size_t offset(size_type x,size_type y) const
{
return (this->width())*y + x;
}
};

DefaultArray2D

DefaultArray2D is a first concrete implementation of Array2D. It holds everything in memory:

template<typename T>
class DefaultArray2D:public Array2D<T>
{
protected:
size_type _width;
size_type _height;
T* matrix;

public:
DefaultArray2D( size_type width, size_type height):_width(width),_height(height)
{
matrix=new T[width*height];
}

virtual ~DefaultArray2D()
{
delete [] matrix;
}

virtual size_type width() const
{
return this->_width;
}

virtual size_type height() const
{
return this->_height;
}

virtual T get(size_type x,size_type y) const
{
return matrix[this->offset(x,y)];
}
virtual void set(size_type x,size_type y,T val)
{
matrix[this->offset(x,y)]=val;
}
};

StoredArray2D

StoredArray2D is second lousy concrete implementation of Array2D. It stores the matrix in a temporary file:

template<typename T>
class StoredArray2D:public Array2D<T>
{
private:
size_type _width;
size_type _height;
/** FILE pointer to the binary file */
mutable FILE* io_ptr;
void move(size_type x, size_type y) const
{
if(std::fseek (io_ptr,this->offset(x,y)*sizeof(T), SEEK_SET)!=0)
{
throw std::runtime_error("Cannot fseek");
}
}
public:
StoredArray2D( size_type width, size_type height):_width(width),_height(height),io_ptr(NULL)
{
io_ptr= std::tmpfile();
if(io_ptr==NULL) throw std::runtime_error("Cannot open tmp file");
T data;
std::memset(&data,0,sizeof(T));
for(std::size_t i=0;i< (size_t)(width*height);++i)
{
if(std::fwrite((void*)&data,sizeof(T),1,io_ptr)!=1)
{
std::fclose(io_ptr);
io_ptr=NULL;
throw std::runtime_error("write matrix");
}
}
}

virtual ~StoredArray2D()
{
if(io_ptr!=NULL) std::fclose(io_ptr);
}

virtual size_type width() const
{
return this->_width;
}

virtual size_type height() const
{
return this->_height;
}

virtual T get(size_type x,size_type y) const
{
T data;
move(x,y);
if(std::fread((void*)&data,sizeof(T),1,io_ptr)!=1)
{
throw std::runtime_error("cannot read");
}
return data;
}
virtual void set(size_type x,size_type y,T data)
{
move(x,y);
if(std::fwrite((void*)&data,sizeof(T),1,io_ptr)!=1)
{
throw std::runtime_error("cannot write");
}
}
};

Array2DFactory

Array2DFactory is a factory for an Array2D. It will allow our program to smoothly switch between any kind of Array2D:

template<typename T>
class Array2DFactory
{
private:
bool stored;
public:
//typedef typename Array2D<T>::size_type size_type;
Array2DFactory():stored(false) {}
virtual ~Array2DFactory() {}

void setStored(bool stored)
{
this->stored=stored;
}
bool isStored() const
{
return this->stored;
}

/** creates a new Array2D */
Array2D<T>* newArray(
size_type width,
size_type height
) const
{
if(isStored())
{
return new StoredArray2D<T>(width,height);
}
return new DefaultArray2D<T>(width,height);
}
};

Penaly

Penalty is an abstract class that will used to fill the edge of the matrix. It returns the score for inserting a symbol in a given sequence at a given position!

template<typename SCORE>
class Penalty
{
public:
Penalty() {}
virtual ~Penalty() {}
virtual SCORE get(const SequencePtr seq,int position) const=0;
};

DefaultPenalty

DefaultPenalty is a concrete implementation of Penalty. It uses a constant value for any position of the sequence.

template<typename SCORE>
class DefaultPenalty:public Penalty<SCORE>
{
private:
SCORE value;
public:
DefaultPenalty(SCORE value):value(value) {}
virtual ~DefaultPenalty() {}
virtual SCORE get(const SequencePtr seq,int position) const
{
return value;
}
};

Aligner

Aligner is an abstract pairwise aligner. It contains two Sequences, an Array2DFactory, a SubstitutionMatrix, etc... Its function align aligns the two sequences and path return the Path for the two aligned sequences:

template<typename T,typename SCORE>
class Aligner
{
private:
/** horizontal sequence */
SequencePtr seqX;
/** vertical sequence */
SequencePtr seqY;
/** array2D factory */
Array2DFactory<T>* array2dfactory;
/** substitution matrix */
SubstitutionMatrix<SCORE>* subsitutionMatrix;
/** penalty X */
Penalty<SCORE>* penaltyX;
/** penalty Y */
Penalty<SCORE>* penaltyY;
protected:
Aligner():seqX(NULL),seqY(NULL),
array2dfactory(NULL),
subsitutionMatrix(NULL)
{
}
public:
(...)
virtual SCORE align()=0;
virtual std::auto_ptr<Path> path()=0;
};

Needleman

Needleman is an implementation of Aligner for the Needleman & Wunsch algorithm.

template<typename T,typename SCORE>
class Needleman:public Aligner<T,SCORE>
{
private:
/** current matrix */
Array2D<T>* matrix;
public:
Needleman( ):matrix(NULL)
{

}

virtual ~Needleman()
{
}

/** clear the internal data if it exists */
virtual void clear()
{
Aligner<T,SCORE>::clear();
if(matrix!=NULL) delete matrix;
matrix=NULL;
}

virtual SCORE align()
{
/* clear the previous Array2D */
clear();
/* ask the factory for a new Array2D */
this->matrix = this->getArray2DFactory()->newArray(this->X()->size()+1, this->Y()->size()+1);

for(Sequence::size_type x=0;x<= this->X()->size();++x)
{
this->matrix->set(x,0,x*this->getPenaltyX()->get(this->X(),x) );
}
for(Sequence::size_type y=0;y<= this->Y()->size();++y)
{
this->matrix->set(0,y,y*this->getPenaltyY()->get(this->Y(),y));
}
for(Sequence::size_type x=1;x<= this->X()->size();++x)
{
for(Sequence::size_type y=1;y<= this->Y()->size();++y)
{
SCORE diag = matrix->get(x-1,y-1)+
this->compare(x-1,y-1);
SCORE delet= matrix->get(x-1,y) + this->getPenaltyX()->get(this->X(),x);
SCORE inser= matrix->get(x,y-1) + this->getPenaltyY()->get(this->Y(),y);
matrix->set(x,y,std::max(diag,std::max(delet,inser)));
}
}
return matrix->get(this->X()->size(),this->Y()->size());
}

public:

virtual std::auto_ptr<Path> path()
{
std::vector<PathItem> items;
Sequence::size_type x= this->X()->size();
Sequence::size_type y= this->Y()->size();
while (x>0 && y>0)
{
SCORE diag= matrix->get(x-1,y-1);
SCORE up = matrix->get(x, y- 1);
SCORE left = matrix->get(x-1,y);

if (diag >= up && diag >= left)
{
items.push_back(PathItem(x-1,y-1));
--x;
--y;
}
else if (left> up)
{
items.push_back(PathItem(x-1,-1));
--x;
}
else
{
items.push_back(PathItem(-1,y-1));
--y;
}
}

while (x > 0)
{
items.push_back(PathItem(x-1,-1));
--x;
}
while (y > 0)
{
items.push_back(PathItem(-1,y-1));
--y;
}
std::reverse(items.begin(),items.end());
return std::auto_ptr<Path>(new Path(this->X(),this->Y(),items));
}


};

SWaterman

SWaterman is an implementation of Aligner for the Smith & Waterman.

template<typename T,typename SCORE>
class SWaterman:public Aligner<T,SCORE>
{
private:
/** current matrix */
Array2D<T>* matrix;
SCORE best_score;
Sequence::size_type best_x;
Sequence::size_type best_y;
public:
SWaterman( ):matrix(NULL),best_score(0),best_x(0),best_y(0)
{

}

virtual ~SWaterman()
{
}

/** clear the internal data if it exists */
virtual void clear()
{
Aligner<T,SCORE>::clear();
if(matrix!=NULL) delete matrix;
matrix=NULL;
best_x=0;
best_y=0;
best_score=0;
}

virtual SCORE align()
{
/* clear the previous Array2D */
clear();
/* ask the factory for a new Array2D */
this->matrix = this->getArray2DFactory()->newArray(this->X()->size()+1, this->Y()->size()+1);

for(Sequence::size_type x=0;x<= this->X()->size();++x)
{
this->matrix->set(x,0,0);
}
for(Sequence::size_type y=0;y<= this->Y()->size();++y)
{
this->matrix->set(0,y,0);
}
for(Sequence::size_type x=1;x<= this->X()->size();++x)
{
for(Sequence::size_type y=1;y<= this->Y()->size();++y)
{
SCORE diag = matrix->get(x-1,y-1)+
this->compare(x-1,y-1);
SCORE delet= matrix->get(x-1,y) + this->getPenaltyX()->get(this->X(),x);
SCORE inser= matrix->get(x,y-1) + this->getPenaltyY()->get(this->Y(),y);

SCORE here= std::max((SCORE)0,std::max(diag,std::max(delet,inser)));
matrix->set(x,y,here);
if(best_score<here)
{
best_score=here;
best_x=x;
best_y=y;
}
}
}
return matrix->get(this->X()->size(),this->Y()->size());
}

public:

virtual std::auto_ptr<Path> path()
{
std::vector<PathItem> items;
Sequence::size_type x= best_x;
Sequence::size_type y= best_y;
while (x>0 && y>0 && matrix->get(x,y)!=0)
{
SCORE diag= matrix->get(x-1,y-1);
SCORE up = matrix->get(x, y- 1);
SCORE left = matrix->get(x-1,y);

if (diag >= up && diag >= left)
{
items.push_back(PathItem(x-1,y-1));
--x;
--y;
}
else if (left> up)
{
items.push_back(PathItem(x-1,-1));
--x;
}
else
{
items.push_back(PathItem(-1,y-1));
--y;
}
}

std::reverse(items.begin(),items.end());
return std::auto_ptr<Path>(new Path(this->X(),this->Y(),items));
}
};

AlignerFactory

AlignerFactory is a factory for an Aligner. It will allow our program to smoothly switch between any kind of Aligner:

template<typename T,typename SCORE>
class AlignerFactory
{
private:
bool local;
public:
AlignerFactory():local(false)
{
}

void setLocal(bool local) { this->local=local;}
bool isLocal() const { return local;}
std::auto_ptr<Aligner<T,SCORE> > newAligner()
{
if(isLocal())
{
return std::auto_ptr<Aligner<T,SCORE> >(new SWaterman<T,SCORE>());
}
return std::auto_ptr<Aligner<T,SCORE> >(new Needleman<T,SCORE>());
}
};

All in one, the full source code

#include <cstring>
#include <cstdlib>
#include <cstdio>
#include <cctype>
#include <algorithm>
#include <cassert>
#include <iostream>
#include <iomanip>
#include <stdexcept>
#include <vector>
#include <memory>
#include <sstream>

typedef signed long size_type;


typedef class Sequence
{
public:
typedef signed long size_type;
Sequence()
{
}
virtual ~Sequence()
{
}

virtual size_type size() const=0;
virtual char at(size_type index) const=0;
char operator[](size_type index) const
{
return this->at(index);
}
virtual std::ostream& print(std::ostream& out) const
{
size_type L=size();
for(size_type i=0;i< L;++i) out << at(i);
return out;
}
}*SequencePtr;

std::ostream& operator << (std::ostream& out, const Sequence& seq)
{
return seq.print(out);
}

class CCSequence:public Sequence
{
private:
std::string str;
public:
CCSequence(const std::string& str):Sequence(),str(str)
{
}
virtual ~CCSequence()
{
}
virtual size_type size() const
{
return (size_type)str.size();
}
virtual char at(size_type index) const
{
return str[index];
}
};

class PtrSequence:public Sequence
{
private:
const char* str;
std::size_t length;
public:
PtrSequence(const char* str,size_t length):Sequence(),str(str),length(length)
{
assert(str!=NULL);
assert(length>=0);
}
PtrSequence(const char* str):str(str),length(0UL)
{
assert(str!=NULL);
length=std::strlen(str);
}
virtual ~PtrSequence()
{
}
virtual size_type size() const
{
return (size_type)this->length;
}
virtual char at(size_type index) const
{
assert(index>=0);
assert(index<size());
return str[index];
}
};

typedef struct PathItem
{
Sequence::size_type x;
Sequence::size_type y;
PathItem(Sequence::size_type x,Sequence::size_type y):x(x),y(y)
{
}
Sequence::size_type get(int side) const { return (side==0?x:y);}
}*PathItemPtr;

typedef class Path: public Sequence
{
private:
const SequencePtr seqX;
const SequencePtr seqY;
std::vector<PathItem> path;

std::auto_ptr<Sequence> _asSeq(int side) const
{
const SequencePtr seq=(side==0?seqX:seqY);
std::vector<PathItem>::const_iterator r=path.begin();
std::vector<PathItem>::const_iterator r_end=path.end();
std::ostringstream os;
while(r!=r_end)
{
Sequence::size_type n=r->get(side);
os << (n==-1?'-':seq->at(n));
++r;
}
return std::auto_ptr<Sequence>(new CCSequence(os.str()));
}
public:
Path(const SequencePtr seqX,const SequencePtr seqY,std::vector<PathItem> path):Sequence(),seqX(seqX),seqY(seqY),path(path)
{
}

virtual ~Path()
{
}

const SequencePtr X() const
{
return seqX;
}
const SequencePtr Y() const
{
return seqY;
}
std::auto_ptr<Sequence> consensusX() const
{
return _asSeq(0);
}
std::auto_ptr<Sequence> consensusY() const
{
return _asSeq(1);
}
std::auto_ptr<std::string> mid() const
{
std::vector<PathItem>::const_iterator r=path.begin();
std::vector<PathItem>::const_iterator r_end=path.end();
std::ostringstream os;
while(r!=r_end)
{
if(r->x==-1 || r->y==-1)
{
os << ' ';
}
else if(std::toupper(X()->at(r->x))==std::toupper(Y()->at(r->y)))
{
os << '|';
}
else
{
os << ' ';
}
++r;
}
return std::auto_ptr<std::string>(new std::string(os.str()));
}
virtual size_type size() const
{
return path.size();
}
virtual char at(size_type index) const
{
const PathItem& i=path.at(index);
if(i.x==-1 && i.y==-1) return '-';
if(i.x==-1) return Y()->at(i.y);
if(i.y==-1) return X()->at(i.x);
char cx= X()->at(i.x);
char cy= Y()->at(i.y);
return cx==cy?cx:'X';
}

std::ostream& print(std::ostream& out) const
{
out << *(consensusX()) << "\n"
<< *(mid()) << "\n"
<< *(consensusY()) << "\n"
;
return out;
}
}*PathPtr;



std::ostream& operator << (std::ostream& out, const Path& path)
{
return path.print(out);
}


template<typename SCORE>
class SubstitutionMatrix
{
public:
virtual SCORE compare(char c1,char c2) const=0;
};



template<typename SCORE>
class Identity:public SubstitutionMatrix<SCORE>
{
public:
virtual SCORE compare(char c1,char c2) const
{
return (SCORE)(std::toupper(c1)==std::toupper(c2)?1:-1);
}
};

/** cf. ftp://ftp.ncbi.nih.gov/blast/matrices */
template<typename SCORE>
class Blosum:public SubstitutionMatrix<SCORE>
{
protected:
virtual int aa2index(char c) const
{
switch(std::toupper(c))
{
case 'A': return 0;
case 'R': return 1;
case 'N': return 2;
case 'D': return 3;
case 'C': return 4;
case 'Q': return 5;
case 'E': return 6;
case 'G': return 7;
case 'H': return 8;
case 'I': return 9;
case 'L': return 10;
case 'K': return 11;
case 'M': return 12;
case 'F': return 13;
case 'P': return 14;
case 'S': return 15;
case 'T': return 16;
case 'W': return 17;
case 'Y': return 18;
case 'V': return 19;
case 'B': return 20;
case 'Z': return 21;
case 'X': return 22;
default: return 23;
}
}
public:
virtual SCORE compare(char c1,char c2) const=0;
};


static const signed char _blosum62[]={
4,-1,-2,-2,0,-1,-1,0,-2,-1,-1,-1,-1,-2,-1,1,0,-3,-2,0,-2,-1,0,-4,
-1,5,0,-2,-3,1,0,-2,0,-3,-2,2,-1,-3,-2,-1,-1,-3,-2,-3,-1,0,-1,-4,
-2,0,6,1,-3,0,0,0,1,-3,-3,0,-2,-3,-2,1,0,-4,-2,-3,3,0,-1,-4,
-2,-2,1,6,-3,0,2,-1,-1,-3,-4,-1,-3,-3,-1,0,-1,-4,-3,-3,4,1,-1,-4,
0,-3,-3,-3,9,-3,-4,-3,-3,-1,-1,-3,-1,-2,-3,-1,-1,-2,-2,-1,-3,-3,-2,-4,
-1,1,0,0,-3,5,2,-2,0,-3,-2,1,0,-3,-1,0,-1,-2,-1,-2,0,3,-1,-4,
-1,0,0,2,-4,2,5,-2,0,-3,-3,1,-2,-3,-1,0,-1,-3,-2,-2,1,4,-1,-4,
0,-2,0,-1,-3,-2,-2,6,-2,-4,-4,-2,-3,-3,-2,0,-2,-2,-3,-3,-1,-2,-1,-4,
-2,0,1,-1,-3,0,0,-2,8,-3,-3,-1,-2,-1,-2,-1,-2,-2,2,-3,0,0,-1,-4,
-1,-3,-3,-3,-1,-3,-3,-4,-3,4,2,-3,1,0,-3,-2,-1,-3,-1,3,-3,-3,-1,-4,
-1,-2,-3,-4,-1,-2,-3,-4,-3,2,4,-2,2,0,-3,-2,-1,-2,-1,1,-4,-3,-1,-4,
-1,2,0,-1,-3,1,1,-2,-1,-3,-2,5,-1,-3,-1,0,-1,-3,-2,-2,0,1,-1,-4,
-1,-1,-2,-3,-1,0,-2,-3,-2,1,2,-1,5,0,-2,-1,-1,-1,-1,1,-3,-1,-1,-4,
-2,-3,-3,-3,-2,-3,-3,-3,-1,0,0,-3,0,6,-4,-2,-2,1,3,-1,-3,-3,-1,-4,
-1,-2,-2,-1,-3,-1,-1,-2,-2,-3,-3,-1,-2,-4,7,-1,-1,-4,-3,-2,-2,-1,-2,-4,
1,-1,1,0,-1,0,0,0,-1,-2,-2,0,-1,-2,-1,4,1,-3,-2,-2,0,0,0,-4,
0,-1,0,-1,-1,-1,-1,-2,-2,-1,-1,-1,-1,-2,-1,1,5,-2,-2,0,-1,-1,0,-4,
-3,-3,-4,-4,-2,-2,-3,-2,-2,-3,-2,-3,-1,1,-4,-3,-2,11,2,-3,-4,-3,-2,-4,
-2,-2,-2,-3,-2,-1,-2,-3,2,-1,-1,-2,-1,3,-3,-2,-2,2,7,-1,-3,-2,-1,-4,
0,-3,-3,-3,-1,-2,-2,-3,-3,3,1,-2,1,-1,-2,-2,0,-3,-1,4,-3,-2,-1,-4,
-2,-1,3,4,-3,0,1,-1,0,-3,-4,0,-3,-3,-2,0,-1,-4,-3,-3,4,1,-1,-4,
-1,0,0,1,-3,3,4,-2,0,-3,-3,1,-1,-3,-1,0,-1,-3,-2,-2,1,4,-1,-4,
0,-1,-1,-1,-2,-1,-1,-1,-1,-1,-1,-1,-1,-1,-2,0,0,-2,-1,-1,-1,-1,-1,-4,
-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,-4,1
};

/** ftp://ftp.ncbi.nih.gov/blast/matrices/BLOSUM62 */
template<typename SCORE>
class Blosum62:public Blosum<SCORE>
{
private:

const static signed char blosum62[];
public:
virtual SCORE compare(char c1,char c2) const
{
return (SCORE)::_blosum62[this->aa2index(c1)*24+ this->aa2index(c2)];
}
};




template<typename T>
class Array2D
{
public:
//typedef unsigned int size_type;
Array2D() {}
virtual ~Array2D() {}

virtual size_type width() const=0;
virtual size_type height() const=0;
virtual T get(size_type x,size_type y) const=0;
virtual void set(size_type x,size_type y,T val)=0;
virtual std::ostream& print(std::ostream& out) const
{
for(size_type j=0;j< height();++j)
{
for(size_type i=0;i< width();++i)
{
if(i>0) out << " ";
out << std::setw(4) << get(i,j);
}
out << std::endl;
}
return out;
}
protected:
std::size_t offset(size_type x,size_type y) const
{
return (this->width())*y + x;
}
};


template<typename T>
class DefaultArray2D:public Array2D<T>
{
protected:
size_type _width;
size_type _height;
T* matrix;

public:
DefaultArray2D( size_type width, size_type height):_width(width),_height(height)
{
matrix=new T[width*height];
}

virtual ~DefaultArray2D()
{
delete [] matrix;
}

virtual size_type width() const
{
return this->_width;
}

virtual size_type height() const
{
return this->_height;
}

virtual T get(size_type x,size_type y) const
{
return matrix[this->offset(x,y)];
}
virtual void set(size_type x,size_type y,T val)
{
matrix[this->offset(x,y)]=val;
}
};

template<typename T>
class StoredArray2D:public Array2D<T>
{
private:
size_type _width;
size_type _height;
/** FILE pointer to the binary file */
mutable FILE* io_ptr;
void move(size_type x, size_type y) const
{
if(std::fseek (io_ptr,this->offset(x,y)*sizeof(T), SEEK_SET)!=0)
{
throw std::runtime_error("Cannot fseek");
}
}
public:
StoredArray2D( size_type width, size_type height):_width(width),_height(height),io_ptr(NULL)
{
io_ptr= std::tmpfile();
if(io_ptr==NULL) throw std::runtime_error("Cannot open tmp file");
T data;
std::memset(&data,0,sizeof(T));
for(std::size_t i=0;i< (size_t)(width*height);++i)
{
if(std::fwrite((void*)&data,sizeof(T),1,io_ptr)!=1)
{
std::fclose(io_ptr);
io_ptr=NULL;
throw std::runtime_error("write matrix");
}
}
}

virtual ~StoredArray2D()
{
if(io_ptr!=NULL) std::fclose(io_ptr);
}

virtual size_type width() const
{
return this->_width;
}

virtual size_type height() const
{
return this->_height;
}

virtual T get(size_type x,size_type y) const
{
T data;
move(x,y);
if(std::fread((void*)&data,sizeof(T),1,io_ptr)!=1)
{
throw std::runtime_error("cannot read");
}
return data;
}
virtual void set(size_type x,size_type y,T data)
{
move(x,y);
if(std::fwrite((void*)&data,sizeof(T),1,io_ptr)!=1)
{
throw std::runtime_error("cannot write");
}
}
};


template<typename T>
class Array2DFactory
{
private:
bool stored;
public:
//typedef typename Array2D<T>::size_type size_type;
Array2DFactory():stored(false) {}
virtual ~Array2DFactory() {}

void setStored(bool stored)
{
this->stored=stored;
}
bool isStored() const
{
return this->stored;
}

/** creates a new Array2D */
Array2D<T>* newArray(
size_type width,
size_type height
) const
{
if(isStored())
{
return new StoredArray2D<T>(width,height);
}
return new DefaultArray2D<T>(width,height);
}
};


template<typename SCORE>
class Penalty
{
public:
Penalty() {}
virtual ~Penalty() {}
virtual SCORE get(const SequencePtr seq,int position) const=0;
};

template<typename SCORE>
class DefaultPenalty:public Penalty<SCORE>
{
private:
SCORE value;
public:
DefaultPenalty(SCORE value):value(value) {}
virtual ~DefaultPenalty() {}
virtual SCORE get(const SequencePtr seq,int position) const
{
return value;
}
};

template<typename T,typename SCORE>
class Aligner
{
private:
/** horizontal sequence */
SequencePtr seqX;
/** vertical sequence */
SequencePtr seqY;
/** array2D factory */
Array2DFactory<T>* array2dfactory;
/** substitution matrix */
SubstitutionMatrix<SCORE>* subsitutionMatrix;
/** penalty X */
Penalty<SCORE>* penaltyX;
/** penalty Y */
Penalty<SCORE>* penaltyY;
protected:


Aligner():seqX(NULL),seqY(NULL),
array2dfactory(NULL),
subsitutionMatrix(NULL)
{
}
public:
/** clear the internal data if it exists */
virtual void clear()
{
}

virtual ~Aligner()
{
clear();
}
/** returns the horizontal sequence */
virtual const SequencePtr X() const
{
return this->seqX;
}
/** returns the vertical sequence */
virtual const SequencePtr Y() const
{
return this->seqY;
}
virtual void setX(SequencePtr seqX) { this->seqX=seqX; clear();}
virtual void setY(SequencePtr seqY) { this->seqY=seqY; clear();}
virtual void setArray2DFactory(Array2DFactory<T>* array2dfactory) { this->array2dfactory=array2dfactory; clear();}
virtual const Array2DFactory<T>* getArray2DFactory() const { return this->array2dfactory;}
virtual void setSubstitutionMatrix(SubstitutionMatrix<SCORE>* subsitutionMatrix) { this->subsitutionMatrix=subsitutionMatrix; clear();}
virtual const SubstitutionMatrix<SCORE>* getSubstitutionMatrix() const { return this->subsitutionMatrix;}
virtual void setPenaltyX(Penalty<SCORE>* penaltyX) { this->penaltyX=penaltyX; clear();}
virtual Penalty<SCORE>* getPenaltyX() const { return this->penaltyX;}
virtual void setPenaltyY(Penalty<SCORE>* penaltyY) { this->penaltyY=penaltyY; clear();}
virtual Penalty<SCORE>* getPenaltyY() const { return this->penaltyY;}

virtual SCORE align()=0;
virtual std::auto_ptr<Path> path()=0;
protected:
/** shortcut to compare X[x] and Y[y] */
int compare(Sequence::size_type posx,Sequence::size_type posy) const
{
return getSubstitutionMatrix()->compare(X()->at(posx),Y()->at(posy));
}
};

template<typename T,typename SCORE>
class Needleman:public Aligner<T,SCORE>
{
private:
/** current matrix */
Array2D<T>* matrix;
public:
Needleman( ):matrix(NULL)
{

}

virtual ~Needleman()
{
}

/** clear the internal data if it exists */
virtual void clear()
{
Aligner<T,SCORE>::clear();
if(matrix!=NULL) delete matrix;
matrix=NULL;
}

virtual SCORE align()
{
/* clear the previous Array2D */
clear();
/* ask the factory for a new Array2D */
this->matrix = this->getArray2DFactory()->newArray(this->X()->size()+1, this->Y()->size()+1);

for(Sequence::size_type x=0;x<= this->X()->size();++x)
{
this->matrix->set(x,0,x*this->getPenaltyX()->get(this->X(),x) );
}
for(Sequence::size_type y=0;y<= this->Y()->size();++y)
{
this->matrix->set(0,y,y*this->getPenaltyY()->get(this->Y(),y));
}
for(Sequence::size_type x=1;x<= this->X()->size();++x)
{
for(Sequence::size_type y=1;y<= this->Y()->size();++y)
{
SCORE diag = matrix->get(x-1,y-1)+
this->compare(x-1,y-1);
SCORE delet= matrix->get(x-1,y) + this->getPenaltyX()->get(this->X(),x);
SCORE inser= matrix->get(x,y-1) + this->getPenaltyY()->get(this->Y(),y);
matrix->set(x,y,std::max(diag,std::max(delet,inser)));
}
}
return matrix->get(this->X()->size(),this->Y()->size());
}

public:

virtual std::auto_ptr<Path> path()
{
std::vector<PathItem> items;
Sequence::size_type x= this->X()->size();
Sequence::size_type y= this->Y()->size();
while (x>0 && y>0)
{
SCORE diag= matrix->get(x-1,y-1);
SCORE up = matrix->get(x, y- 1);
SCORE left = matrix->get(x-1,y);

if (diag >= up && diag >= left)
{
items.push_back(PathItem(x-1,y-1));
--x;
--y;
}
else if (left> up)
{
items.push_back(PathItem(x-1,-1));
--x;
}
else
{
items.push_back(PathItem(-1,y-1));
--y;
}
}

while (x > 0)
{
items.push_back(PathItem(x-1,-1));
--x;
}
while (y > 0)
{
items.push_back(PathItem(-1,y-1));
--y;
}
std::reverse(items.begin(),items.end());
return std::auto_ptr<Path>(new Path(this->X(),this->Y(),items));
}


};


template<typename T,typename SCORE>
class SWaterman:public Aligner<T,SCORE>
{
private:
/** current matrix */
Array2D<T>* matrix;
SCORE best_score;
Sequence::size_type best_x;
Sequence::size_type best_y;
public:
SWaterman( ):matrix(NULL),best_score(0),best_x(0),best_y(0)
{

}

virtual ~SWaterman()
{
}

/** clear the internal data if it exists */
virtual void clear()
{
Aligner<T,SCORE>::clear();
if(matrix!=NULL) delete matrix;
matrix=NULL;
best_x=0;
best_y=0;
best_score=0;
}

virtual SCORE align()
{
/* clear the previous Array2D */
clear();
/* ask the factory for a new Array2D */
this->matrix = this->getArray2DFactory()->newArray(this->X()->size()+1, this->Y()->size()+1);

for(Sequence::size_type x=0;x<= this->X()->size();++x)
{
this->matrix->set(x,0,0);
}
for(Sequence::size_type y=0;y<= this->Y()->size();++y)
{
this->matrix->set(0,y,0);
}
for(Sequence::size_type x=1;x<= this->X()->size();++x)
{
for(Sequence::size_type y=1;y<= this->Y()->size();++y)
{
SCORE diag = matrix->get(x-1,y-1)+
this->compare(x-1,y-1);
SCORE delet= matrix->get(x-1,y) + this->getPenaltyX()->get(this->X(),x);
SCORE inser= matrix->get(x,y-1) + this->getPenaltyY()->get(this->Y(),y);

SCORE here= std::max((SCORE)0,std::max(diag,std::max(delet,inser)));
matrix->set(x,y,here);
if(best_score<here)
{
best_score=here;
best_x=x;
best_y=y;
}
}
}
return matrix->get(this->X()->size(),this->Y()->size());
}

public:

virtual std::auto_ptr<Path> path()
{
std::vector<PathItem> items;
Sequence::size_type x= best_x;
Sequence::size_type y= best_y;
while (x>0 && y>0 && matrix->get(x,y)!=0)
{
SCORE diag= matrix->get(x-1,y-1);
SCORE up = matrix->get(x, y- 1);
SCORE left = matrix->get(x-1,y);

if (diag >= up && diag >= left)
{
items.push_back(PathItem(x-1,y-1));
--x;
--y;
}
else if (left> up)
{
items.push_back(PathItem(x-1,-1));
--x;
}
else
{
items.push_back(PathItem(-1,y-1));
--y;
}
}

std::reverse(items.begin(),items.end());
//matrix->print(std::cout);
//std::cout << best_score << "=" << best_x<<"," << best_y << std::endl;
return std::auto_ptr<Path>(new Path(this->X(),this->Y(),items));
}


};

template<typename T,typename SCORE>
class AlignerFactory
{
private:
bool local;
public:
AlignerFactory():local(false)
{
}

void setLocal(bool local) { this->local=local;}
bool isLocal() const { return local;}
std::auto_ptr<Aligner<T,SCORE> > newAligner()
{
if(isLocal())
{
return std::auto_ptr<Aligner<T,SCORE> >(new SWaterman<T,SCORE>());
}
return std::auto_ptr<Aligner<T,SCORE> >(new Needleman<T,SCORE>());
}
};




typedef signed long penalty_t;

/**
*
* main
*
*/
int main(int argc,char **argv)
{
try
{
int optind=1;
bool align_local=false;
bool stored_matrix=false;
bool use_blosum62=true;
penalty_t gX=-5;
penalty_t gY=-5;
while(optind < argc)
{
if(std::strcmp(argv[optind],"-h")==0)
{
std::cout << argv[0] << ". Pierre Lindenbaum PhD. Compiled on " << __DATE__ << " at " << __TIME__ << "\n" <<
"options\n" <<
" -h Help (this screen)\n" <<
" -L local-alignement\n" <<
" -s stored matrix (slower)\n" <<
" -i identity matrix (instead of blosum62)\n" <<
" -x penalty X (" <<gX <<")\n" <<
" -y penalty Y (" <<gY <<")\n" <<
"seq1 seq2\n" <<
std::endl;
std::exit(EXIT_FAILURE);
}
else if(std::strcmp(argv[optind],"-L")==0)
{
align_local=true;
}
else if(std::strcmp(argv[optind],"-s")==0)
{
stored_matrix=true;
}
else if(std::strcmp(argv[optind],"-i")==0)
{
use_blosum62=false;
}
else if(std::strcmp(argv[optind],"-x")==0)
{
gX=(penalty_t)atol(argv[++optind]);
}
else if(std::strcmp(argv[optind],"-y")==0)
{
gY=(penalty_t)atol(argv[++optind]);
}
else if(std::strcmp(argv[optind],"--")==0)
{
++optind;
break;
}
else if(argv[optind][0]=='-')
{
std::cerr << "unknown option " << argv[optind] << std::endl;
std::exit(EXIT_FAILURE);
}
else
{
break;
}
++optind;
}
if(optind+2!=argc)
{
std::cerr << "expected only two args"<< std::endl;
std::exit(EXIT_FAILURE);
}
Blosum62<penalty_t> blosum;
Identity<penalty_t> identity;


Array2DFactory<int> factory;
factory.setStored(stored_matrix);

DefaultPenalty<penalty_t> penalty_x(gX);
DefaultPenalty<penalty_t> penalty_y(gY);
AlignerFactory<int,penalty_t> algofactory;
algofactory.setLocal(align_local);

std::auto_ptr<Aligner<int,penalty_t> > algo = algofactory.newAligner();

algo->setArray2DFactory(&factory);
algo->setSubstitutionMatrix(use_blosum62?(SubstitutionMatrix<penalty_t>*) &blosum :(SubstitutionMatrix<penalty_t>*) &identity);
algo->setPenaltyX(&penalty_x);
algo->setPenaltyY(&penalty_y);


PtrSequence seq1(argv[optind++]);
PtrSequence seq2(argv[optind++]);

algo->setX(&seq1);
algo->setY(&seq2);
algo->align();

std::cout << (*(algo->path())) << std::endl;

}
catch(std::exception& err)
{
std::cerr << err.what() << std::endl;
}
catch(...)
{
std::cerr << "BOUM" << std::endl;
}
return EXIT_SUCCESS;
}

Compilation

g++ -Wall -O3 needle.cpp

Examples

./a.out -h
./a.out. Pierre Lindenbaum PhD. Compiled on Jul 8 2010 at 23:03:05
options
-h Help (this screen)
-L local-alignement
-s stored matrix (slower)
-i identity matrix (instead of blosum62)
-x penalty X (-5)
-y penalty Y (-5)
seq1 seq2

time ./a.out MERQKRKADIEKGLQFIQSTLPLKQEEYEAFLLKLVQNLFAEGN MSSVSEERRKRQQNIKEGLQFIQSPLSYPGTQEQYAVYLHALVR
MERQK-----RKADIEKGLQFIQSTLP--LKQEEYEAFLLKLVQNLFAEGN
| | | ||||||| | || | | ||
MSSVSEERRKRQQNIKEGLQFIQSPLSYPGTQEQYAVYLHALVR-------


real 0m0.006s
user 0m0.004s
sys 0m0.004s

time ./a.out -s MERQKRKADIEKGLQFIQSTLPLKQEEYEAFLLKLVQNLFAEGN MSSVSEERRKRQQNIKEGLQFIQSPLSYPGTQEQYAVYLHALVR

MERQK-----RKADIEKGLQFIQSTLP--LKQEEYEAFLLKLVQNLFAEGN
| | | ||||||| | || | | ||
MSSVSEERRKRQQNIKEGLQFIQSPLSYPGTQEQYAVYLHALVR-------


real 0m0.037s
user 0m0.004s
sys 0m0.032s

./a.out -i MERQKRKADIEKGLQFIQSTLPLKQEEYEAFLLKLVQNLFAEGN MSSVSEERRKRQQNIKEGLQFIQSPLSYPGTQEQYAVYLHALVR

MERQKRKADIEKGLQFIQSTLPLKQEEYEAFLLKLVQNLFAEGN
| | | |
MSSVSEERRKRQQNIKEGLQFIQSPLSYPGTQEQYAVYLHALVR

./a.out -L MERQKRKADIEKGLQFIQSTLPLKQEEYEAFLLKLVQNLFAEGN MSSVSEERRKRQQNIKEGLQFIQSPLSYPGTQEQYAVYLHALVR

ME--RQKRKADIEKGLQFIQSTLP--LKQEEYEAFLLKLVQ
| || | ||||||| | || | | ||
VSEERRKRQQNIKEGLQFIQSPLSYPGTQEQYAVYLHALVR

./a.out -x -5 -y -20 -L MERQKRKADIEKGLQFIQSTLPLKQEEYEAFLLKLVQNLFAEGN MSSVSEERRKRQQNIKEGLQFIQSPLSYPGTQEQYAVYLHALVR

ME--RQKRKADIEKGLQFIQSTL
| || | ||||||| |
VSEERRKRQQNIKEGLQFIQSPL

That's it

Pierre