This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Issues generating bcf file with samtools

I am using samtools to make a bcf file but I am having some issues. I am not really sure what is going on here and I was hoping to get some help! To generate my bam file I used trimmomatic, bwa, and picard. I am not too sure what is going on here because none of my other files are giving me any issues.

So the issue I am having here is that when I run the following command, the bcf file that is generated never gets larger over time. It just stays at 0 bytes. I am not sure why this is happening? I had no issues with my other samples.

samtools mpileup --redo-BAQ --min-BQ 30 --per-sample-mF -C 50 --output-tags DP,DV,DPR,INFO/DPR,DP4,SP -u -f Felis_catus.Felis_catus_9.0.dna_sm.toplevel.fa --BCF C04893_L001.sorted.dedup.bam > C04893.bcf

the output from samtools flagstat from the file that is giving me issues

690090027 + 0 in total (QC-passed reads + QC-failed reads)
7861265 + 0 secondary
0 + 0 supplementary
52950324 + 0 duplicates
688355595 + 0 mapped (99.75% : N/A)
682228762 + 0 paired in sequencing
341114381 + 0 read1
341114381 + 0 read2
657843484 + 0 properly paired (96.43% : N/A)
679649948 + 0 with itself and mate mapped
844382 + 0 singletons (0.12% : N/A)
19554958 + 0 with mate mapped to a different chr
16483263 + 0 with mate mapped to a different chr (mapQ>=5)

An example of a file that is working correctly

619420497 + 0 in total (QC-passed reads + QC-failed reads)
8721163 + 0 secondary
0 + 0 supplementary
89789840 + 0 duplicates
617477073 + 0 mapped (99.69% : N/A)
610699334 + 0 paired in sequencing
305349667 + 0 read1
305349667 + 0 read2
577735322 + 0 properly paired (94.60% : N/A)
608184428 + 0 with itself and mate mapped
571482 + 0 singletons (0.09% : N/A)
27685608 + 0 with mate mapped to a different chr
23986111 + 0 with mate mapped to a different chr (mapQ>=5)
assembly genome alignment next-gen
• 222 views • link
written by williamsbrian5064 williamsbrian5064

What is the problem?

• 0 views • link
written by ATpoint ATpoint

Sorry, when I go to generate the bcf file, there is just no output. I the size of the file stays at 0

• 0 views • link
written by williamsbrian5064 williamsbrian5064

samtools is deprecated for calling variants : use bcftools mpileup.

do you have any output without ' --redo-BAQ --min-BQ 30 --per-sample-mF -C 50 '

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

I can give bcftools a try but I am still not sure why I was able to get all my other samples to run.

When I run the following

samtools mpileup  --output-tags DP,DV,DPR,INFO/DPR,DP4,SP -u -f Felis_catus.Felis_catus_9.0.dna_sm.toplevel.fa --BCF C04893_L001.sorted.dedup.bam > C04893.bcf

I still get no output.

• 0 views • link
written by williamsbrian5064 williamsbrian5064

Why have you got the --BCF flag? - you don't need that. Pileup format cannot be BCF, in any case, so you should not have .bcf as your output file extension.

• 0 views • link
written by Kevin Blighe Kevin Blighe

I think it can be in bcf format? Link "-o, --output FILE Write pileup or VCF/BCF output to FILE, rather than the default of standard output"

Regardless, I am not sure why this method has worked with other samples but not this particular sample. I am really trying to figure out what the issue with this one sample. I don't think it is samtools/bcftools because other samples finished fine. The flagstat seems normal enough to me

• 0 views • link
written by williamsbrian5064 williamsbrian5064

Did you not try BCFtools?

I have samtools v1.9 and it says:

Note that using "samtools mpileup" to generate BCF or VCF files is now deprecated. To output these formats, please use "bcftools mpileup" instead.

• 0 views • link
written by Kevin Blighe Kevin Blighe

I noticed that there was an extra number added to the read group. It said 

@RG\tID:C04893_S47_0001\tSM:C04893\tPL:illumina\tLB:C04893_lib1\tPU:barcode_C048937

but it should be 

@RG\tID:C04893_S47_0001\tSM:C04893\tPL:illumina\tLB:C04893_lib1\tPU:barcode_C04893

Do you think that could be the issue?

• 0 views • link
written by williamsbrian5064 williamsbrian5064 
bwa mem -M -t 16 -R '@RG\tID:C04893_S47_0001\tSM:C04893\tPL:illumina\tLB:C04893_lib1\tPU:barcode_C048937' Felis_catus.Felis_catus_9.0.dna_sm.toplevel.fa C04893_R1_001.trim.paired.fastq.gz C04893_R2_001.trim.paired.fastq.gz > C04893_001_aligned.sam

java -jar -Xmx8g -Djava.io.tmpdir=tmp_C04893_L001_sort picard.jar SortSam INPUT=C04893_001_aligned.sam OUTPUT=C04893_L001.sorted.bam SORT_ORDER=coordinate MAX_RECORDS_IN_RAM=5000000

java -jar -Xmx8g -Djava.io.tmpdir=tmp_C04893_L001_dedup picard.jar MarkDuplicates INPUT=C04893_L001.sorted.bam OUTPUT=C04893_L001.sorted.dedup.bam CREATE_INDEX=true METRICS_FILE=C04893_L001.metrics MAX_FILE_HANDLES_FOR_READ_ENDS_MAP=1000
• 0 views • link
written by williamsbrian5064 williamsbrian5064

don't use picard for sorting when you can use 'samtools sort'

• 0 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

You seems to use a quite old version of samtools. I would suggest upgrade first and use bcftools mpileup as Pierre suggested.

Are you aware that the output of mpileup is not your final variant calling?

• 0 views • link
written by finswimmer finswimmer

0 answers

No answers yet.

Log in to answer this question.

More posts like this