| Previous changeset 5:d5725351bf22 (2015-11-02) Next changeset 7:b617219a70b5 (2015-11-06) |
|
Commit message:
new script updates |
|
modified:
vcf_snp.py |
| b |
| diff -r d5725351bf22 -r 3a352cb57117 vcf_snp.py --- a/vcf_snp.py Mon Nov 02 19:42:54 2015 -0500 +++ b/vcf_snp.py Fri Nov 06 15:12:17 2015 -0500 |
| [ |
| @@ -1,13 +1,14 @@ #!/usr/bin/env python +##!/usr/local/anaconda/bin/python ######################################################################################### # # # Name : vcf_snp.py # -# Version : 0.2 # +# Version : 0.3 # # Project : extract snp from vcf # # Description : Script to exctract snps # # Author : Brigida Rusconi # -# Date : November 02 2015 # +# Date : November 06 2015 # # # ######################################################################################### @@ -30,40 +31,54 @@ parser.add_argument('-o', '--output', help="snp tab") parser.add_argument('-s', '--snp_table', help="vcf") +parser.add_argument('-p', '--positions', help="positions") + args = parser.parse_args() output_file = args.output input_file = args.snp_table +input_file2=args.positions #------------------------------------------------------------------------------------------ #read in file as dataframe df =read_csv(input_file,sep='\t', dtype=object) +df3=read_csv(input_file2,sep='\t', dtype=object, names=['#CHROM','POS']) +#pdb.set_trace() +#get the positions that are missing +if df.index.size!= df3.index.size: + df4=df3[~df3.POS.isin(df.POS)] +#pdb.set_trace() + df4=df4.set_index('POS') + df=df.set_index('POS') + df2=concat([df,df4], axis=1) +#if there is no information for a loci the position is dropped so it fills the nan with N + df5=df2.iloc[:, 2:4].fillna('N') + df5=df5.reindex(df3.POS) +else: + df5=df.iloc[:,3:5] #------------------------------------------------------------------------------------------ -# only columns with qbase and refbase in table -#count_qbase=list(df.columns.values) -#qindexes=[] -#for i, v in enumerate(count_qbase): -# if 'ALT' in v: -# qindexes.append(i) -df2=df.iloc[:,3:5] + #pdb.set_trace() #------------------------------------------------------------------------------------------ +#replaces the . with the nucleotide call of the reference also deals with multiallelic states calling them N ref_list=[] -for i in range(0,df2.index.size): - if df2.iloc[i,1]==".": - ref_list.append(df2.iloc[i,0][0]) +for i in range(0,df5.index.size): + if df5.iloc[i,1]==".": + ref_list.append(df5.iloc[i,0][0]) + elif len(df5.iloc[i,1])>1: + ref_list.append('N') else: - ref_list.append(df2.iloc[i,1][0]) + ref_list.append(df5.iloc[i,1][0]) #pdb.set_trace() # ##------------------------------------------------------------------------------------------ # #save file with output name for fasta -o option and removes header and index with open(output_file,'w') as output: - output.write(df.columns.values[4] + '\t' + ''.join([str(i) for v,i in enumerate(ref_list)])) + output.write('ALT' + '\t' + ''.join([str(i) for v,i in enumerate(ref_list)])) ##------------------------------------------------------------------------------------------ \ No newline at end of file |