Dear Friends,

I mostly script in shell/bash as am comfortable with it and also use a lot of linux commands and Entrez utilities for various big file analysis; however, when I stared working on VCF files, dealing with million rows and columns I find issues in operating these files with bash. For instance, below is an example script to count number of SNPs for samples(TCGA...barcode...) in a VCF file:

    for i in {10..10000}
    do
    snp_count=$(zcat xxx.vcf.gz | grep -v "##" | awk "{print \$$i}" | grep -v "TCGA" | grep "0/1\|1/1" | wc -l)
    sample=$(zcat xxx.vcf.gz | grep -v "##" | awk "{print \$$i}" | grep "TCGA")
    echo -e "$sample\t$snp_count\n"
    done

    Output:
    TCGA...barcode1     60123
    TCGA...barcode2      45245
    .
    .


The script works but it has been running for 2 days now, and will take time. I would really appreciate if you could share/suggest on the best ways to deal with VCF files using bash/shell or if another programming language could make a difference? However, if possible, I would prefer bash scripting ways to tackle VCF files.

Thanks,
DK

Here is oneliner using bioalcidaejdk  http://lindenb.github.io/jvarkit/BioAlcidaeJdk.html

    $ wget  -O - "ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ALL.chrY.phase3_integrated_v2a.20130502.genotypes.vcf.gz" |\
    gunzip -c |\
    java -jar dist/bioalcidaejdk.jar -e 'stream().flatMap(V->V.getGenotypes().stream()).filter(G->!(G.isHomRef()||G.isNoCall())).map(G->G.getSampleName()).collect(Collectors.groupingBy(Function.identity(),Collectors.counting())).forEach((k,v)->println(k+" "+v));' -F VCF
    (...)
    NA20852 441
    HG02489 1233
    HG02008 112
    NA11919 905
    HG02002 119
    HG01398 133
    HG03577 837
    HG02009 1217
    NA18549 819

detail:

* `stream().` get the stream of variants
* `flatMap(V->V.getGenotypes().stream()).` convert  it to a stream of genotypes
* `filter(G->!(G.isHomRef()||G.isNoCall())).` remove the hom_ref and the nocall 
* `map(G->G.getSampleName()).` extract the sample name from the genotype
* `collect(Collectors.groupingBy(Function.identity(),Collectors.counting())).` count the number of distinct names
* `forEach((k,v)->println(k+" "+v));` print it

