Question: filtering pileup file with samtools for SomaticSniper
gravatar for TriS
24 months ago by
United States, Buffalo
TriS3.0k wrote:

hi all

I'm working with matched tumor-normal samples and I want to call somatic variants using Somatic Sniper

So far I got the main function working with

bam-somaticsniper -Q 40 -G -L -f genome/hg19.fa Results/kayrotypicT.bam Results/kayrotypicN.bam SomSniperResults.txt

and now I need to filter the results.

I'm following what the SomaticSniper guide says and running the pipeline following this guide on Synapse.

first question is:

the SomaticSniper guide says explicitly NOT to use mpileup

You will also need to generate a samtools pileup (not mpileup) indel file.

so, an old post explains the difference..while I think samtools replaced pileup with mpileup, so, since mpileup should output pileup still, is mpileup fine then?

since the example uses a -vcf output I'm creating the pileup and doing the filtering in the following way

# create pileup file 
samtools mpileup -A -B -vuf genomes/hg19.fa Results/kayrotypic.bam > SomSniperResults/var.raw.pileup 

a snapshot of my output looks like this:

##FORMAT=<ID=PL,Number=G,Type=Integer,Description="List of Phred-scaled genotype likelihoods">
#CHROM  POS     ID      REF     ALT     QUAL    FILTER  INFO    FORMAT  RS-01926692_JC590_RS-01925823
chr1    10008   .       A       <X>     0       .       DP=1;I16=0,1,0,0,33,1089,0,0,0,0,0,0,0,0,0,0;QS=1,0;MQ0F=1      PL      0,3,4
chr1    10009   .       A       <X>     0       .       DP=1;I16=0,1,0,0,30,900,0,0,0,0,0,0,1,1,0,0;QS=1,0;MQ0F=1       PL      0,3,4
chr1    10010   .       C       <X>     0       .       DP=1;I16=0,1,0,0,33,1089,0,0,0,0,0,0,2,4,0,0;QS=1,0;MQ0F=1      PL      0,3,4

however, when i call varFilter test.pileup | awk '$6>=20'

I get the following error:

Argument "DP=1;I16=0,1,0,0,33,1089,0,0,0,0,0,0,0,0,0,0;QS=1,0;MQ0F=1" isn't numeric in numeric lt (<) at line 117, <> line .

which comes from the fact that the 7th column should be numeric while what I have is the info from the vcf file...

I also tried not to output a vcf file and tried to use bcftools as:

samtools mpileup -A -B -guf /genomes/hg19.fa chr20_bam.bam > test.pileup
bcftools/bcftools view -Hgc test.pileup

[bcf_sync] incorrect number of fields (0 != 5) at 0:0

I looked at the post here but it didn't solve my problem since i'm not piping and I am testing only one small chr (chr20) from the bam file


so I'm not totally sure how to get it to work.


-- txt was edited with some updated info --

pileup snp samtools somaticsniper • 1.2k views
ADD COMMENTlink modified 10 months ago by guangxujin0 • written 24 months ago by TriS3.0k
gravatar for TriS
23 months ago by
United States, Buffalo
TriS3.0k wrote:

alright, not sure ppl will have the same problem but here's what I think can be the solution:

1) I had IT install the newer samtools 1.2 + bcftools 1.2 on the server

2) using the following command line:

samtools mpileup -A -B -uf genomes/hg19.fa Results_Norm/norm.bam | bcftools call -c | varFilter -Q 20 | awk 'NR > 55 {print}'> indel_pileup_Norm.pileup
samtools mpileup -A -B -uf genomes/hg19.fa Results_Tum/Tum.bam | bcftools call -c | varFilter -Q 20 | awk 'NR > 55 {print}' > indel_pileup_Tum.pileup

I was able to get (hopefully) the correct files for normal and tumor samples to then use in SomaticSniper

ADD COMMENTlink modified 23 months ago • written 23 months ago by TriS3.0k

It works well. Thanks.

ADD REPLYlink written 10 months ago by guangxujin0
gravatar for guangxujin
10 months ago by
guangxujin0 wrote:

It works well. Thanks.

ADD COMMENTlink written 10 months ago by guangxujin0
Please log in to add an answer.


Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1670 users visited in the last hour