Drawing a SVG timeline with the http://data.bnf.fr data

The National French Library has recently started to release
its data as RDF/XML. Here, I've played with the biographies of the famous French writers to create a simple timeline.

Unfortunately, this timeline is too large to be displayed here. :-)

However, here is the java code I used to generate this map as a SVG document:


See also: Freebase and the History of Sciences

That's it,

Pierre

An interactive 'dialog' generating a SAM Flag

The following shell script uses the dialog utility to generate a numeric Sam flag.

$ sh selectsam.sh 
      ┌───────────────────────SAM FLAGS──────────────────────────┐
      │ SELECT FLAGS                                             │  
      │ ┌──────────────────────────────────────────────────────┐ │  
      │ │ [X] 1     Read Paired                                │ │  
      │ │ [ ] 2     Read mapped in proper pair                 │ │  
      │ │ [X] 4     Read unmapped                              │ │  
      │ │ [ ] 8     Mate unmapped                              │ │  
      │ │ [ ] 16    Read reverse strand                        │ │  
      │ │ [X] 32    Mate reverse strand                        │ │  
      │ │ [ ] 64    First in pair                              │ │  
      │ │ [X] 128   Second in pair                             │ │  
      │ │ [ ] 256   Not primary alignment                      │ │  
      │ │ [ ] 512   Read fails platform/vendor quality checks  │ │  
      │ │ [ ] 1024  Read is PCR or optical duplicate           │ │  
      │ └─────↓(+)─────────────────────────────────────────────┘ │  
      │                                                          │  
      │                                                          │  
      ├──────────────────────────────────────────────────────────┤  
      │                         <  OK  >                         │  
      └──────────────────────────────────────────────────────────┘  
        
$ 165

Source available at: https://github.com/lindenb/samtools-utilities/blob/master/script/selectflag.sh

That's it,

Pierre

I can't save you, but I can help you: reverse-complementing a DNA in OpenOffice

People in my lab are using "Microsoft Word" to edit and annotate their sequences. Just like Neil, I can't save them but I can help them: here is an OpenOffice / LibreOffice macro reverse-complementing the current selected text.

• open libreoffice/openoffice
• Open menu "Tools / Macros / Organize macros / Basic...
• Select "My Macros/Standard" and click on "New"
• Copy+paste my macro here
• Click on the top button "compile"
• Close the dialog
• Menu "Tools/Customize" , "Keyboard" tab, "Toolbar content" . Button "Add... " / Category "Libre Office Macros" . Select "My Macro/Standard/Modulexxx"/ReverseComplementSelection". Click on "Modify" to assign an icon to this new button
• There is now a new button in the toolbar. Clicking on that button reverse-complements the selected text.

That's it,

Pierre

PS: Hum ? what did you say, they use Microsoft Word ? Not OpenOffice ? pfff....

A text alignment viewer using the samtools API

I wrote a modified version of samtools tview. The original code uses the curses API to display an interactive sequence viewer.

I wanted to generate this kind of screen for a large number of positions and I wanted to be able to redirect the output to a unix pipeline so I changed the original code in order to only handle an extensible 2-dimensional array of characters. The code I wrote is available on github at:

https://github.com/lindenb/samtools-utilities/blob/master/src/ttview.c

.

Compilation

cd samtools/
make #compile samtools
gcc -o bamttview -g -Wall -O2 -DSTANDALONE_VERSION  -I. -Lbcftools  bam_ttview.c  bam2bcf.o   errmod.o  bam_color.o libbam.a -lbcf  -lm -lz

Example

print a position

$ ./bamttview  -g "ref:5" examples/toy.bam examples/toy.fa |\
cat -n
     1       11              21        31         41        51        61          
     2 TGTTAGATAA****GATA**GCTGTGCTAGTAGGCAG*TCAGCGCCATNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
     3   ........    ....  ......K.K......K. ..........                                
     4   ........AGAG....***...      ,,,,,    ,,,,,,,,,                                
     5     ......GG**....AA                                                            
     6     ..C...**** ...**...>>>>>>>>>>>>>>T.....                                     

print a list of positions:

$ cat positions.txt 
ref2:10-100
ref:10-15
ref2:11

$ ./bamttview  -f positions.txt examples/toy.bam examples/toy.fa |\
cat -n | head -n 25
     1 
     2 
     3 > ref2:10-100
     4 
     5  11            21        31        41        51        61        71             
     6 aaaac****aattaagtctacagagcaactaNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
     7 ....Y    ..W...................                                                 
     8 .....****..A...                                                                 
     9 .....****..A...T.                                                               
    10 .....AAAT.............                                                          
    11 C...T****....................                                                   
    12   ..T****.....................                                                  
    13     T****......................                                                 
    14 
    15 
    16 
    17 > ref:10-15
    18 
    19  11              21        31         41        51        61        71          
    20 GATAA****GATA**GCTGTGCTAGTAGGCAG*TCAGCGCCATNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
    21 .....    ....  ......K.K......K. ..........                                     
    22 .....AGAG....***...      ,,,,,    ,,,,,,,,,                                     
    23 .....GG**....AA                                                                 
    24 .C...**** ...**...>>>>>>>>>>>>>>T.....                                          
    25 

./bamttview  -f positions.txt -d examples/toy.bam examples/toy.fa |\
cat -n | head -n 25
     1 
     2 
     3 > ref2:10-100
     4 
     5  11            21        31        41        51        61        71             
     6 aaaac****aattaagtctacagagcaactaNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
     7 ....Y    ..W...................                                                 
     8 AAAAC****AAATAA                                                                 
     9 AAAAC****AAATAATT                                                               
    10 AAAACAAATAATTAAGTCTACA                                                          
    11 CAAAT****AATTAAGTCTACAGAGCAAC                                                   
    12   AAT****AATTAAGTCTACAGAGCAACT                                                  
    13     T****AATTAAGTCTACAGAGCAACTA                                                 
    14 
    15 
    16 
    17 > ref:10-15
    18 
    19  11              21        31         41        51        61        71          
    20 GATAA****GATA**GCTGTGCTAGTAGGCAG*TCAGCGCCATNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
    21 .....    ....  ......K.K......K. ..........                                     
    22 GATAAAGAGGATA***CTG      taggc    cagcgccat                                     
    23 GATAAGG**GATAAA                                                                 
    24 GCTAA**** ATA**GCT>>>>>>>>>>>>>>TTCAGC                                          
    25 

That's it,

Pierre

Storing SNPs in a HDF5 file: my notebook

I discovered the benefits of using twitter in 2008, the day Deepak Singh replied to one of my tweets related to the storage of a large number of genotypes.

Since that day, I've tried to use the HDF5 library, without any success (there's a large disheartening documentation/API on the HDF5 site and the API seems to be only used by a small number of geeks). Furthermore, HDF5 is a technology used by the IGV genome browser.

In that post, I'll describe a C program loading a set of SNPs defined by the following C structure:

typedef struct structSnp
 {
 char rsId[RS_LENGTH];//rs##
 char chrom[CHROM_LENGTH];//chromosome name
 int chromStart;//genomic start index
 int chromEnd;//genomic end index
 }Snp;

The SNPs will be read from stdin. Four columns are expected: rs-id, chrom, chromStart and chromEnd. Here is the command line I used to get a small input file:

curl -s "http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/snp132.txt.gz" |\
 gunzip -c |\
 cut -d ' ' -f 2-5 |\
 head -n 10000 > sample.tsv

Source code

The C program is described below. I hope my comments will make the code readable.

The Makefile

Compile and run

gcc -Wall -o a.out -I  /path/to/hdf5/include -L /path/to/hdf5/lib dbsnp2hdf5.c  -lhdf5
./a.out < sample.tsv

Dump the data

At the en, the program creates a structured binary file containing our SNPs. The HDF5 utility h5dump can be used to display the data as text:

$ h5dump storage.h5

HDF5 "storage.h5" {
GROUP "/" {
   GROUP "variations" {
      DATASET "dbSNP" {
         DATATYPE  H5T_COMPOUND {
            H5T_STRING {
               STRSIZE 13;
               STRPAD H5T_STR_NULLTERM;
               CSET H5T_CSET_ASCII;
               CTYPE H5T_C_S1;
            } "rs";
            H5T_STRING {
               STRSIZE 7;
               STRPAD H5T_STR_NULLTERM;
               CSET H5T_CSET_ASCII;
               CTYPE H5T_C_S1;
            } "chrom";
            H5T_STD_I32LE "chromStart";
            H5T_STD_I32LE "chromEnd";
         }
         DATASPACE  SIMPLE { ( 10000 ) / ( H5S_UNLIMITED ) }
         DATA {
         (0): {
               "rs112750067",
               "chr1",
               10326,
               10327
            },
         (1): {
               "rs56289060",
               "chr1",
               10433,
               10433
            },
         (2): {
               "rs112155239",
               "chr1",
               10439,
               10440
            },
         (3): {
               "rs112766696",
               "chr1",
               10439,
               10440
            },
          (...)
               450425,
               450426
            }
         }
         ATTRIBUTE "rdfs:comment" {
            DATATYPE  H5T_STRING {
                  STRSIZE 15;
                  STRPAD H5T_STR_NULLTERM;
                  CSET H5T_CSET_ASCII;
                  CTYPE H5T_C_S1;
               }
            DATASPACE  SCALAR
            DATA {
            (0): "My list of SNPs"
            }
         }
      }
   }
}
}

That's it !
Pierre