|
0
|
1 #!/usr/bin/env python
|
|
6
|
2 ##!/usr/local/anaconda/bin/python
|
|
0
|
3
|
|
|
4 #########################################################################################
|
|
|
5 # #
|
|
|
6 # Name : vcf_snp.py #
|
|
6
|
7 # Version : 0.3 #
|
|
0
|
8 # Project : extract snp from vcf #
|
|
|
9 # Description : Script to exctract snps #
|
|
|
10 # Author : Brigida Rusconi #
|
|
6
|
11 # Date : November 06 2015 #
|
|
0
|
12 # #
|
|
|
13 #########################################################################################
|
|
|
14
|
|
2
|
15 # If the position is identical to the reference it does not print the nucleotide. I have to retrieve it from the ref column.
|
|
0
|
16
|
|
|
17
|
|
|
18 #------------------------------------------------------------------------------------------
|
|
|
19
|
|
|
20
|
|
|
21 import argparse, os, sys, csv, IPython
|
|
|
22 import pandas
|
|
|
23 import pdb
|
|
|
24 from pandas import *
|
|
4
|
25
|
|
0
|
26 #------------------------------------------------------------------------------------------
|
|
|
27
|
|
|
28
|
|
|
29 #output and input file name to give with the script
|
|
|
30 parser = argparse.ArgumentParser()
|
|
|
31
|
|
|
32 parser.add_argument('-o', '--output', help="snp tab")
|
|
|
33 parser.add_argument('-s', '--snp_table', help="vcf")
|
|
6
|
34 parser.add_argument('-p', '--positions', help="positions")
|
|
|
35
|
|
0
|
36
|
|
|
37
|
|
|
38 args = parser.parse_args()
|
|
|
39 output_file = args.output
|
|
|
40 input_file = args.snp_table
|
|
6
|
41 input_file2=args.positions
|
|
0
|
42 #------------------------------------------------------------------------------------------
|
|
|
43
|
|
|
44
|
|
|
45 #read in file as dataframe
|
|
|
46 df =read_csv(input_file,sep='\t', dtype=object)
|
|
6
|
47 df3=read_csv(input_file2,sep='\t', dtype=object, names=['#CHROM','POS'])
|
|
|
48 #pdb.set_trace()
|
|
|
49 #get the positions that are missing
|
|
|
50 if df.index.size!= df3.index.size:
|
|
|
51 df4=df3[~df3.POS.isin(df.POS)]
|
|
|
52 #pdb.set_trace()
|
|
|
53 df4=df4.set_index('POS')
|
|
|
54 df=df.set_index('POS')
|
|
|
55 df2=concat([df,df4], axis=1)
|
|
|
56 #if there is no information for a loci the position is dropped so it fills the nan with N
|
|
|
57 df5=df2.iloc[:, 2:4].fillna('N')
|
|
|
58 df5=df5.reindex(df3.POS)
|
|
|
59 else:
|
|
|
60 df5=df.iloc[:,3:5]
|
|
0
|
61
|
|
|
62 #------------------------------------------------------------------------------------------
|
|
|
63
|
|
6
|
64
|
|
0
|
65 #pdb.set_trace()
|
|
|
66
|
|
|
67 #------------------------------------------------------------------------------------------
|
|
6
|
68 #replaces the . with the nucleotide call of the reference also deals with multiallelic states calling them N
|
|
2
|
69 ref_list=[]
|
|
6
|
70 for i in range(0,df5.index.size):
|
|
|
71 if df5.iloc[i,1]==".":
|
|
|
72 ref_list.append(df5.iloc[i,0][0])
|
|
|
73 elif len(df5.iloc[i,1])>1:
|
|
|
74 ref_list.append('N')
|
|
2
|
75 else:
|
|
6
|
76 ref_list.append(df5.iloc[i,1][0])
|
|
2
|
77 #pdb.set_trace()
|
|
|
78 #
|
|
|
79 ##------------------------------------------------------------------------------------------
|
|
|
80 #
|
|
0
|
81 #save file with output name for fasta -o option and removes header and index
|
|
|
82 with open(output_file,'w') as output:
|
|
6
|
83 output.write('ALT' + '\t' + ''.join([str(i) for v,i in enumerate(ref_list)]))
|
|
2
|
84 ##------------------------------------------------------------------------------------------ |
