Hi all,
I have a pipeline which I have written in python to do SNP calling and annotation. The pipeline links Samtools, SomaticSniper, and SnpSift/SnpEff together to produce annotated VCF files for human cancer studies. The inputs are two sorted indexed BAM files, and a human reference fasta file, in this case GRCh38.
The first step in my pipeline splits the bams by chromosome. Then each tumor/normal chromosomal pair is run through SomaticSniper as described in the documentation. The pipeline has worked on various BAM files in the past, but my current data set is not working properly. When I run my current BAMs through SomaticSniper, the program runs, exits cleanly (no errors) and produces empty output files every time.
Here is the command I am using to call SomaticSniper:
bam-somaticsniper -F vcf -f /path/GRCh38/hg38.fa /path/TumorBams/Tumor_B02077_CHR1.bam /path/NormalBams/Normal_B02077_CHR1.bam /path/SomSniper/B02077_CHR1Snipe.vcf
Here is a head from one of the BAM files (converted to SAM for display)
@HD VN:1.0 GO:none SO:none
@SQ SN:CHR5 LN:181538259 UR:>CHR5
@SQ SN:CHRM LN:16569 UR:>CHRM
@SQ SN:CHR22 LN:50818468 UR:>CHR22
@SQ SN:CHR14 LN:107043718 UR:>CHR14
@SQ SN:CHR1 LN:248956422 UR:>CHR1
@SQ SN:CHR4 LN:190214555 UR:>CHR4
@SQ SN:CHR6 LN:170805979 UR:>CHR6
@SQ SN:CHR10 LN:133797422 UR:>CHR10
@SQ SN:CHR12 LN:133275309 UR:>CHR12
@SQ SN:CHRX LN:156040895 UR:>CHRX
@SQ SN:CHR19 LN:58617616 UR:>CHR19
@SQ SN:CHR7 LN:159345973 UR:>CHR7
@SQ SN:CHR18 LN:80373285 UR:>CHR18
@SQ SN:CHR13 LN:114364328 UR:>CHR13
@SQ SN:CHR2 LN:242193529 UR:>CHR2
@SQ SN:CHR15 LN:101991189 UR:>CHR15
@SQ SN:CHR3 LN:198295559 UR:>CHR3
@SQ SN:CHR21 LN:46709983 UR:>CHR21
@SQ SN:CHR11 LN:135086622 UR:>CHR11
@SQ SN:CHR20 LN:64444167 UR:>CHR20
@SQ SN:CHR17 LN:83257441 UR:>CHR17
@SQ SN:CHR16 LN:90338345 UR:>CHR16
@SQ SN:CHR9 LN:138394717 UR:>CHR9
@SQ SN:CHR8 LN:145138636 UR:>CHR8
@SQ SN:CHRY LN:57227415 UR:>CHRY
r 0 CHR10 10001 25 76M * 0 0 CTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAGCCCTAACCCAAACCCTAACCCTGA @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
r 0 CHR10 10001 25 76M * 0 0 CTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAA @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
r 0 CHR10 10001 25 76M * 0 0 CTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCAAACCCAAACCCAAACCCTAA @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
r 0 CHR10 10001 25 76M * 0 0 CTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCCAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAA @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
r 0 CHR10 10001 25 76M * 0 0 CTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTCA @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
r 16 CHR10 10001 25 76M * 0 0 CTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAA @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
r 16 CHR10 10001 25 76M * 0 0 CTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTACCCCTAACCCTAACCCTAA @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
r 16 CHR10 10001 25 76M * 0 0 CTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTCACCCTAACCCTAACACTCA @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
r 0 CHR10 10002 25 76M * 0 0 TAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCCAAC @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
r 0 CHR10 10002 25 76M * 0 0 TAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAAC @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
You may notice that the quality scores are all @. That is because these are "homemade" BAM files (which I think maybe causing the issues). Any thoughts/ideas as to why I am getting empty files would be appreciated.
Thanks,
Matt
I'm curious about the capitlisation of CHR in the bam file. I've not used GRCh38 but for hg19, the chromosomes do not have a CHR prefixes, but chr prefixes, and GRCh37 has no prefix at all, just the chromosome numbers. Does your reference FASTA file have CHR entries for chromosome names?
It does, we had to modify the the GRCh38 reference to add the capital CHR. Figured that was easier than modifying all of the SAM files we had already created.
Your bam file is not sorted. Did you try sorting it and then run?
What makes you think its not sorted?
The header file says so.
SO
should saysorted
instead ofnone
.Interesting, I sorted the BAM and then used reheader, which would explain why the header says "none". Samtools and the other tools I have used still accept this file as sorted (they don't tell me it needs to be sorted), so they must check some other aspect of the BAM to determine if it is sorted or not.
Good catch though Jordan. That might have been causing the issues my group was experiencing. We ended up rebuilding all of the files from scratch to get it them to work.