This is a test version of Biostars. For the public version, visit https://www.biostars.org.
bash loop to count variants using vcftools

I can run the below vcftools command to count individual variant types in vcf.gz files

zcat /home/cmccabe/Desktop/vcf/file1.vcf.gz | vcf-annotate --fill-type | grep -oP "TYPE=\w+" | sort | uniq -c > /home/cmccabe/Desktop/vcf/file1_variant_counts.bed

However when I try a bash loop the command does not run at all and I can not seem to figure out why, I think it looks right? Thank you :).

for f in /home/cmccabe/Desktop/vcf/*.vcf.gz ; do
     bname=`basename $f`
     pref=${bname%%.vcf.gz
     zcat /home/cmccabe/Desktop/vcf/$f | vcf-annotate --fill-type | grep -oP "TYPE=\w+" | sort | uniq -c > /home/cmccabe/Desktop/vcf/${pref}_variant_counts.bed
done

I think I see it now, I missed a closing brace in the pref=.... Thank you :).

vcftools bash
• 3,928 views • link
updated by Ram Ram • written by bioguy24 bioguy24

Exactly. The closing brace was the problem. Nice use of %% BTW :)

• 0 views • link
written by Ram Ram

In your loop, f=/home/cmccabe/Desktop/vcf/*.vcf.gz; in your zcat command you open

"/home/cmccabe/Desktop/vcf/home/cmccabe/Desktop/vcf/*.vcf.gz"
• 0 views • link
updated by Ram Ram • written by michael.ante michael.ante

It must be zcat $f

• 0 views • link
written by venu venu

2 answers 

for f in /home/cmccabe/Desktop/vcf/*.vcf.gz
> do
> zcat $f | vcf-annotate --fill-type | grep -oP "TYPE=\w+" | sort | uniq -c > "$f"_variant_counts.bed
> done

Then change all the file names with rename command.

Ex: 

rename "s/.vcf.gz_variant_counts.bed/.bed/" *.bed
• 1 views • link
updated by Ram Ram • written by venu venu

Thank you all :).

• 0 views • link
written by bioguy24 bioguy24

Try to use GNU Parallel for these kind of tasks:

parallel --jobs <int> "zcat {} |  vcf-annotate --fill-type | grep -oP \"TYPE=\\w+\" | sort | uniq -c > {.}_variant_counts.bed" ::: /home/cmccabe/Desktop/vcf/*.vcf.gz
• 0 views • link
updated by Ram Ram • written by GouthamAtla GouthamAtla

Log in to answer this question.

More posts like this
  • gatk vqsr --recal-file error
    written by bioguy24 bioguy24
    In GATK 4.0.6.0 I am getting the below error using the command: gatk VariantRecalibrator \ -R /home/cmccabe/Desktop/fastq/resources/ucsc.hg19.fasta \ -V NA12878.vcf.gz \ --resource hapmap,known=false,training=true,truth=true,prior=15.0:/home/cmccabe/Desktop/fastq/resources/hapmap_3.3.hg19.sites.vcf.gz \ --resource …
  • error using bqsr and incompatible contigs
    written by bioguy24 bioguy24
    Using the below command in GATK 4.0.6.0, I am getting: A USER ERROR has occurred: Input files reference and features have incompatible contigs: No overlapping …
  • GATK MarkDuplicates output and bam
    written by bioguy24 bioguy24
    The below command executes and produces the dup_metrics below, I am trying to interpret them and don't really understand the metrics and in the bam …
  • multithread support for gatk
    written by bioguy24 bioguy24
    The below command is meant to use bqsr with multithreading: gatk BaseRecalibrator --num_cpu_threads_per_data_thread=8 \ -I /home/cmccabe/Desktop/fastq/NA12878_marked_duplicates.bam \ -R /home/cmccabe/Desktop/fastq/resources/ucsc.hg19.fasta \ --known-sites /home/cmccabe/Desktop/fastq/fastq/resources/dbsnp_138.b37.vcf \ --known-sites /home/cmccabe/Desktop/fastq/resources/Mills_and_1000G_gold_standard.indels.hg19.vcf …
  • using bash variable in picard MarkDuplicates command
    written by bioguy24 bioguy24
    In the below `bash` the two variables `$bname` and `$sample` are extracted correctly however when I pass them to the `java` command `$bname` is set …
  • Bash to match and store numerical prefix to update files
    written by bioguy24 bioguy24
    In the bash below the unique headers of each vcf.gz are stored in a text file with the same name. There can be multiple vcf.gz …
  • add gene name from bed file to vcf in vcftools
    written by bioguy24 bioguy24
    I am trying to add gene information to a vcf file using vcftools: vcftools annotate \ -a genes.bed.gz \ -c CHROM,FROM,TO,GENE \ -h <(echo '##INFO=<ID=GENE,Number=1,Type=String,Description="Gene …
  • coverageBed with -g option ?
    written by bioguy24 bioguy24
    I am running a coverageBed and getting the below error. ``` cmccabe@DTV-A5211QLM:~/Desktop/NGS$ coverageBed -d -sorted -a /home/cmccabe/Desktop/NGS/bed/bedtools/xgen_targets_sorted.bed -b /home/cmccabe/Desktop/NGS/pool_I_090215/IonXpress_008_150902_newheader.bam > /home/cmccabe/Desktop/NGS/pool_I_090215/IonXpress_008_150902_output.txt ERROR: Database file /home/cmccabe/Desktop/NGS/pool_I_090215/IonXpress_008_150902_newheader.bam …
  • coverage bed error but produces output file
    written by bioguy24 bioguy24
    The below command produces a file (674 MB). ``` coverageBed \ -d \ -sorted \ -a /home/cmccabe/Desktop/NGS/bed/sorted_unix_5column_xgen_targets.bed \ -b /home/cmccabe/Desktop/NGS/pool_I_090215/IonXpress_008_150902_newheader.bam \ > /home/cmccabe/Desktop/NGS/pool_I_090215/IonXpress_008_150902_output.txt ``` Output …
  • run picard on all bam files in directory
    written by bioguy24 bioguy24
    I am trying to run a bash loop on all bam files in a directory and having a bit of trouble. The bash so far …