In [125]:
#Toggle code
from IPython.display import HTML
HTML('''<script>
code_show=true;
function code_toggle() {
if (code_show){
$('div.input').hide();
} else {
$('div.input').show();
}
code_show = !code_show
}
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')
Out[125]:
In [1]:
#load modules
import os
import subprocess
import fnmatch
from IPython.core.display import Image
In [2]:
#go to working directory, where the data is
os.chdir("/nagyvinyok/adat83/sotejedlik/ribli/dt40/snp/varscan")
In [18]:
#create sql util function
import subprocess
def run_sqlilte3(command,db,remove_db=False,output=None):
tempf=open('tempf.sql','w')
tempf.write(command)
tempf.close()
if (remove_db):
subprocess.check_output('rm '+db, shell=True)
if (output == None):
try:
print subprocess.check_output('sqlite3 '+ db + ' < tempf.sql ', shell=True,
stderr=subprocess.STDOUT)
except subprocess.CalledProcessError, e:
print e.output
else:
try:
print subprocess.check_output('sqlite3 '+ db + ' < tempf.sql > '+ output, shell=True,
stderr=subprocess.STDOUT)
except subprocess.CalledProcessError, e:
print e.output
subprocess.check_output(['rm','tempf.sql'])
In [88]:
#load the vcf files into sqlite
#headers??
run_sqlilte3('''
.separator "\t"
CREATE TABLE indels(
chrom TEXT,
position INTEGER, --!!!
ref TEXT,
var TEXT,
normal_reads1 TEXT,
normal_reads2 TEXT,
normal_var_freq TEXT,
normal_gt TEXT,
tumor_reads1 TEXT,
tumor_reads2 TEXT,
tumor_var_freq TEXT,
tumor_gt TEXT,
somatic_status TEXT,
variant_p_value TEXT,
somatic_p_value TEXT,
tumor_reads1_plus TEXT,
tumor_reads1_minus TEXT,
tumor_reads2_plus TEXT,
tumor_reads2_minus TEXT,
normal_reads1_plus TEXT,
normal_reads1_minus TEXT,
normal_reads2_plus TEXT,
normal_reads2_minus TEXT
);
CREATE TABLE snps(
chrom TEXT,
position INTEGER, --!!!
ref TEXT,
var TEXT,
normal_reads1 TEXT,
normal_reads2 TEXT,
normal_var_freq TEXT,
normal_gt TEXT,
tumor_reads1 TEXT,
tumor_reads2 TEXT,
tumor_var_freq TEXT,
tumor_gt TEXT,
somatic_status TEXT,
variant_p_value TEXT,
somatic_p_value NUMERIC, --!!!
tumor_reads1_plus TEXT,
tumor_reads1_minus TEXT,
tumor_reads2_plus TEXT,
tumor_reads2_minus TEXT,
normal_reads1_plus TEXT,
normal_reads1_minus TEXT,
normal_reads2_plus TEXT,
normal_reads2_minus TEXT
);
.import DS14_DS15.vsc.snp snps
.import DS14_DS15.vsc.indel indels
--delete headers
DELETE FROM indels WHERE chrom='chrom';
DELETE FROM snps WHERE chrom='chrom';
--create index to make next join faster?
CREATE INDEX indels_chrom_pos_index ON indels (chrom,position);
CREATE INDEX snps_chrom_pos_index ON snps (chrom,position);
''',db='variantdb',remove_db=True)
The number of somatic mutations¶
- using p_value < 0.01
- normal is clean
In [129]:
run_sqlilte3('''
.separator "\t"
SELECT COUNT(*)
FROM snps AS s
WHERE
s.somatic_status='Somatic' AND
s.somatic_p_value < 0.01 AND
s.normal_reads2 = '0';
''',db='variantdb',remove_db=False)
The number of somatic mutations after indel proximity filtering¶
- using p_value < 0.01
- normal is clean
- indel proximity filtering
In [128]:
run_sqlilte3('''
.separator "\t"
--EXPLAIN QUERY PLAN SELECT COUNT(*) --i.chrom,i.position,s.chrom,s.position
SELECT COUNT(*) --i.chrom,i.position,s.chrom,s.position
FROM snps AS s
LEFT JOIN indels AS i ON (
i.chrom = s.chrom AND i.position BETWEEN s.position-150 AND s.position+150 )
WHERE
i.chrom IS NULL AND
s.somatic_status='Somatic' AND
s.somatic_p_value < 0.01 AND
s.normal_reads2 = '0';
''',db='variantdb',remove_db=False)
List some of them in the varscan output format¶
- headers can be peeked from the code
In [130]:
run_sqlilte3('''
.separator "\t"
--EXPLAIN QUERY PLAN SELECT COUNT(*) --i.chrom,i.position,s.chrom,s.position
--SELECT COUNT(*) --i.chrom,i.position,s.chrom,s.position
SELECT s.*
FROM snps AS s
LEFT JOIN indels AS i ON (
i.chrom = s.chrom AND i.position BETWEEN s.position-150 AND s.position+150 )
WHERE
i.chrom IS NULL AND
s.somatic_status='Somatic' AND
s.somatic_p_value < 0.01 AND
s.normal_reads2 = '0'
LIMIT 20;
''',db='variantdb',remove_db=False)
In [116]:
#collect filenames for samplenames
bam_link_dir='/nagyvinyok/adat83/sotejedlik/ribli/dt40/bam_links/'
#collect filenames
fnames=[]
for fname in os.listdir(bam_link_dir):
if (fnmatch.fnmatch(fname, '*.bam') and
not fnmatch.fnmatch(fname,"*.bai")):
fnames.append(fname)
fnames=sorted(fnames)
In [117]:
#define function to run pileup command
def run_pileup(chrom,pos):
samtools="/nagyvinyok/adat87/home/ribli/tools/samtools-0.1.19/samtools"
ref_fa="/nagyvinyok/adat87/home/ribli/input/index/gallus/Gallus_gallus.Galgal4.74.dna.toplevel.fa"
cmd_mpileup = samtools + ' mpileup -Q 30 -B -f ' + ref_fa
cmd_mpileup += ' -r ' + str(chrom) +':'+ str(pos) +'-' + str(pos) + ' '
cmd_mpileup += ' '.join([bam_link_dir+x for x in fnames])
try:
pup_line=subprocess.check_output(cmd_mpileup,executable='/bin/bash',
stderr=subprocess.STDOUT,shell=True)
except subprocess.CalledProcessError,e:
print e
pup_list=pup_line.split('\t')
print pup_list[0],pup_list[1],pup_list[2]
for i in xrange((len(pup_list)-3)/3):
print fnames[i].split('.')[0][:-21],pup_list[3+3*i],pup_list[4+3*i]
Check a pileup¶
In [137]:
chrom,pos=1,634212
run_pileup(chrom,pos)
Check it in IGV¶
In [153]:
#connect to IGV socket
import socket
SERVER_IP = 'localhost'
PORT_NUMBER = 60151
SIZE = 1024
my_socket = socket.socket( socket.AF_INET, socket.SOCK_STREAM )
my_socket.connect((SERVER_IP, PORT_NUMBER))
#Load the bams, and genome
def load_bams(my_socket,bams,genome):
my_socket.sendall('new \n') #clear it
print my_socket.recv(1024),
my_socket.sendall('genome ' + genome + ' \n') #load genome
print my_socket.recv(1024),
for bam in bams:
my_socket.sendall('load ' + bam + '\n' ) # load the bams
print my_socket.recv(1024),
return
def make_snapshot(my_socket, snapshor_dir,chromosome,posmin,posmax):
my_socket.sendall('goto '+chromosome+':'+str(posmin)+'-'+str(posmax)+'\n')
print my_socket.recv(1024),
my_socket.sendall('snapshotDirectory ' + snapshor_dir + '\n')
print my_socket.recv(1024),
my_socket.sendall('snapshot \n')
print my_socket.recv(1024),
return
bam1='/nagyvinyok/adat83/sotejedlik/ribli/dt40/bam_links/DS014_RMdup_picard_realign.bam'
bam2='/nagyvinyok/adat83/sotejedlik/ribli/dt40/bam_links/DS015_RMdup_picard_realign.bam'
genome='/home/ribli/input/index/gallus/Gallus_gallus.Galgal4.74.dna.toplevel.fa'
load_bams(my_socket,[bam1,bam2],genome)
snapshot_dir='/home/ribli/DT40/SNP/varscan/indel_prox/'
In [154]:
chromosome='1'
pos=634212
posmin,posmax=pos-100,pos+100
make_snapshot(my_socket,snapshot_dir,chromosome,posmin,posmax)
In [155]:
Image('/home/ribli/DT40/SNP/varscan/indel_prox/1_634,112_634,312.png')
Out[155]:
Check another pos¶
In [121]:
chrom,pos=1,670994
run_pileup(chrom,pos)
In [164]:
chromosome='1'
pos=670994
posmin,posmax=pos-100,pos+100
make_snapshot(my_socket,snapshor_dir,chromosome,posmin,posmax)
Weird region¶
In [166]:
Image('/home/ribli/DT40/SNP/varscan/indel_prox/1_670,894_671,094.png')
Out[166]:
Check another one¶
In [132]:
chrom,pos=1,17064922
run_pileup(chrom,pos)
In [161]:
chromosome='1'
pos=17064922
posmin,posmax=pos-100,pos+100
make_snapshot(my_socket,snapshor_dir,chromosome,posmin,posmax)
In [162]:
Image('/home/ribli/DT40/SNP/varscan/indel_prox/1_17,064,822_17,065,022.png')
Out[162]:
