In this notebook i will try to explore VarScan 'SNP' positions in a R1 untreated sample: Sample 4¶
- Im trying to look at alignment context information using IGV, and all sample information using mpileup
In [1]:
#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 Sample4_Sample1.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 [4]:
#select some very unlikely fals positives ( p< 0.0001)
best_pos_cmd='''
head Sample4_Sample1.vsc.snp -n 100000 | 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 [5]:
#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 [6]:
#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 [2]:
Image('/home/ribli/DT40/SNP/varscan/igv_snapshots/sample4_hits/1-4938485.png')
Out[2]:
In [7]:
i=0
print best_pos_list[i]
print
run_pileup(1,best_pos_list[i].split('\t')[1])
In [3]:
Image('/home/ribli/DT40/SNP/varscan/igv_snapshots/sample4_hits/1-4940483.png')
Out[3]:
In [8]:
i=1
print best_pos_list[i]
print
run_pileup(1,best_pos_list[i].split('\t')[1])
In [4]:
Image('/home/ribli/DT40/SNP/varscan/igv_snapshots/sample4_hits/1-5501200.png')
Out[4]:
In [9]:
i=4
print best_pos_list[i]
print
run_pileup(1,best_pos_list[i].split('\t')[1])
In [5]:
Image('/home/ribli/DT40/SNP/varscan/igv_snapshots/sample4_hits/1-5696972.png')
Out[5]:
In [10]:
i=5
print best_pos_list[i]
print
run_pileup(1,best_pos_list[i].split('\t')[1])
In [6]:
Image('/home/ribli/DT40/SNP/varscan/igv_snapshots/sample4_hits/1-5966182.png')
Out[6]:
In [11]:
i=6
print best_pos_list[i]
print
run_pileup(1,best_pos_list[i].split('\t')[1])
In [7]:
Image('/home/ribli/DT40/SNP/varscan/igv_snapshots/sample4_hits/1-6008344.png')
Out[7]:
In [12]:
i=7
print best_pos_list[i]
print
run_pileup(1,best_pos_list[i].split('\t')[1])
In [8]:
Image('/home/ribli/DT40/SNP/varscan/igv_snapshots/sample4_hits/1-7860898.png')
Out[8]:
In [13]:
i=8
print best_pos_list[i]
print
run_pileup(1,best_pos_list[i].split('\t')[1])
And another one¶
- 1:8672283
- Deletion, sometimes misaligned -> SNP
- why is it consequently misaligned in 2 samples, and not is others? A nearby SNP? The local realigner make it like this?
In [9]:
Image('/home/ribli/DT40/SNP/varscan/igv_snapshots/sample4_hits/1-8672283.png')
Out[9]:
In [14]:
i=9
print best_pos_list[i]
print
run_pileup(1,best_pos_list[i].split('\t')[1])
In [10]:
Image('/home/ribli/DT40/SNP/varscan/igv_snapshots/sample4_hits/1-10936409.png')
Out[10]:
In [15]:
i=10
print best_pos_list[i]
print
run_pileup(1,best_pos_list[i].split('\t')[1])
Conclusions:¶
Almost all indel errors¶
- Why are the indel errors so homoegenous? Local realignement make them homogeneous?
All positions could be filtered looking at all the samples¶
- More robust germline escapes false positive calls because of fluctuation, or germline LOH-s.
- At strange positions, more than one samples shows the stange behaviour.