Intersecting VCF

We can use vcftools (together with bgzip and tabix) to intersect VCF files. The first 4 columns for each VCF file(chr, position, ref and alt) can be easily overlapped. However it is undefined on how to "intersecting scores" in VCF files because it is not a simple arithmetic operation. To generate a VCF format compatible file, vcftools used the scores from the first VCF file only when intersecting multiple VCF files. This is simple but not a perfect solution.

Another possible solution is to merge the BAM firstly before calling variances. Therefore we will have only one VCF file after calling variances.

Is this perfect solution? Maybe. Merging alignment BAM files from different samples (individuals) is not well studied. We need a well known dataset to benchmark these two methods.

MarkDuplicate is very necessary

Recently I did an analysis which basically comparing the positions of nucleosome in two samples to see if the pattern differs.

The target region is only 1500bp in length. HiSeq generated big raw read sequence files. After aligning to hg19, I can see many alignemnts are incorrect because they were mapped to other chromosomes and positions.

The first analysis without MarkDuplicate from Picard in BAM files result in huge peaks of depth of coverage which basically drowns any useful signals. After removing duplicates, the peaks were totally removed and I got strong signals.

Later I found another problem in GATK's DepthOfCoverage which reports 1000 as maximum depth of coverage. For example, even the actual depth of coverage on two positions are 5000 and 1020 respectively, GATK reports the same 1000.

As an alternative solution, samtools can report the actual depth of coverage without truncating.

samtools mpileup -ABr chr1:12345-123456 -d 1000000 x1.bam | awk '{print $4}' > x1.cov
samtools mpileup -ABr chr1:12345-123456 -d 1000000 x2.bam | awk '{print $4}' > x2.cov

Extracting the 4th column that has the depth of coverage.

The final step is plotting with R.

x1=as.numeric(unlist(read.table('x1.cov')))
x2=as.numeric(unlist(read.table('x2.cov')))
pdf(file="1.pdf") 
op = par( mfrow = c( 3, 1 ))
barplot(x1)
title("X1")
barplot(x2)
title("X2")
op = par(op) 
dev.off() 
pdf(file="2.pdf") 
plot(x1,type='l',col='green',xlab='chr1:12345-123456',ylab='Depth of Coverage')
lines(x2,type='l',col='red')
legend(10, 8000, c("X1", "X2"), col = c("green","red"), text.col = "green4", lty = c(3,4), pch = c(3,4), merge = TRUE, bg = 'gray90')
dev.off()

Remove any variances that were in CEU population

#use CEU as common variance database, not dbSNP
import sys
def filter_ceu(target_file,background_file='CEU.all.ids'):
    x = []
    for line in file(background_file):
        x.append(line.strip())
    mx = {}.fromkeys(x)
    for line in file(target_file):
        line = line.strip()
        if line.startswith('#'):
            pass
            #print line
        else:
            toks = line.split('\t')
            if not mx.has_key(toks[2]):
                print line
            else:
                pass
    
if __name__=='__main__':
    filter_ceu(sys.argv[1])

1000Genome Population

CHB Han Chinese Han Chinese in Beijing, China
CHS Southern Han Chinese Han Chinese South
CDX Dai Chinese Chinese Dai in Xishuangbanna, China
CHD Denver Chinese Chinese in Denver, Colorado (pilot 3 only)
JPT Japanese Japanese in Tokyo, Japan
KHV Kinh Vietnamese Kinh in Ho Chi Minh City, Vietnam

CEU CEPH Utah residents (CEPH) with Northern and Western European ancestry
TSI Tuscan Toscani in Italia
GBR British British in England and Scotland
FIN Finnish Finnish in Finland
IBS Spanish Iberian populations in Spain

YRI Yoruba Yoruba in Ibadan, Nigeria
LWK Luhya Luhya in Webuye, Kenya
GWD Gambian Gambian in Western Division, The Gambia
GHN Ghanaian Ghanaian in Navrongo, Ghana
MAB Malawian Malawian in Blantyre, Malawi

ASW African-American SW African Ancestry in Southwest US
AJM African-American MS African American in Jackson, Mississippi
ACB African-Caribbean African Caribbean in Barbados

MXL Mexican-American Mexican Ancestry in Los Angeles, California
CLM Colombian Colombian in Medellin, Colombia
PEL Peruvian Peruvian in Lima, Peru
PUR Puerto Rican Puerto Rican in Puerto Rico

Hemizygous

Describes an individual who has only one member of a chromosome pair or chromosome segment rather than the usual two;
refers in particular to X-linked genes in males who under usual circumstances have only one X chromosome.

Minus and plus strand

ACTGGCGGCCATGGAACTCACCG

X - 15353623 15353645

SO this mapped to the minus strand

To match the plus strand, first reverse the sequence to:

GCCACTCAAGGTACCGGCGGTCA

Then compliment the reversed sequence to:

CGGTGAGTTCCATGGCCGCCAGT

Then it will be mapped to:

X + 15353623 15353645

Trac and Mercurial

#mkdir tool; cd tool

#Python-2.7.2
wget http://www.python.org/ftp/python/2.7.2/Python-2.7.2.tgz
tar -zxvf Python-2.7.2.tgz
cd Python-2.7.2
./configure --prefix=/home/u0592675/tool/python-2.7.2
make
make install


#apache
wget http://www.eng.lsu.edu/mirrors/apache//httpd/httpd-2.2.21.tar.gz
tar -zxvf httpd-2.2.21.tar.gz
cd httpd-2.2.21
./configure --prefix=/home/u0592675/tool/apache-2.2.21
make
make install

#mod_wsgi
wget http://modwsgi.googlecode.com/files/mod_wsgi-3.3.tar.gz
tar -zxvf mod_wsgi-3.3.tar.gz
cd mod_wsgi-3.3
./configure --with-apxs=/home/u0592675/tool/apache-2.2.21/bin/apxs --with-python=/home/u0592675/tool/python_2.7.2/bin/python
make
make install


#Trac
wget http://ftp.edgewall.com/pub/trac/Trac-0.12.2.tar.gz
wget http://peak.telecommunity.com/dist/ez_setup.py
~/tool/python-2.7.2/bin/python ez_setup.py



#
trac-admin /home/ying/mercurial_database/tomato initenv
trac-admin /home/ying/mercurial_database/tomato deploy /home/ying/web/tomato_trac


{'mod_wsgi.listener_port': '80',
'HTTP_COOKIE': 'trac_form_token=0b8daa13a5ec1c5f4819083c; trac_session=80547f6793f2d13cbfea5b1c; trac_form_token=a30953232455e34b94f82aa5; trac_session=d5bdc2b3269ef9a4001807dd',
'mod_wsgi.listener_host': '',
'SERVER_SOFTWARE': 'Apache/2.2.21 (Unix) mod_wsgi/3.3 Python/2.7.2+',
'SCRIPT_NAME': '/tomato',
'mod_wsgi.handler_script': '',
'SERVER_SIGNATURE': '',
'REQUEST_METHOD': 'GET',
'PATH_INFO': '/',
'SERVER_PROTOCOL': 'HTTP/1.1',
'QUERY_STRING': '',
'HTTP_ACCEPT_CHARSET': 'ISO-8859-1,utf-8;q=0.7,*;q=0.3',
'HTTP_USER_AGENT': 'Mozilla/5.0 (Windows NT 6.1; WOW64) AppleWebKit/535.2 (KHTML, like Gecko) Chrome/15.0.874.106 Safari/535.2', 'HTTP_CONNECTION': 'keep-alive', 'SERVER_NAME': '155.100.213.215', 'REMOTE_ADDR': '155.100.213.215', 'mod_wsgi.request_handler': 'wsgi-script', 'wsgi.url_scheme': 'http', 'PATH_TRANSLATED': '/home/ying/tool/apache-2.2.21/htdocs/index.html', 'SERVER_PORT': '80', 'wsgi.multiprocess': True, 'mod_wsgi.input_chunked': '0', 'SERVER_ADDR': '192.168.247.129', 'DOCUMENT_ROOT': '/home/ying/tool/apache-2.2.21/htdocs', 'mod_wsgi.process_group': '', 'SCRIPT_FILENAME': '/home/ying/trac/cgi-bin/trac.wsgi', 'SERVER_ADMIN': 'you@example.com', 'wsgi.input': , 'HTTP_HOST': '155.100.213.215', 'wsgi.multithread': False, 'mod_wsgi.callable_object': 'application', 'REQUEST_URI': '/tomato/', 'HTTP_ACCEPT': 'text/html,application/xhtml+xml,application/xml;q=0.9,*/*;q=0.8', 'wsgi.version': (1, 1), 'GATEWAY_INTERFACE': 'CGI/1.1', 'wsgi.run_once': False, 'wsgi.errors': , 'REMOTE_PORT': '6064', 'HTTP_ACCEPT_LANGUAGE': 'en-US,en;q=0.8', 'mod_wsgi.version': (3, 3), 'mod_wsgi.application_group': '',
'mod_wsgi.script_reloading': '1',
'wsgi.file_wrapper': ,
'HTTP_ACCEPT_ENCODING': 'gzip,deflate,sdch'}

[Tue Nov 08 11:15:24 2011] [error]