SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
problems with bam file and mpileup Danielbenitezr Bioinformatics 5 03-06-2013 09:07 AM
Corrupt mpileup file? JackieBadger Bioinformatics 3 10-22-2012 12:38 PM
VarScan input error with mpileup file ??? JackieBadger Bioinformatics 3 10-22-2012 09:15 AM
large BAM, but very small mpileup file mrfox Bioinformatics 5 10-16-2012 12:16 PM
mpileup on whole genome bam file adrian Bioinformatics 3 08-14-2012 02:42 AM

Reply
 
Thread Tools
Old 08-26-2015, 03:59 AM   #1
Momina
Member
 
Location: uk

Join Date: Aug 2015
Posts: 19
Lightbulb mpileup file

hi all
I am trying to create a mpileup file so I can get vcf file which I can run on IGV and can use for SNP calling
but i am facing following problems:
i gave following command for mpile up file:
source samtools
samtools mpileup -E -uf reference.fa file.sorted.bam > file.mpileup
it gave output:
[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000
and computer stucks here

when i restart, a mpileup is created
then i try to creat vcf file from that mpileup by giving following command:
source bcftools
bcftools view -cg file.mpileup > file.vcf
error occurs:
Error: Could not parse --min-ac g

can anyone help me in that i am stuck at SNP calling
Momina is offline   Reply With Quote
Old 08-26-2015, 05:02 AM   #2
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,061
Default

What version of samtools are you using?
GenoMax is offline   Reply With Quote
Old 08-26-2015, 05:13 AM   #3
Momina
Member
 
Location: uk

Join Date: Aug 2015
Posts: 19
Default

samtools-0.1.19
Momina is offline   Reply With Quote
Old 08-26-2015, 05:23 AM   #4
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,061
Default

Pardon me for the obvious question but are you using your own file names in place of "reference.fa/file.sorted.bam"?
GenoMax is offline   Reply With Quote
Old 08-26-2015, 05:36 AM   #5
Momina
Member
 
Location: uk

Join Date: Aug 2015
Posts: 19
Default

yes
this is what I am commanding:
samtools mpileup -E -uf Triticum_aestivum.IWGSC1.0+popseq.28.dna_sm.genome.fa 606.sorted.bam > 606.mpileup
Momina is offline   Reply With Quote
Old 08-26-2015, 05:42 AM   #6
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,061
Default

What do you mean by "computer stuck here"? Depending on how large your 606.sorted.bam file it can take a while for the pileup file to be created.

Can you post the output of

Code:
$ ls -lh 606*
GenoMax is offline   Reply With Quote
Old 08-26-2015, 05:46 AM   #7
Momina
Member
 
Location: uk

Join Date: Aug 2015
Posts: 19
Default

this is what I ma doing:
source samtools-0.1.19
samtools mpileup -E -uf Triticum_aestivum.IWGSC1.0+popseq.28.dna_sm.genome.fa 604.sorted.bam > 604.mpileup
output:
[mpileup] 1 samples in 1 input files
<mpileup> Set max per-file depth to 8000
after this I press Ctrl C
it starts again this time it does not stuck
but when I write:
source bcftools-1.2
[hussainm@NBI-HPC analysis]$ bcftools view -cg 604.mpileup > 604.vcf
following error occurs
Error: Could not parse --min-ac g
Momina is offline   Reply With Quote
Old 08-26-2015, 05:51 AM   #8
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,061
Default

It sounds like you are not allowing the first command to complete its work (samtools mpileup). Like I said before it will take some time (could be a few hours if you have a large bam file) for this command to finish. Until that command completes normally you should not start the bcftools process.
GenoMax is offline   Reply With Quote
Old 08-26-2015, 05:55 AM   #9
Momina
Member
 
Location: uk

Join Date: Aug 2015
Posts: 19
Default

[hussainm@NBI-HPC analysis]$ ls -lh 606*
-rwxrwx--- 1 hussainm JIC_a1 12G Aug 21 12:36 606.bam
-rwxrwx--- 1 hussainm JIC_a1 32M Aug 21 14:08 606.sorted.bai
-rwxrwx--- 1 hussainm JIC_a1 7.5G Aug 21 13:30 606.sorted.bam
-rwxrwx--- 1 hussainm JIC_a1 413 Aug 21 15:12 606.sorted_flagstat.txt
Momina is offline   Reply With Quote
Old 08-26-2015, 06:03 AM   #10
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,061
Default

Looks like you have more than one file that starts with 60* names.

In any case if they are all similar size (e.g. 606.sorted.bam is 7.5G) then it is going to take some time for the mpileup to do its work.

Are you using a job scheduler and/or working on a cluster or is this a standalone server? If this is a standalone server wait for the command prompt to return (do not cancel the job by ctrl+c) after you issue this command:

Code:
[hussainm@NBI-HPC analysis]$ samtools mpileup -E -uf Triticum_aestivum.IWGSC1.0+popseq.28.dna_sm.genome.fa 604.sorted.bam > 604.mpileup

Last edited by GenoMax; 08-26-2015 at 06:06 AM.
GenoMax is offline   Reply With Quote
Old 08-26-2015, 06:04 AM   #11
Momina
Member
 
Location: uk

Join Date: Aug 2015
Posts: 19
Default

ok I will wait now for few hours
and following is the output:
[hussainm@NBI-HPC analysis]$ ls -lh 606*
-rwxrwx--- 1 hussainm JIC_a1 12G Aug 21 12:36 606.bam
-rwxrwx--- 1 hussainm JIC_a1 32M Aug 21 14:08 606.sorted.bai
-rwxrwx--- 1 hussainm JIC_a1 7.5G Aug 21 13:30 606.sorted.bam
-rwxrwx--- 1 hussainm JIC_a1 413 Aug 21 15:12 606.sorted_flagstat.txt
Momina is offline   Reply With Quote
Old 08-26-2015, 06:08 AM   #12
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,061
Default

Looks like you are using a high performance compute cluster (based on your command prompt). Get some help from your IT support to see if you can submit this job via a job scheduler (so you do not need to keep an interactive session going). If you are running this job on the "head node" that is probably not a good idea anyway.
GenoMax is offline   Reply With Quote
Old 08-26-2015, 06:11 AM   #13
Momina
Member
 
Location: uk

Join Date: Aug 2015
Posts: 19
Default

ya I have 7 files and they all are near 12 -13 GB
and this was a standalone server
now I ma thinking of submitting job
yes I can submit job
and I do that too
I will let you know what there response
following is output of bjobs:
JOBID USER STAT QUEUE FROM_HOST EXEC_HOST JOB_NAME SUBMIT_TIME
62876 hussain RUN NBI-Prod25 v0292.hpccl ncn-768-01. *5.mpileup Aug 26 15:08
Momina is offline   Reply With Quote
Old 08-26-2015, 06:17 AM   #14
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,061
Default

Great. You can start jobs for other files in parallel so that all 7 get done at the same time. Once those jobs finish go on to the bcftools for the 7 mpileup files.

For future reference you can combine these two processes into one by using unix pipes. See example here: http://samtools.sourceforge.net/mpileup.shtml
GenoMax is offline   Reply With Quote
Old 08-26-2015, 06:23 AM   #15
Momina
Member
 
Location: uk

Join Date: Aug 2015
Posts: 19
Default

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)?
from where var.raw.bcf comes?
and what is vcfutils.pl??
if you can tell me this with one good example?
Momina is offline   Reply With Quote
Old 08-26-2015, 08:11 AM   #16
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,061
Default

Quote:
Originally Posted by Momina View Post
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.
Quote:
from where var.raw.bcf comes?
Output of samtoools mpileup is passed to bcftools command which ultimately produces the "var.raw.bcf" file.
Quote:
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.

Code:
$ samtools mpileup -E -uf reference.fa file.sorted.bam | bcftools view -cg - > file.vcf
GenoMax is offline   Reply With Quote
Old 08-26-2015, 08:47 AM   #17
Momina
Member
 
Location: uk

Join Date: Aug 2015
Posts: 19
Default

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
Momina is offline   Reply With Quote
Old 08-26-2015, 08:55 AM   #18
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,061
Default

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
GenoMax is offline   Reply With Quote
Old 08-26-2015, 09:34 AM   #19
Momina
Member
 
Location: uk

Join Date: Aug 2015
Posts: 19
Default

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?
Momina is offline   Reply With Quote
Old 08-26-2015, 09:55 AM   #20
GenoMax
Senior Member
 
Location: East Coast USA

Join Date: Feb 2008
Posts: 7,061
Default

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)
Code:
$ bcftools call -cv -Ov file.mpileup > calls.vcf
OR if you want to get a compressed calls file

Code:
$ bcftools call -cv -Oz file.mpileup > calls.vcf.gz
GenoMax is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 07:45 AM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO