search

samtools / bcftools

This is the official development repository for BCFtools. See installation instructions and other documentation here http://samtools.github.io/bcftools/howtos/install.html
http://samtools.github.io/bcftools/
Other
633 stars 241 forks source link

Apart from the header, bcftools mpileup did not output any other content #2189

Closed dongaoqian closed 1 month ago

dongaoqian commented 1 month ago

Hello, I encountered an issue while using bcftools. The problem is that when I use bcftools mpileup, the output file only contains a header starting with '#', but no other content. I have referred to issues #1336 and #2069, checked my BAM file, and tried removing the -v option, but none of these solutions have resolved my problem. Here is the command I am using: /share/home/stu_dongaoqian/miniconda3/envs/GAEP/bin/bcftools mpileup -r Chr10A -a SP -f ../genome.fa var_calling_output.markdup.bam | /share/home/stu_dongaoqian/miniconda3/envs/GAEP/bin/bcftools call -m -f GQ -o var_calling_output.Chr10A.vcf My bcftools version is 1.9. The content of my output file is as follows (with some chromosomes and contigs omitted in the middle): image image

Looking forward to your response!

dongaoqian commented 1 month ago

I'm sorry, I don't know why I can't upload picture. I can only display my BAM file in this way. The output of checking my BAM file according to #2069 is as follows: (samtools view var_calling_output.markdup.bam Chr10A | less)

SRR23891451.193666211 113 Chr10A 1 0 151M = 60885302 60885301 AAGATATGGTAATCGAATATAAGTTTTGAGCAGGAATGTACATATATGAGGCCATGATGAAGATCCGGGAAAATGGGTGTCTCATGGCCGCCCACCAGTTTTGGTAAGTCCCAAAGATACTTGGAACATATATAGCAACTGAGATCATGAC A<-A<AAAJFAF7FFJAJJJFJFF-JFFJFJJJJJF<JJAFJJJFFAFAJ7FA-AFJJAJJFF-AAAFFJJFJJFFFJFAA-JJJJJFA7JJ<JJFJFAAFF-JAF<<FJJJJJJFJJJJJJFFJJJJJJFJFJAFFJFJFFJJJJFFF-A NM:i:0 MD:Z:151 AS:i:151 XS:i:151 XA:Z:Chr10A,-60885302,151M,0;Chr10E,-35404132,151M,2;Chr10B,-55594745,151M,3;Chr10C,-57919232,151M,3;Chr10F,-25923861,151M,3; MQ:i:0 MC:Z:151M ms:i:5363 SRR23891451.182525645 129 Chr10A 16 0 151M = 60885317 60885301 AATATAAGTTTTGAGCAGGAATGTACATATATGAGGCCATGATGAAGATCCGGGAAAATGGGTGTCTCATGGCCGCCCACCAGTTTTGGTAAGTCCCAAAGATACTTGGAACATATATAGCAACTGAGATCATGACTTGTCATATTTGAAC AAFFFJJJJJJJJJFJJJJAJJJJJJJJJJJJJJJJJJJJFFJJJJJJJJJJFJJFJFFJJJ<AAJJJJJJJJJFFJFJFJFJFJFFJJJJ-JJFJJ<FJF-7-FFJFFJFFJJJ<JAJ--F<<<JJFJFFFAFJJFAFF7-<7<AJ7AF- NM:i:0 MD:Z:151 AS:i:151 XS:i:151 XA:Z:Chr10A,+60885317,151M,0;Chr10E,+35404147,151M,2;Chr10B,+55594760,151M,3;Chr10F,+25923876,151M,3;Chr10C,+57919247,151M,3; MQ:i:0 MC:Z:151M ms:i:5478 SRR23891451.182337362 65 Chr10A 20 0 151M = 60885321 60885301 TAAGTTTTGAGCAGGAATGTACATATATGAGGCCATGATGAAGATCCGGGAAAATGGGTGTCTCATGGCCGCCCACCAGTTTTGGTAAGTCCCAAAGATACTTGGAACATATATAGCAACTGAGATCATGACTTGTCATATTTGAACATTT AAFFFJJJJJJJJJJJJFJJJJFJJJJJJAJJJJJJJJJJJJJFJJJJJJAJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJFJJJJJFF7FAJJJAFJJFJJJFJJJFJJJJFJFJ7F<AJF7<FJFFJJJ<<FAJAA<AJJFF- NM:i:0 MD:Z:151 AS:i:151 XS:i:151 XA:Z:Chr10A,+60885321,151M,0;Chr10E,+35404151,151M,2;Chr10F,+25923880,151M,3;Chr10C,+57919251,151M,3;Chr10B,+55594764,151M,3; MQ:i:0 MC:Z:151M ms:i:5828

pd3 commented 1 month ago

This requires a test case to reproduce and debug the problem. This script can be used to create a small slice of the bam and the reference https://github.com/pd3/mpileup-tests/blob/main/misc/create-bam-test

dongaoqian commented 1 month ago

Thank you for your reply. This is a small portion of the BAM file I extracted using this script. slice.tar.gz

pd3 commented 1 month ago

Thank you for the test case. The cause of the problem is that all reads are silently filtered because they have the PAIRED flag set, but not PROPER_PAIR. The program has the option to work with these anyway

-A, --count-orphans     Include anomalous read pairs, i.e. with flag PAIRED but not PROPER_PAIR

I just pushed a commit which clarifies what anomalous read pairs are, as this was not described anywhere.

Hopefully this resolves the problem!

dongaoqian commented 1 month ago

Thank you for your help,It has indeed solved my problem!