I am running out of disk space

again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again again

Multiple hits? All or None?

#-r All
57G Oct 14 09:44 1.sam
65G Oct 14 10:06 2.sam
65G Oct 14 10:26 3.sam
45G Oct 14 10:46 4.sam
56G Oct 14 11:09 5.sam


#-r None
56G Oct 11 02:10 1.sam
64G Oct 11 02:32 2.sam
63G Oct 11 02:53 3.sam
44G Oct 11 03:07 4.sam
54G Oct 11 03:25 5.sam

pre-processing of masked FASTA for novoalign

1. Convert lower case to upper case
sed -i 's/\(.*\)/\U\1/' hg19.fasta

> chr1
NNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNN
ggggggttttTGaaaaaaaCCC

Will be


> CHR1
NNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNN
GGGGGGGTTTTTTGAAAAACCC

2. covert '>CHR' to '>chr'
sed -i -e 's/>CHR/>chr/g' hg19.fasta

3. convert N to n
sed -i -e 's/N/n/g' hg19.masked.fasta

4. mask these UTR
novoindex -m hg19.nix hg19.fasta

Keyless ssh connection between A and B

A:
ssh-keygen -t rsa
...
...

Then scp the "id_rsa.pub" to B


B:
cat id_rsa.pub >> .ssh/authorized_keys
chmod 644 .ssh/authorized_keys


Now connection from A to B is keyless.

ACL Setup

On the job side, I need a root-like user who has write permission on any sub-folders under PATH_JOB while the individual user has only 744 like permissions on his/her own job folders. However, Unix/Linux does not have fine control over file/directory by default, so you can not assign permissions to different users. ACL (Access Control Lists) can be used to solve this problem.

Here is a simple test on my VM_Ubuntu_10.10:

1. List partitions, the /dev/sda1 is the root directory "/"
# blkid
/dev/sda1: UUID="xxxx" TYPE="ext4"
/dev/sda5: UUID="yyyy" TYPE="swap"

2. Modify the /etc/fstab to support ACL by inserting "acl" into options.
#vim /etc/fstab
UUID=xxxx / ext4 errors=remount-ro,acl 0 1

3. Remount the root directory
#mount / -o remount

4. Make a test
#setfacl -m u:abc:rw hello
#getfacl hello

# file: hello.txt
# owner: root
# group: root
user::rw-
user:abc:rw-
group::rw-
mask::rw-
other::---

5. Now user "abc" can read and write "hello"!

#ls
rw-rw----+ 1 root root 0 2011-09-16 11:41 hello.txt

Note the new "+" character appears the end of permission field. This shows this file/folder supports ACL

Merge the SAM/BAM file

Suppose I have two alignments in SAM (1.sam and 2.sam). Since these two alignments are two lanes from the same sample, I should merge these two alignments into one.


>wc -l 1.sam
10000

>wc -l 2.sam
10000

So both the SAM files has 10000 lines (include the header)

First using picard's "MergeSamFiles.jar"

>java -jar MergeSamFiles.jar I=1.sam I=2.sam O=12.sam
>wc -l 12.sam
19973

The missing 27 lines are actually duplicated header lines and they were removed by PICARD automatically. However, there is a subtle bug here.

>less 12.sam

@HD VN:1.0 GO:none SO:coordinate
@SQ SN:chr1 LN:249250621 AS:hg19.nov.illumina.nix
@SQ SN:chr2 LN:243199373 AS:hg19.nov.illumina.nix
@SQ SN:chr3 LN:198022430 AS:hg19.nov.illumina.nix
@SQ SN:chr4 LN:191154276 AS:hg19.nov.illumina.nix
@SQ SN:chr5 LN:180915260 AS:hg19.nov.illumina.nix
@SQ SN:chr6 LN:171115067 AS:hg19.nov.illumina.nix
@SQ SN:chr7 LN:159138663 AS:hg19.nov.illumina.nix
@SQ SN:chr8 LN:146364022 AS:hg19.nov.illumina.nix
@SQ SN:chr9 LN:141213431 AS:hg19.nov.illumina.nix
@SQ SN:chr10 LN:135534747 AS:hg19.nov.illumina.nix
@SQ SN:chr11 LN:135006516 AS:hg19.nov.illumina.nix
@SQ SN:chr12 LN:133851895 AS:hg19.nov.illumina.nix
@SQ SN:chr13 LN:115169878 AS:hg19.nov.illumina.nix
@SQ SN:chr14 LN:107349540 AS:hg19.nov.illumina.nix
@SQ SN:chr15 LN:102531392 AS:hg19.nov.illumina.nix
@SQ SN:chr16 LN:90354753 AS:hg19.nov.illumina.nix
@SQ SN:chr17 LN:81195210 AS:hg19.nov.illumina.nix
@SQ SN:chr18 LN:78077248 AS:hg19.nov.illumina.nix
@SQ SN:chr19 LN:59128983 AS:hg19.nov.illumina.nix
@SQ SN:chr20 LN:63025520 AS:hg19.nov.illumina.nix
@SQ SN:chr21 LN:48129895 AS:hg19.nov.illumina.nix
@SQ SN:chr22 LN:51304566 AS:hg19.nov.illumina.nix
@SQ SN:chrX LN:155270560 AS:hg19.nov.illumina.nix
@SQ SN:chrY LN:59373566 AS:hg19.nov.illumina.nix
@SQ SN:chrM LN:16571 AS:hg19.nov.illumina.nix
@RG ID:TOMATO PL:ILLUMINA LB:LIBTMP SM:SAMPLE
@PG ID:novoalign VN:V2.07.10 CL:novoalign
@PG ID:novoalign.1 VN:V2.07.10 CL:novoalign



The last two lines are both @PG which differ the ID. Before merging, they are both "novoalign", after merging, the duplicated ID was renamed as "novoalign.1" (if merging 3 files, you will have "novoalign.2" as well). When you check the individual read alignment, you will see the PG:Z field referenced the PG line. When you call variances later, the genotyper application will process merged file as two files which is not what you want.


To fix this:

1. Keep the first @PG line and delete all other @PG lines
2. Replace all "novoalign.X" to "novoalign".

Or you can use samtools.

>java -jar SortSam.jar INPUT=1.sam OUTPUT=1.bam SO=unsorted
>java -jar SortSam.jar INPUT=2.sam OUTPUT=2.bam SO=unsorted
>samtools merge 3.bam 1.bam 2.bam
>samtools view -h -o 3.sam 3.bam
>less 3.sam

Samtools only use the first SAM's header in merged SAM, so it is OK in this case.

mouse dbSNP VCF

OK I spend few hours on searching data and writing a script to make a dbSNP VCF for mouse. The steps are simple and dirty. Test of this VCF file has not been done yet. Put it here for reference only.

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

#1. download dbSNP mouse
wget --continue --no-host-directories --no-directories --preserve-permissions ftp://ftp.ncbi.nih.gov/snp/organisms/mouse_10090/chr_rpts/*

#2. unzip
gunzip *.gz

#3. exclude all SNPs that are not mapped to mm9 (MGSCv37) genome.
for i in `ls`; do grep 'MGSCv37' "$i" > mm9."$i" ; done


#However, you can NOT find REF/ALT field in above files. So we must get it from other files.
#4. download dbSNP mouse with genotype
wget --continue --no-host-directories --no-directories --preserve-permissions ftp://ftp.ncbi.nih.gov/snp/organisms/mouse_10090/genotype/*

#5. unzip
gunzip *.gz

#6. check out these files, should end with ".xml"
ls

Black6GtyFromContig gt_chr13.xml gt_chr17.xml gt_chr2.xml gt_chr6.xml gt_chrAltOnly.xml gt_chrUn.xml
gt_chr10.xml gt_chr14.xml gt_chr18.xml gt_chr3.xml gt_chr7.xml gt_chrMT.xml gt_chrX.xml
gt_chr11.xml gt_chr15.xml gt_chr19.xml gt_chr4.xml gt_chr8.xml gt_chrMulti.xml gt_chrY.xml
gt_chr12.xml gt_chr16.xml gt_chr1.xml gt_chr5.xml gt_chr9.xml gt_chrNotOn.xml README

#7. make a 3 column file: "id REF ALT"
grep '$//g' | sed 's/rsId=//g' | sed 's/observed=//g' | awk -F '/' '{print $1, $2}' > id.allele.txt



#8. write a simple python script that mapping id from mm9.chr_N.txt to make the final VCF file - "mm9.dbsnp.vcf"

def mm9():
    m = {}
    for line in file('id.allele.txt'):
        line = line.strip()
        toks=line.split(' ')
        if len(toks)==3:
            m[toks[0]] = (toks[1],toks[2])
    
    #chr1 - chr19, chr X,Y and MT only.
    nfs = [str(i) for i in range(1,20)]
    nfs.append('X')
    nfs.append('Y')
    nfs.append('MT')
    
    fo = file('mm9.dbsnp.vcf','w')
    print >>fo,'##fileformat=VCFv4.0'
    print >>fo,'##fileDate=20110712'
    print >>fo,'##source=dbSNP'
    print >>fo,'##dbSNP_BUILD_ID=132'
    print >>fo,'##reference=mm9'
    print >>fo,'#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO'
    
    total = 0
    for n in nfs:
        sm = {}
        fn = os.path.join('mm9.chr_%s.txt'%n)
        for line in file(fn):
            line = line.strip()
            toks=line.split('\t')
            try:
                rs_id = toks[0]
                mapped = toks[1]
                withdraw =toks[2]  
                chr =toks[6]  
                position =toks[11]
                if m.has_key(rs_id) and mapped=='2' and withdraw=='0':
                    ref,alt = m[rs_id]
                    #use an arbitratry "ms"+ID as prefix refID (instead of "rs")
                    sm[int(position)] = [chr,position,'ms'+rs_id,ref,alt,'.','PASS','.'] 
            except Exception,e:
                pass
        sk = sorted(sm.keys())
        for k in sk:
            total += 1
            print >>fo,'\t'.join(sm[k])
    fo.close()
    print total
mm9()    

#9. make sure no duplicated IDs
awk '{print $3}' mm9.dbsnp.vcf | sort| uniq -d



13745646
13501994