diff --git a/compile.sh b/compile.sh new file mode 100755 index 0000000..9e55095 --- /dev/null +++ b/compile.sh @@ -0,0 +1,2 @@ +#!/bin/sh +gcc -Wall -pthread vcf_to_sfs.c -lm -lz -std=c99 -Wextra -o vcf_to_sfs diff --git a/vcf_to_sfs.c b/vcf_to_sfs.c new file mode 100644 index 0000000..8ea4c8d --- /dev/null +++ b/vcf_to_sfs.c @@ -0,0 +1,62 @@ +# include +# include +# include +#include +#include + +bool StartsWith(const char *a, const char *b) +{ + if(strncmp(a, b, strlen(b)) == 0) return 1; + return 0; +} + +void slice_str(const char * str, char * buffer, size_t start, size_t end) +{ + size_t j = 0; + for ( size_t i = start; i <= end; ++i ) { + buffer[j++] = str[i]; + } + buffer[j] = 0; +} + + +# define LL 8192 /* line length maximum */ + +int main ( int argc, char *argv[] ){ + if ( argc < 3) { + printf("Need 2 args!\n"); + return 1; + } + gzFile fp; + char line[LL]; + char delim[] = "\t"; + fp = gzopen( argv[1], "r" ); + + gzgets( fp, line, LL ); + while ( ! gzeof( fp ) ){ + int k = 0; + if ( StartsWith(line, "##") || ( StartsWith(line, "#") ) || (strstr(line, "./.:.") != NULL)){ + gzgets( fp, line, LL ); + continue; + } + + char *vcf_field = strtok(line, delim); + while(vcf_field != NULL){ + k++; + if (k > 9) { + const size_t len = strlen(vcf_field); + char buffer[len + 1]; + //printf("'%s'\n", ptr); + slice_str(vcf_field, buffer, 0, 0); + printf("%s ", buffer); + } + vcf_field = strtok(NULL, delim); + } + // printf("%s", line ); + // loads the next line + gzgets( fp, line, LL ); + } + + gzclose( fp ); + return 0; +} diff --git a/vcf_to_sfs.py b/vcf_to_sfs.py index 34ca424..ae1df7d 100755 --- a/vcf_to_sfs.py +++ b/vcf_to_sfs.py @@ -71,10 +71,11 @@ with gzip.open(sys.argv[1], "rb") as inputgz: for k in set(snp_genotypes): allele_counts[snp_genotypes.count(k)] = k allele_counts_list.append(snp_genotypes.count(k)) - if folded : - for al in range(polyallelic-1): + if folded and len(ALT) >= 2: + for al in range(len(ALT)-1): SFS_values[min(allele_counts_list)-1] += 1/len(ALT) allele_counts_list.remove(min(allele_counts_list)) + else: + SFS_values[min(allele_counts_list)-1] += 1 line = inputgz.readline() print(SFS_values) -print(polycount)