Counting lines in "fastq.gz"

Before splitting the "fastq.gz" into fragments, we need to counting how many lines in the raw file, then we can calculate how many lines per fragment.

If the raw file is quite big, the counting itself will take a long time to process because the file has to be decompressed firstly.

The common approach is using "zcat" with "wc". Is it the fastest?

#1. use zcat - 9.381s 

$time zcat A_R1.fq.gz | wc -l

 20000000

 real    0m9.381s

 user    0m7.423s

 sys     0m1.919s

#2. use zgrep - 16.258s

 $time zgrep -Ec "$" A_R1.fq.gz

 20000000

 real    0m16.258s

 user    0m13.109s

 sys     0m0.490s

#3. use pigz - 7.227s

 $sudo apt-get install pigz

 $time pigz -d -p 1 -c A_R1.fq.gz  | wc -l

 20000000

 real    0m7.227s

 user    0m5.544s

 sys     0m0.386s

#4. use pigz with 4 threads - 5.973s

 $time pigz -d -p 4 -c A_R1.fq.gz  | wc -l

 20000000

 real    0m5.973s

 user    0m5.599s

 sys     0m2.200s

By default, pigz will use all available processors so "-p" is not necessary. The above command can be simplified as 

 $pigz -dc A_R1.fq.gz  | wc -l

Clearly, our winner is pigz.

I am going to write another blog on splitting the "fastq.gz" file with only one goal - as fast as possible. 

 1. Count lines

 2. Determine how many lines per splitted fastq file

 3. Unzip and split the paired "fastq.gz" files

 4. Zip splitted "fastq" files again.

Reference:

http://superuser.com/questions/135329/count-lines-in-a-compressed-file

Build hadoop-2.4.1-src

1. http://apache.cs.utah.edu/hadoop/common/hadoop-2.4.1/hadoop-2.4.1-src.tar.gz

2. http://protobuf.googlecode.com/files/protobuf-2.5.0.tar.bz2

./configure && make && sudo make install && sudo ldconfig

mvn compile && mvn package

VirtualBox 4.3.12 configuration

File->"Preferences"->Click "Network" in left panel -> Add a new NAT network

Name -> "MyNAT"
CIDR: 192.168.0.0/24
Click the "Port forwarding" add port "22" in Ipv4.
Name "SSH"
Protocol "TCP"
Host IP "10.2.2.72" (run "cmd"->"ipconfig" to get it)
Host Port "22"
Guest IP "192.168.0.4" (The IP of the primary virtual machine)
Guest Port "22"


For each virtual machine, in "Settings"->"network"->"Adapter1"->"NAT Network" choose "MyNAT"

Windows+VirtualBox+Ubuntu+Xen+Ubuntu

What we have:

    Native machine installed with Windows 7.

Now I want to play with some virtual machines:

    1. Install VirtualBox into Windows 7, as a VM
    2. Install Ubuntu(1) into VirtualBox, as a guest OS
    3. Install Xen into Ubuntu, as a VM
    4. Install Ubuntu(2) into Xen, as a guest OS

Versions:
    VirtualBox 4.3.10 r93012
    Ubuntu 14.04 Server 
    Xen 4.4


It is quite easy to install Ubuntu into VirtualBox (step1 and step2). One thing to remember: "Do not assign all allocated space to Ubuntu(1)". Supposedly we decide to give 100GB disk space to VMs inside the VirtualBox, then the Ubuntu(1), as a dom0 OS, should use as less space as possible, say, 2GB. The remaining 98GB disk space should be reserved for one or more guest Ubuntu(2) guest OS. Otherwise your Xen will hang up when trying to "brought up CPU".


#1. get the Xen package
$sudo apt-get install xen-hypervisor-4.4-amd64

#2. Set Xen as the default loading OS
$sudo sed -i 's/GRUB_DEFAULT=.*\+/GRUB_DEFAULT="Xen 4.4-amd64"/' /etc/default/grub
sudo update-grub

#3. use xl as toolstack
$sudo sed -i 's/TOOLSTACK=.*\+/TOOLSTACK="xl"/' /etc/default/xen
$sudo reboot

$sudo xl list

#4. Create primary partition on remaining disk space
sudo fdisk /dev/sda

... /dev/sda3

#5. restart 

sudo reboot

#6. create physical storage

sudo pvcreate /dev/sda3

#7. create volume group "xen-vg" on /dev/sda3

sudo vgcreate xen-vg /dev/sda3

#8. take a look at PVs. 

sudo pvs

  PV         VG     Fmt  Attr PSize  PFree

  /dev/sda3  xen-vg lvm2 a--  51.62g 51.62g

#9. create logical volume "ubuntu" guest OS with 10G disk space

sudo lvcreate -L 10G -n ubuntu /dev/xen-vg

#10. Tell it where to get the image - Ubuntu 14.40 (a.k.a utopic)

sudo mkdir -p /var/lib/xen/images/ubuntu-netboot

$cd /var/lib/xen/images/ubuntu-netboot

$sudo wget http://ubuntu.cs.utah.edu/ubuntu/dists/utopic/main/installer-amd64/current/images/netboot/xen/initrd.gz

$sudo wget http://ubuntu.cs.utah.edu/ubuntu/dists/utopic/main/installer-amd64/current/images/netboot/xen/vmlinuz

#11. Configuration

$sudo vim /etc/xen/ubuntu.cfg

name = "ubuntu"

memory = 512

disk = ['phy:/dev/xen-vg/ubuntu,xvda,w']

vif = [' ']

kernel = "/var/lib/xen/images/ubuntu-netboot/vmlinuz"

ramdisk = "/var/lib/xen/images/ubuntu-netboot/initrd.gz"

extra = "debian-installer/exit/always_halt=true -- console=hvc0"

#12. Start installing the Ubuntu(2) 

sudo xl create /etc/xen/ubuntu.cfg -c  

Reference 

https://help.ubuntu.com/community/Xen

LocalCache in hadoop MRv2 aka YARN

The old DistributedCache is deprecated in the new API of hadoop 2.2.3. Now the sugguested method is Job.addCacheFile

By default, the cached file was add to a special folder on each slave node.

Assume we have two slave nodes, namely "n1" and "n2". In the configuration file of "yarn-site.xml" which is under /HADOOP_PATH/etc/hadoop/, you will find a property
like this:

 yarn.nodemanager.local-dirs
 /home/hadoop/localdirs


Here "/home/hadoop/localdirs" is the home path for all cached files. 

>ls /home/hadoop/localdirs
filecache  nmPrivate  usercache


If we add cache files per application, the cache files will be put under "usercache".

Job job = Job.getInstance(conf, "HELLO");
FileStatus[] fileStatus = fs.listStatus(new Path("/user/hadoop/data/"));
for (FileStatus f : fileStatus)
{
job.addCacheFile(new URI(f.getPath().toUri().toString() + "#_" + f.getPath().getName()));
}

Now switch to machine n1 or n2, have a look at the local cache directory you will find it somewhere:

>find /hadoop/localdirs/ -name "hg19.fa"

/hadoop/localdirs/usercache/hadoop/appcache/application_1398026166795_0007/container_1398026166795_0007_01_000010/hg19.fa

To use the cache files in the MR application is interesting. Our Mapper will looks like this

public static class MyMapper extends Mapper
{

@Override
                private ListcacheFiles;
protected void setup(Context context) throws IOException
{
/*
URI[] uris = context.getCacheFiles();
for (URI u : uris)
{
System.out.println("CACHED:" + u.getPath());
}
*/
Path[] uris = context.getLocalCacheFiles();
for (Path u : uris)
{
System.out.println("CACHED:" + u.toString());

}

}

//static Log LOG = LogFactory.getLog(MyMapper.class);
@Override
public void map(Text keySampleName, Text valueSplitFileName, Context context) throws IOException, InterruptedException
{
                        //use the cachedFiles here.
//context.write(keySampleName, valueSplitFileName);
}

}



The new API to retrieve cached files is  "context.getCacheFiles()" while the "context.getLocalCacheFiles()" is deprecated.  However, context.getCacheFiles() returns a URI array which each element is a HDFS path(hdfs://master:50011/user/hadoop/...),
while context.getLocalCacheFiles() returns a Path array which each element is a local path (//hadoop/localdirs/usercache/...)

Besides, each container in every slave node will have a copy of all cached files. Which means in slave node "n1", if "n1" has 10 containers, then you get 10 copies of cached files in n1.

Gene Fusion Detection - deFuse

##########################################
# Input Data, from the 1000G NA12878 RNASeq
##########################################
FASTQ1=ERR356372_1.fastq.gz 
FASTQ2=ERR356372_2.fastq.gz

#######################################
# deFuse
#######################################
#you have to install dependencies include "samtools", "bowtie", "blat", "faToTwoBits", "R" 
#(check out if library "ada" was installed) and "gmap"

#download app
#http://sourceforge.net/apps/mediawiki/defuse/index.php?title=DeFuse_Version_0.6.1

#path
DEFUSE=/home/ysun/app/defuse-0.6.1/scripts

#download lib data
cd /home/ysun/data/
$DEFUSE/create_reference_dataset.pl -c defuse_cbbl.config

#the default create_reference_dataset.pl has one bug on creating the gmap indicies.
#you need to remove "$gmap_index_directory" in "make" command from "sub #create_gmap_indices".
#In my system , the total size of all lib files is 48GB.

#unzip the FASTQ. deFuse will not work on zipped FASTQ?
gunzip *.fastq.gz

#let us start with a small dataset with 100K reads for testing

head -n 400000 ERR356372_1.fastq > A_1.fastq
head -n 400000 ERR356372_2.fastq > A_2.fastq

#now start the testing job
$DEFUSE/defuse.pl -c defuse_cbbl.config -1 A_1.fastq -2 A_2.fastq -o output_defuse

#if everything goes well  you will see outputs like this:
Importing fastq files
Splitting fastq files
Discordant alignments
Read Stats
        Fragment mean 173.540675389827 stddev 59.0075428839599
        Read length min 49 max 49
Generating discordant alignment clusters
Remove mitochondrial-genomic clusters
Generating maximum parsimony solution
Selecting fusion clusters
Preparing sequences for local realignment
Performing local realignment
Filtering concordant clusters
Generating spanning alignment regions file
Initializing split read alignments
Calculating split read alignments
Evaluating split reads
Calculating spanning stats
Calculating spanning p-values
Calculating split read pvalues
Creating fastas
Splitting fastas
Breakpoint alignments
Annotating fusions
Coallating fusions
Running adaboost classifier
Filtering fusions
Success

#now let us go full scale

$time $DEFUSE/defuse.pl -c defuse_cbbl.config -1 ERR356372_1.fastq -2 ERR356372_2.fastq -o output_defuse 

real    1703m53.046s
user    1642m55.346s
sys     83m2.979s

#Kind of slow

The NA12878 RNASeq dataset is very small with 18 million reads and read length is 49. A rough calculation indicates the average coverage is ~60x. Even so, it took more than one day to get it done!

#What's next
Have a glimpse on the result on deFuse. 
Repeat the analysis with Tophat-Fusion and FusionMap.

Learning CoffeeScript

Supposedly I want to declare a closure nested method in JavaScript like this:

###Sample1
//declaration
 function addInternal(x)
 { 
    //return anonymous nested function
     return function(y) { console.log(x+y);};
 }

//use
var addByTen  = addInternal(10);
addByTen(15);

.

or wrap it up as a variable

###Sample2, assign the parameter to outer function addInternal in declaration
var addByTen  = (function addInternal(x)
{
    //return anonymous nested function
    return function(y){console.log(x+y);};
})(10);

//use
addByTen(15);


The corresponding CoffeeScript version is very tight:

###sample1
addInternal = (x) ->
 (y) -> console.log x+y

addByTen = addInternal 10
addByTen 15

###sample2
addByTen = ((x) ->
 (y) -> console.log x+y ) (10)

addByTen 15