Friday, May 25, 2012

Show all bases with IGV


Today someone asked me if I can make a revised version of snpview to support "Show all bases".

Manual operations shows that it required to right-click the alignment and select "Show all bases" from popped menu.

http://www.broadinstitute.org/igv/PortCommands

The supported socket commands for IGV does not include this operation. (In fact it only supports a very few commands)

Since IGV is open resource I decide to check out the source code to see if I can change the code then build a "Show all bases by default" IGV.

1. Find which source file is responsible for this option

>grep -r -n "Show all bases" IGVDistribution_2.1.14/src/

src/org/broad/igv/sam/AlignmentTrack.java:1533:            final JMenuItem item = new JCheckBoxMenuItem("Show all bases");

2. OK let us check this file "src/org/broad/igv/sam/AlignmentTrack.java"

vim +1533 src/org/broad/igv/sam/AlignmentTrack.java


The context shows:
*************************************************************
public void addShowAllBasesMenuItem() {
            // Change track height by attribute
            final JMenuItem item = new JCheckBoxMenuItem("Show all bases");

            if (renderOptions.colorOption == ColorOption.BISULFITE || renderOptions.colorOption == ColorOption.NOMESEQ) {
                //    item.setEnabled(false);
            } else {
                item.setSelected(renderOptions.showAllBases);
            }
            item.addActionListener(new ActionListener() {
                public void actionPerformed(ActionEvent aEvt) {
                    renderOptions.showAllBases = item.isSelected();
                    refresh();
                }
            });
            add(item);
        }
*************************************************************


It looks like "renderOptions.showAllBases" controls the switch on/off.

3. Check the  "renderOptions.showAllBases"
grep -r -n "renderOptions.showAllBases" IGVDistribution_2.1.14/src/

src/org/broad/igv/sam/AlignmentTrack.java:1538:                item.setSelected(renderOptions.showAllBases);
src/org/broad/igv/sam/AlignmentTrack.java:1542:                    renderOptions.showAllBases = item.isSelected();
src/org/broad/igv/sam/AlignmentRenderer.java:556:                if (renderOptions.showMismatches || renderOptions.showAllBases) {
src/org/broad/igv/sam/AlignmentRenderer.java:635:        boolean showAllBases = renderOptions.showAllBases &&



Notice the last line defines a boolean variable "boolean showAllBases ..."

Is this what I am looking for?

4. Further investigation
vim +635 src/org/broad/igv/sam/AlignmentRenderer.java


The context shows:
*************************************************************
// Disable showAllBases in bisulfite mode
        boolean showAllBases = renderOptions.showAllBases &&
                !(colorOption == ColorOption.BISULFITE || colorOption == ColorOption.NOMESEQ);

*************************************************************

Now it is quite clear: "showAllBases" controls the on/off of "Show all bases".

5. Now let us make a dirty patch

*************************************************************
// Disable showAllBases in bisulfite mode
//        boolean showAllBases = renderOptions.showAllBases &&
//                !(colorOption == ColorOption.BISULFITE || colorOption == ColorOption.NOMESEQ);

boolean showAllBases = true;

*************************************************************

We just set "true" to it for always "on"! Like I said, it is dirty because it completely ignore all options there.
It may have potential problem, but at now, let us make a test.

6. Go back to the source and build

>ant


We will got a new "igv.jar". rename it as "igv_showallbases.jar" to discriminate the original one.
>mv igv.jar igv_showallbases.jar

7. Test!

java -jar igv_showallbases.jar

It works! Now by default all bases in alignment are showed. 

The user is happy because he got what he wants!

Tuesday, May 22, 2012

microRNA differential expression

Goal: 
Find the microRNAs that are differentially expressed (DE) between two groups. Each group consists of 6 ~ 7 biological samples which go through illumina HiSeq2000 sequencer. Single-end with 50bp length.

Method:
1. download gene table for microRNA from miRbase
>wget ftp://mirbase.org/pub/mirbase/CURRENT/genomes/hsa.gff3


2. convert it from GFF into BED format(and filter non-microRNA)
>awk '/^#/{print;next;} $3 == "miRNA" {print "chr"$1"\t"$4"\t"$5"\t"$9}' hsa.gff3 > hsa.bed


3. Align reads to hg19 whole genome using novoalign
>novoalign -r None -d hg19.nix -f A_1.fq.gz A_2.fq.gz | gzip > A.sam.gz
>java -jar ~/picard/SortSam.jar INPUT=A.sam.gz OUTPUT=A.bam CREATE_INDEX=true SO=coordinate VERBOSITY=ERROR QUIET=true
>mv A.bai A.bam.bai 


4. Counting with bedtools
>~/BEDTools-Version-2.16.2/bin/bedtools multicov -bams A.bam -bed hsa.bed | awk '$4 !~ 0 {print $0}' > A.txt

(to be continued)