search

varlociraptor / varlociraptor

Flexible, uncertainty-aware variant calling with parameter free filtration via FDR control.
https://varlociraptor.github.io
GNU General Public License v3.0
118 stars 16 forks source link

Unable to open a SAM file, varlociraptor 2.3.0 #151

Closed nh13 closed 3 years ago

nh13 commented 3 years ago

varlociraptor 2.3.0

command

$ varlociraptor preprocess variants tmp.fasta --bam tmp.normal.sam --candidates tmp.vcf
Error: unable to open SAM/BAM/CRAM file at tmp.normal.sam
$ file tmp.normal.sam 
tmp.normal.sam: Sequence Alignment/Map (SAM), with header version 1.6
>chr1
ggctAcggannnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
nnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnnn
##fileformat=VCFv4.2
##ALT=<ID=NON_REF,Description="Represents any non-reference allele.">
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Depth per allele.">
##FORMAT=<ID=FREQ,Number=1,Type=Float,Description="Non-ref frequency.">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype quality.">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype string.">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Phred scaled genotype likelihoods.">
##FORMAT=<ID=PS,Number=1,Type=Integer,Description="Phase set.">
##INFO=<ID=AC,Number=A,Type=Integer,Description="Alternate allele counts in genotypes.">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Depth across all samples.">
##contig=<ID=chr1,length=1000>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  normal  tumor
chr1    5   .   A   C   .   .   .   GT:AD:FREQ:PS   0|0:50,0:0.0:5  0/1:80,20:0.2
##fileformat=VCFv4.2
##ALT=<ID=NON_REF,Description="Represents any non-reference allele.">
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Depth per allele.">
##FORMAT=<ID=FREQ,Number=1,Type=Float,Description="Non-ref frequency.">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype quality.">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype string.">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Phred scaled genotype likelihoods.">
##FORMAT=<ID=PS,Number=1,Type=Integer,Description="Phase set.">
##INFO=<ID=AC,Number=A,Type=Integer,Description="Alternate allele counts in genotypes.">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Depth across all samples.">
##contig=<ID=chr1,length=1000>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  normal  tumor
chr1    5   .   A   C   .   .   .   GT:AD:FREQ:PS   0|0:50,0:0.0:5  0/1:80,20:0.2
(testing-varlociraptor) nils-fg:varlociraptor nhomer$ cat tmp.normal.sam 
@HD VN:1.6  SO:unsorted GO:none
@SQ SN:chr1 LN:1000
@RG ID:normal   SM:normal
0000    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0001    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0002    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0003    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0004    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0005    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0006    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0007    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0008    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0009    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0010    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0011    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0012    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0013    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0014    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0015    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0016    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0017    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0018    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0019    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0020    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0021    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0022    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0023    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0024    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0025    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0026    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0027    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0028    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0029    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0030    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0031    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0032    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0033    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0034    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0035    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0036    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0037    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0038    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0039    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0040    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0041    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0042    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0043    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0044    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0045    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0046    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0047    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0048    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0049    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0050    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0051    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0052    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0053    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0054    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0055    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0056    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0057    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0058    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0059    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0060    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0061    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0062    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0063    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0064    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0065    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0066    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0067    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0068    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0069    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0070    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0071    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0072    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0073    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0074    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0075    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0076    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0077    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0078    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0079    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0080    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0081    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0082    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0083    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0084    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0085    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0086    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0087    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0088    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0089    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0090    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0091    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0092    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0093    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0094    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0095    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0096    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0097    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0098    0   chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
0099    16  chr1    1   60  9M  *   0   0   GGCTACGGA   ?????????   RG:Z:normal
nh13 commented 3 years ago

varlociraptor 3.0.1

Caused by:
    unable to open SAM/BAM/CRAM index for tmp.normal.sam

well I don't have a SAM index ;P, it's a SAM file, so I was hoping it would just scan the file?

dlaehnemann commented 3 years ago

This is actually rust-htslib behaviour when using a bam::IndexedReader:

https://github.com/rust-bio/rust-htslib/blob/4ee5dd16c5c67b198782dbbc5bc31f72f073ee22/src/bam/mod.rs#L591-L595

As it will fetch reads piling up above candidate sites and those may be apart by quite some distance (or you might have parallelized preprocessing over different parts of the genome), the default behavior is index-based skipping to those sites instead of streaming to them.

So would you want this to be a feature request to allow streaming instead of index-based skipping in some way? Or is the most recent error message (and then simply creating the index) fine as a solution for you?

nh13 commented 3 years ago

Either the error message should say "file must be indexed" or "support streaming". Either is fine, but sometimes I use SAM files to produce test cases (unindexed) and debug odd variants, so that'd be great.

johanneskoester commented 3 years ago

Since varlociraptor starts from candidate variants, I think we won't add support for just streaming through the bam. It will always have to jump, and just for testcases, I don't want specialized code. In particular, since it is pretty easy to convert into a BAM+index. Also, note that because of varlociraptors estimation of sampling bias related variables (see #154), it is anyway usually not optimal to have a handcrafted testcase. We have derived all of the ~60 included testcases from synthetic WGS, and real data.

johanneskoester commented 3 years ago

But I agree about the error message. It should definitely be improved. Hence, leaving it open for now.

johanneskoester commented 3 years ago

xref: https://github.com/rust-bio/rust-htslib/pull/303