In this notebook i will try to explore VarScan 'SNP' positions in a R2 untreated sample: DS14¶
- Im trying to look at alignment context information using IGV, and all sample information using mpileup
In [24]:
#load modules
import os
import subprocess
import fnmatch
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
import gc
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 [3]:
%%bash
#Check out the schema of the vsc file
head DS14_DS15.vsc.snp -n 1 | awk '{for (i=1; i<=NF;i++) {print i" " $i}}'
Look at positions with very low p-value¶
p < 0.0001
It looks like there is some cluestering: e.g.: 4.9M, 10.9M
- LOH in germline pair?
In [10]:
#select some very unlikely fals positives ( p< 0.0001)
best_pos_cmd='''
head DS14_DS15.vsc.snp -n 600000 | awk ' $13 == "Somatic" && $15 < 0.0001 { print }'
'''
best_pos=subprocess.check_output(best_pos_cmd,executable='/bin/bash',shell=True).strip()
best_pos_list=best_pos.split('\n')
print '\n'.join(best_pos_list)
In [11]:
#collect filenames for samplenames
rm_dup_dir='/nagyvinyok/adat83/sotejedlik/dt40/rmdup/'
#collect filenames
fnames=[]
for fname in os.listdir(rm_dup_dir):
if (fnmatch.fnmatch(fname, '*.bam') and
not fnmatch.fnmatch(fname,"*.bam.bai")):
fnames.append(fname)
fnames=sorted(fnames)
In [12]:
#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([rm_dup_dir+x for x in fnames])
pup_line=subprocess.check_output(cmd_mpileup,executable='/bin/bash',shell=True)
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]
In [28]:
Image('/home/ribli/DT40/SNP/varscan/igv_snapshots/ds14__hits/1-3089329.png')
Out[28]:
Mpileup line:
In [13]:
i=0
print best_pos_list[i]
print
run_pileup(1,best_pos_list[i].split('\t')[1])
In [30]:
Image('/home/ribli/DT40/SNP/varscan/igv_snapshots/ds14__hits/1-7971336.png')
Out[30]:
In [14]:
i=1
print best_pos_list[i]
print
run_pileup(1,best_pos_list[i].split('\t')[1])
In [31]:
Image('/home/ribli/DT40/SNP/varscan/igv_snapshots/ds14__hits/1-15680078.png')
Out[31]:
In [15]:
i=2
print best_pos_list[i]
print
run_pileup(1,best_pos_list[i].split('\t')[1])
In [32]:
Image('/home/ribli/DT40/SNP/varscan/igv_snapshots/ds14__hits/1-23852983.png')
Out[32]:
In [16]:
i=3
print best_pos_list[i]
print
run_pileup(1,best_pos_list[i].split('\t')[1])
And another one¶
- 1:35282831
- strange again, two snp just disappeared next to each other, I think indel is too far to be the problem
In [33]:
Image('/home/ribli/DT40/SNP/varscan/igv_snapshots/ds14__hits/1-35282831.png')
Out[33]:
In [18]:
i=5
print best_pos_list[i]
print
run_pileup(1,best_pos_list[i].split('\t')[1])
In [34]:
Image('/home/ribli/DT40/SNP/varscan/igv_snapshots/ds14__hits/1-37277423.png')
Out[34]:
In [20]:
i=7
print best_pos_list[i]
print
run_pileup(1,best_pos_list[i].split('\t')[1])
And another one¶
- 1:52021699
- indel alignment problem, bit strange, why so homogenoues in samples? Local realignement makes it homogeneous?
In [35]:
Image('/home/ribli/DT40/SNP/varscan/igv_snapshots/ds14__hits/1-52021699.png')
Out[35]:
In [21]:
i=8
print best_pos_list[i]
print
run_pileup(1,best_pos_list[i].split('\t')[1])
In [36]:
Image('/home/ribli/DT40/SNP/varscan/igv_snapshots/ds14__hits/1-60425138.png')
Out[36]:
In [22]:
i=9
print best_pos_list[i]
print
run_pileup(1,best_pos_list[i].split('\t')[1])
In [37]:
Image('/home/ribli/DT40/SNP/varscan/igv_snapshots/ds14__hits/1-77729845.png')
Out[37]:
In [23]:
i=10
print best_pos_list[i]
print
run_pileup(1,best_pos_list[i].split('\t')[1])
Conclusions:¶
Some positions are still misterious.¶
Some positions only beacause of fluctuation.¶
- Thats good news, thats what we can very easily avoid looking at a robust germline.
- Different from the high coverage R1!