Friday, October 2, 2015

Add and delete users in batch

#add some users to group "research" with initial password "123"

g=research
for u in t1 t2 t3
do 
sudo useradd -d /home/${u} -m ${u} -g ${g}
echo ${u}:${u}123 | sudo chpasswd
echo -ne "${u}123\n${u}123\n" | sudo smbpasswd -a  ${u}
done


#delete these users as well as their home directory

for u in t1 t2 t3
do 
sudo userdel -r $u 2>/dev/null
done


#configure samba

sudo apt-get install -y samba
cd /etc/samba/
sudo cp smb.conf smb.conf.bak
sudo vim smb.conf

############################
[global]
security = user

[store]
path = /store/public
valid users = @research
browsable =yes
writable = yes
guest ok = no
############################

sudo mkdir -p /store/public
sudo chmod -R 0755 /store/public
sudo chown -R hadoop:research /store/public
sudo service smbd restart

Thursday, September 10, 2015

Redirect nohup's output

Some novice Linux users like to execute a long time running pipeline using nohup 

"nohup mycommand"

which by default will direct stdout and stderr to the current working directory with name "nohup.out". If we have many of these files then we can easily get lost due to the same names under different folders.

While we can tell the user to redirect their output to another managed path, it is hard to let them run something like 

"nohup 1>/log/a.log 2>&1" 

It is just too complicated.

TO solve this problem without changing their old habit, we can re-define the nohup. Add following code to the end of "/etc/bash.bashrc"

function nohup () {
UUID=$(cat /proc/sys/kernel/random/uuid)
logfile=/log/tmp/${UUID}.txt
printf  "log:\n${logfile}\n"
/usr/bin/nohup $@ 1> ${logfile} 2>&1}

Next time when user launched a nohup, the system will print a message like this:

log:

/log/tmp/16170f02-7413-49fb-a8ea-0b5c621cdd93.txt


Tuesday, August 18, 2015

Evaluation on speedseq


Recently read an interesting article from here:

http://www.nature.com/nmeth/journal/vaop/ncurrent/full/nmeth.3505.html

The author is kind enough to open source it
https://github.com/hall-lab/speedseq


So I checked out the code and make some tests. Why is it faster than popular GATK package? After spending few hours , my conclusions:


1. The speedseq use simplified process.

Basically it use works like:

"mapping"->"remove duplicates"->"Call variants"

while traditional GATK works like:

"mapping"->"remove duplicates"->"realign"->"recalibrate"->"Call variants" plus lots of sorting and indexing between steps.


As indicated in freebayes' website, freebayes internally handle realign and recalibrate.



2. User GNU parallel to parallel the variants calling, because freebayes does not support threading or sub-processing.


So supposedly you assign 8 CPU cores to the speedseq, it will launch a 8-processes pool to call variants using freebayes on each chromosome independently.

3. Other improvements.

For example, use sambamba instead of picard (in my experience sambamba is at least 3x faster than picard). use samblaster to get discordant and split reads simultaneously for later SV and CNV callings right after the alignment.


So what are the catches? Nothing can be perfect.

1. The components are tightly wrapped into the pipeline. Changing of individual parameters are no easy tasks.

2. The variant calling process did not split the BAM actually. If fact it only splits the BAM header for a chromosome based region to accelerate the processing because the caller will not need to go through the whole genome. Each process has to load the full size reference genome and full size BAM file(s). As a result, the memory requirement will be a problem for some machines. Also the disk I/O could be a bottleneck if you have many processes trying to read/write large data on the same time.

No wonder why the authors did the test on a pretty good AWS EC2 c3.8xlarge instance with big memory and large SSD.



Final words: Even it is not perfect, it is a wonderful solution for some people. The best of all, it is free! No need to pay greedy and lovely Broad Institute $$$$$$$$$$$$$$$$$$ for the full GATK package.


Increasing the amount of inotify watchers

echo fs.inotify.max_user_watches=524288 | sudo tee -a /etc/sysctl.conf && sudo sysctl -p

Saturday, August 15, 2015

The famous NA12878 WGS at 50x

Due to the lacking of true variants from real genome sequencing, many tools use NA12878 for benchmarking. For whole genome sequencing at 50x depth, we can download these two files: ftp.sra.ebi.ac.uk/vol1/fastq/ERR194/ERR194147/ERR194147_1.fastq.gz ftp.sra.ebi.ac.uk/vol1/fastq/ERR194/ERR194147/ERR194147_2.fastq.gz
Some statistics of these two files

NA12878 WGS 50x
Facts ERR194147_1.fastq.gz ERR194147_2.fastq.gz
URI ERR194147_1.fastq.gz ERR194147_2.fastq.gz
Size 48G 49G
Reads# 787265109 787265109
Reads Length 101 101
Coverage 101*787265109/3000000000=26.5 101*787265109/3000000000=26.5

Friday, June 12, 2015

Download sequencing data from illumina's basespace

Spent 30 mins creating a script to download sequencing data from illumina's BaseSpace, by using its REST API

Here is the code:

https://gist.github.com/anonymous/91c5accc3988cb233347

Tuesday, May 12, 2015

Simple SAMBA

sudo apt-get install -y samba
mkdir /project/share
chmod -R a+rxw /project/share
###############################
#/etc/samba/smb.conf
[share]
path = /project/share
available = yes
guest only=yes
read only = no
browseable = yes
public = yes
writable = yes
###############################
sudo service smbd restart

\\10.2.5.212\share