[GenoMax]

Thanku
I saw this before but I didn't get it:
samtools mpileup -uf ref.fa aln1.bam aln2.bam | bcftools view -bvcg - > var.raw.bcf
bcftools view var.raw.bcf | vcfutils.pl varFilter -D100 > var.flt.vcf
aln.bam= these are bam files (probably sorted bam files)?

Correct.

from where var.raw.bcf comes?

Output of samtoools mpileup is passed to bcftools command which ultimately produces the "var.raw.bcf" file.

and what is vcfutils.pl??
if you can tell me this with one good example?

vcfutils.pl is part of the samtools package. Your admins hopefully have made it available in same path as the samtools executables.

$ samtools mpileup -E -uf reference.fa file.sorted.bam | bcftools view -cg - > file.vcf

[Momina]

hi
I submit this job:
bsub -q NBI-Prod256 "bcftools view -bvcg 607.mpileup >607.var.raw.bcf"

The output (if any) follows:

view: invalid option -- 'b'

About: VCF/BCF conversion, view, subset and filter VCF/BCF files.
Usage: bcftools view [options] <in.vcf.gz> [region1 [...]]

Output options:
-G, --drop-genotypes drop individual genotype information (after subsetting if -s option set)
-h/H, --header-only/--no-header print the header only/suppress the header in VCF output
-l, --compression-level [0-9] compression level: 0 uncompressed, 1 best speed, 9 best compression [-1]
-o, --output-file <file> output file name [stdout]
-O, --output-type <b|u|z|v> b: compressed BCF, u: uncompressed BCF, z: compressed VCF, v: uncompressed VCF [v]
-r, --regions <region> restrict to comma-separated list of regions
-R, --regions-file <file> restrict to regions listed in a file
-t, --targets [^]<region> similar to -r but streams rather than index-jumps. Exclude regions with "^" prefix
-T, --targets-file [^]<file> similar to -R but streams rather than index-jumps. Exclude regions with "^" prefix

Subset options:
-a, --trim-alt-alleles trim alternate alleles not seen in the subset
-I, --no-update do not (re)calculate INFO fields for the subset (currently INFO/AC and INFO/AN)
-s, --samples [^]<list> comma separated list of samples to include (or exclude with "^" prefix)
-S, --samples-file [^]<file> file of samples to include (or exclude with "^" prefix)
--force-samples only warn about unknown subset samples

Filter options:
-c/C, --min-ac/--max-ac <int>[:<type>] minimum/maximum count for non-reference (nref), 1st alternate (alt1), least frequent
(minor), most frequent (major) or sum of all but most frequent (nonmajor) alleles [nref]
-f, --apply-filters <list> require at least one of the listed FILTER strings (e.g. "PASS,.")
-g, --genotype [^]<hom|het|miss> require one or more hom/het/missing genotype or, if prefixed with "^", exclude sites with hom/het/missing genotypes
-i/e, --include/--exclude <expr> select/exclude sites for which the expression is true (see man page for details)
-k/n, --known/--novel select known/novel sites only (ID is not/is '.')
-m/M, --min-alleles/--max-alleles <int> minimum/maximum number of alleles listed in REF and ALT (e.g. -m2 -M2 for biallelic sites)
-p/P, --phased/--exclude-phased select/exclude sites where all samples are phased
-q/Q, --min-af/--max-af <float>[:<type>] minimum/maximum frequency for non-reference (nref), 1st alternate (alt1), least frequent
(minor), most frequent (major) or sum of all but most frequent (nonmajor) alleles [nref]
-u/U, --uncalled/--exclude-uncalled select/exclude sites without a called genotype
-v/V, --types/--exclude-types <list> select/exclude comma-separated list of variant types: snps,indels,mnps,other [null]
-x/X, --private/--exclude-private select/exclude sites where the non-reference alleles are exclusive (private) to the subset samples

[GenoMax]

That seems to be bcftools program from new version of samtools (v.1.2). You should not mix old (v.0.1.19) and new version (v.1.2) samtools programs.

If you are going to use the new version of samtools then you should use the bcftools call option: https://github.com/samtools/bcftools...pileup-calling

[Momina]

I am getting confuse now
I got mpileup file, correctly build it now I need to convert it in to file.vcf so I can call a SNP
my samtool version is v0.1.19 and bcftools version is v1.2
I tried different commands and get following outputs:
1. bcftools view -cg file.mpileup > file.vcf
output:
Error: Could not parse --min-ac g
(I don't know what is this??)
2. bcftools view -bvcg - > var.raw.bcf
output:
invalid option -b
now if I am not getting it wrong -bvcg is a command which needs a type witten next to it and it will convert the file into var.raw.bcf format
but in this command we are writing name of piledup file
so what to do?

[GenoMax]

Options for the bcftools command have changed between the old and new version.

I am not sure why your system administrators have mixed old samtools with new bcftools (unless you have sourced the programs incorrectly during startup).

Since you are using the new bcftools you will want to use something like (not sure if this will work with your file)

$ bcftools call -cv -Ov file.mpileup > calls.vcf

OR if you want to get a compressed calls file

$ bcftools call -cv -Oz file.mpileup > calls.vcf.gz

[Momina]

this is the output of command:
[hussainm@NBI-HPC analysis]$ source bcftools-1.2
[hussainm@NBI-HPC analysis]$ bcftools call -cv -Ov WT.mpileup > WT.calls.vcf
[vcf.c:759 bcf_hdr_read] invalid BCF2 magic string: only BCFv2.2 is supported.
Failed to open WT.mpileup: could not parse header

[GenoMax]

Either find bcftools from the old version (v. 0.1.19) or re-do mpileup using samtools from new.

Out of curiosity what is the outcome of these two?

1. After you "source" old samtools what is the output for the following command (if it does run or do you get a "not found" error)?

$ bcftools

Do you see anything printed to screen (specially version number)?

2. Are you able to source new samtools?

$ source samtools-1.2

[Momina]

in cluster I have only these:
[hussainm@NBI-HPC analysis]$ source samtools-
samtools-0.1.19 samtools-1.0
[hussainm@NBI-HPC analysis]$ source bcftools-1.2
now I am doing mpileup with different command:
busb -q NBI-Prod256 "samtools mpileup -g -f Triticum_aestivum.IWGSC1.0+popseq.28.dna_sm.genome.fa WT.sorted.bam 604.sorted.bam 605.sorted.bam 606.sorted.bam 607.sorted.bam 608.sorted.bam 609.sorted.bam > wheat.raw.bcf"
job is still running
after this I will go for bcftools:
bcftools view -bvcg wheat.raw.bcf > wheat-var.bcf
but question is "I have to source bcftools but version is not compatible
so what you suggest should I download new one and from where I can do that?"

[GenoMax]

You are going to run in to the same problem if you are using old samtools for this job. You should re-do mpileup using samtools from samtools-1.0, if you don't have other immediate options.

Perhaps there is a misunderstanding with the admins (they are likely not biologists and may not understand the nuances). Does bcftools-1.2 package contain corresponding samtools by any chance?

Ideally you should contact your system admins and get them to build a corresponding samtools-1.2 to keep everything in sync. You sound like a new user (and this is a cluster) so I advise against trying to download and compile samtools yourself.

[Momina]

as I convert sam to bam following old version so I have to repeat all steps?

[GenoMax]

Start with new samtools at mpileup step. Hopefully that should be enough.

[Momina]

I did that too
following output comes:
The output (if any) follows:
/usr/users/JIC_a1/hussainm/.lsbatch/1440751455.63113: line 8: samtools: command not found

[GenoMax]

Are you sourcing the new samtools (v.1.0) before running that command? What do you see when you do

$ which samtools

[GenoMax]

How do you select a particular software version to use (are you using "modules" software)?

[Momina]

this is what I did:
[hussainm@NBI-HPC analysis]$ samtools-
samtools-0.1.19 samtools-1.0 samtools-1.2
[hussainm@NBI-HPC analysis]$ samtools-1.2
[hussainm@NBI-HPC analysis]$ busb -q NBI-Prod256 "samtools mpileup -g -f Triticum_aestivum.IWGSC1.0+popseq.28.dna_sm.genome.fa WT.sorted.bam 604.sorted.bam 605.sorted.bam 606.sorted.bam 607.sorted.bam 608.sorted.bam 609.sorted.bam > wheat.raw.bcf"
I got output:
The output (if any) follows:
/usr/users/JIC_a1/hussainm/.lsbatch/1440751455.63113: line 8: samtools: command not found
I install the software on cluster by using the standard method as recommended by admin