Bam file doesn't show on IGV and pileup VCF has only header
0
0
Entering edit mode
8.8 years ago
umn_bist ▴ 390

I have PE RNA-seq fastq files that I trimmed, aligned, and deduped, split N_CIGAR, recalibrated. Once I got my bam file, I wanted to make sure it was done correctly, but I see nothing in IGV. I zoomed in, still nothing. I printed the header and it looked fine.

samtools mpileup -Buf /Documents/gencode/GRCh38.p5.genome.fa \
    -l /Documents/plist_hg38.bed /Documents/Projects/120529_0165_L007_tstaidsr.bam |\
    bcftools call -Amv > /Documents/Projects/120529_0165_L007_tstaidsr.vcf

samtools view -H /Documents/Projects/120529_0165_L007_tstaidsr.bam

Output:

@HD    VN:1.5    GO:none    SO:coordinate
@SQ    SN:chr1    LN:248956422
@SQ    SN:chr2    LN:242193529
@SQ    SN:chr3    LN:198295559
@SQ    SN:chr4    LN:190214555
@SQ    SN:chr5    LN:181538259
@SQ    SN:chr6    LN:170805979
@SQ    SN:chr7    LN:159345973
@SQ    SN:chr8    LN:145138636
@SQ    SN:chr9    LN:138394717
@SQ    SN:chr10    LN:133797422
@SQ    SN:chr11    LN:135086622
@SQ    SN:chr12    LN:133275309
@SQ    SN:chr13    LN:114364328
@SQ    SN:chr14    LN:107043718
@SQ    SN:chr15    LN:101991189
@SQ    SN:chr16    LN:90338345
@SQ    SN:chr17    LN:83257441
@SQ    SN:chr18    LN:80373285
@SQ    SN:chr19    LN:58617616
@SQ    SN:chr20    LN:64444167
@SQ    SN:chr21    LN:46709983
@SQ    SN:chr22    LN:50818468
@SQ    SN:chrX    LN:156040895
@SQ    SN:chrY    LN:57227415
@SQ    SN:chrM    LN:16569
@SQ    SN:GL000008.2    LN:209709
@SQ    SN:GL000009.2    LN:201709
@SQ    SN:GL000194.1    LN:191469
@SQ    SN:GL000195.1    LN:182896
@SQ    SN:GL000205.2    LN:185591
@SQ    SN:GL000208.1    LN:92689
@SQ    SN:GL000213.1    LN:164239
@SQ    SN:GL000214.1    LN:137718
@SQ    SN:GL000216.2    LN:176608
@SQ    SN:GL000218.1    LN:161147
@SQ    SN:GL000219.1    LN:179198
@SQ    SN:GL000220.1    LN:161802
@SQ    SN:GL000221.1    LN:155397
@SQ    SN:GL000224.1    LN:179693
@SQ    SN:GL000225.1    LN:211173
@SQ    SN:GL000226.1    LN:15008
... ... ...

As always, thank you for your time and help.
RNA-Seq • 3.3k views
ADD COMMENT •
1
Entering edit mode

Please follow the steps outlined in the answer to this question. Briefly, pipe the output of samtools view into head to get the first few alignments. Then zoom to those regions in IGV. Let me know if you still don't see anything.

ADD REPLY •
0
Entering edit mode

Thanks for your reply. That was actually a post I was looking at. The above output was with

samtools view -H
ADD REPLY •
1
Entering edit mode

Right, so that's just showing your header. We want to look at the first few alignments without the header. Instead, type this:

samtools view [your.bam] | head
ADD REPLY •
0
Entering edit mode

Hey, Dan. This is what I get:

DGR4KXP1:165:C0UMTACXX:7:1107:11496:46307    163    chr1    12691    1    47M    =    12782    138    CACTCTTGTGAGTGTCCCCAGTGTTGCAGAGGTGAGAGGAGAGTAGA    BF:==?DF?DBD@CBDBBBDDBDDADCEECBC2AFEDD?CEEE@CFE    BD:Z:PPQOWWVTUUQRPPNOOLLOOMNLKJNOONOONPNQRUXXYXYQOPQ    MD:Z:47    PG:Z:MarkDuplicates    RG:Z:120529_UNC16-SN851_0165_BC0UMTACXX_CAGATC_L007    NH:i:3    BI:Z:MMMNVWTRPNMLNMMKJIHMKLKKIIIMLKLKONOPRUVXWWWPNLN    HI:i:1    jI:B:i,-1    NM:i:0    jM:B:c,-1    nM:i:0    AS:i:92
DGR4KXP1:165:C0UMTACXX:7:1107:11496:46307    83    chr1    12782    1    47M    =    12691    -138    CGTCTCCTGTCTCCTGGAGAGGCTTCGATGCCCCTCCACACCCTCTT    ?EDDBBDDDBDCBCDBCCCDBCEBC>CDDBBBBDEDECD@A<@=<DB    BD:Z:QQQRXYYWWUSPPPOONNNNOOLLONMMONLMPPPRSSWWUYYQOPP    MD:Z:47    PG:Z:MarkDuplicates    RG:Z:120529_UNC16-SN851_0165_BC0UMTACXX_CAGATC_L007    NH:i:3    BI:Z:PLNNXVXXUURQMNPKLNKMIJHIMLKMLJIIKLMLMNPUSTVNLMM    HI:i:1    jI:B:i,-1    NM:i:0    jM:B:c,-1    nM:i:0    AS:i:92
DGR4KXP1:165:C0UMTACXX:7:1104:19969:30599    419    chr1    12963    0    22M2I23M    =    13252    336    CTTCACTCCCAGCTCAGAGCCCCCAGGCCAGGGGCCCCCAAGAAAGG    CD=>AACFCCEFEDCDDCDCBBBBDDBCBDDBBBBBBBBDADBACFD    BD:Z:PPOOWWVWWRSRQPONONOOOKKKNNONOOOOLMOQPRUYVVXOIOR    MD:Z:45    PG:Z:MarkDuplicates    RG:Z:120529_UNC16-SN851_0165_BC0UMTACXX_CAGATC_L007    NH:i:7    BI:Z:MMMMXTUSNKOMLNLMLKJJIHHHMKIJJNMKKKMOPSUZVUWMDLM    HI:i:6    jI:B:i,-1    NM:i:2    jM:B:c,-1    nM:i:0    AS:i:84
DGR4KXP1:165:C0UMTACXX:7:1104:19969:30599    339    chr1    13252    0    47M    =    12963    -336    CAGCATAGTGCTCCTGGACCAGTGATACACCCGGCACCCTGTCCTGG    FFEDCBDCDCDCBDDBDBBDDCDCCBDECBB>BCFCCCEHGB:?@@C    BD:Z:RRQQUUVXXVSPPPOOONOOLNNMKLMONLONPOQQPUWVWXYRRPP    MD:Z:47    PG:Z:MarkDuplicates    RG:Z:120529_UNC16-SN851_0165_BC0UMTACXX_CAGATC_L007    NH:i:7    BI:Z:PPMOXXWVYURQMNPKLNKMJILKKJKKLINMKKMNJLPTSVUOQMM    HI:i:6    jI:B:i,-1    NM:i:0    jM:B:c,-1    nM:i:0    AS:i:84
DGR4KXP1:165:C0UMTACXX:7:2104:20021:121056    355    chr1    13355    0    47M    =    13509    201    GCCATTGCTGCTGTGTGGAAGTTCACTCCTGCCTTTTCCTTTCCCTA    CE9@@@FEECDFFCDBDBDBDCBCDCDCBDCCBDAAACBDAABBDFD    BD:Z:PPRRWUTWWURRQNPNOONLLNLLNNLNONONOOMIKSXXVQVROQN    MD:Z:47    PG:Z:MarkDuplicates    RG:Z:120529_UNC16-SN851_0165_BC0UMTACXX_CAGATC_L007    NH:i:8    BI:Z:MMMQYSTRQOKNNMMLLKKIHLIILIJKJNNKLOMFHTVYVMVMLPO    HI:i:5    jI:B:i,-1    NM:i:0    jM:B:c,-1    nM:i:0    AS:i:92
DGR4KXP1:165:C0UMTACXX:7:2104:20021:121056    403    chr1    13509    0    47M    =    13355    -201    GCCCTTCCTTTGCTCTGCCCGCTGGAGACGGTGTTTGTCATGGGCCT    EDDD@DBDBBDCDDDDCBB>CDDBCDCC>BCDCBCFDDEGG=<=9FC    BD:Z:QOROVXYVORSQPOOONLOOONNNMMOOMNNMLGNPRSVVYUXQRPP    MD:Z:47    PG:Z:MarkDuplicates    RG:Z:120529_UNC16-SN851_0165_BC0UMTACXX_CAGATC_L007    NH:i:8    BI:Z:MLMLVXVUMTTPNMNNKINMJKMIJLJMLJJLHBKMJNPVXSUMMMM    HI:i:5    jI:B:i,-1    NM:i:0    jM:B:c,-1    nM:i:0    AS:i:92
DGR4KXP1:165:C0UMTACXX:7:2102:15520:8787    419    chr1    13554    1    47M    =    13662    154    CTGGTCTGCAGGGATCCTGCTACAAAGGTGAAACCCAGGAGAGTGTG    CD=9:@CFDEEDDCCCBDCCDBCDAADBCEBBACBBDDBBDDD@FDF    BD:Z:PPRRWXXXVUSRNONMONONOJLLKELONOMLFLONSUXXYXYQRPR    MD:Z:47    PG:Z:MarkDuplicates    RG:Z:120529_UNC16-SN851_0165_BC0UMTACXX_CAGATC_L007    NH:i:4    BI:Z:MMPNXUWUNOMKJMNLJMLILKGKI@HJMLMKCIMNUVVXWWWPOOO    HI:i:3    jI:B:i,-1    NM:i:0    jM:B:c,-1    nM:i:1    AS:i:89
DGR4KXP1:165:C0UMTACXX:7:2207:10816:90234    419    chr1    13615    0    47M    =    13748    180    GCCAGGACCCAGGCACAGGCATTAGTGCCCGTTGGAGAAAACAGGGG    CD9@@;FGCCEFDCDBDDBCDCBBDCEBBA>DBDBDDCAABCCDEDD    BD:Z:PPRRXXXWVRSRQOPNMOONOLJHLMONOLMOLLPPSTUPPUWRROO    MD:Z:47    PG:Z:MarkDuplicates    RG:Z:120529_UNC16-SN851_0165_BC0UMTACXX_CAGATC_L007    NH:i:5    BI:Z:MMMQWTVQNKOMKKNJLLIILMHIHLKJJINOLLNQRUVMMSXOMLL    HI:i:4    jI:B:i,-1    NM:i:0    jM:B:c,-1    nM:i:0    AS:i:92
DGR4KXP1:165:C0UMTACXX:7:2102:15520:8787    339    chr1    13662    1    46M1S    =    13554    -154    AATCCCGAAGAAATGGTGGGTCTTGGCCATCCGTGAGATCTTCCCAA    DEGBA>CADCAACDBCDBBCCDADBCBDCCB>CDEFBDDGGB:9=DC    BD:Z:NPQOYXVVWSKMNPNNOLNNNLLOONONMNOOOOPQQSVUUXVRNPP    MD:Z:22C23    PG:Z:MarkDuplicates    RG:Z:120529_UNC16-SN851_0165_BC0UMTACXX_CAGATC_L007    NH:i:4    BI:Z:LPOLZYVVYTHNQPLKOJJIJHIMIIJKLLNMKNLNMOOQTVTNMMM    HI:i:3    jI:B:i,-1    NM:i:1    jM:B:c,-1    nM:i:1    AS:i:89
DGR4KXP1:165:C0UMTACXX:7:2207:10816:90234    339    chr1    13748    0    47M    =    13615    -180    TGCGTGGCCGAGGGCCAGGCTTCTCACTGGGCCTCTGCAGGAGGCTG    FE@CDBCB>CDBBCBEDBCDBCCDFBECBBCBDCFFBEBBC?<=?GC    BD:Z:RSQQYYXYWURNPOOONOOLLNOMONOOLONPPOQRRUVXWXYRRPP    MD:Z:47    PG:Z:MarkDuplicates    RG:Z:120529_UNC16-SN851_0165_BC0UMTACXX_CAGATC_L007    NH:i:5    BI:Z:PPPMZVVZYUTNMLMNNKKIIJJKKLKMHJJJLLMNKOQRUWUNOMM    HI:i:4    jI:B:i,-1    NM:i:0    jM:B:c,-1    nM:i:0    AS:i:92
ADD REPLY •
1
Entering edit mode

OK, so if you load the BAM file and reference in IGV, and you zoom to chr1:12900-14000, do any alignments appear?

ADD REPLY •
0
Entering edit mode

I see nothing. I even tried zooming in as much as I can.

EDIT: My posting limit has reached max for 6 hours. I'll have to wait... I hope this edit reaches you.

Note that it's only after I try recalibrating my bases (in GATK) this happens.

java -jar -Xmx32g ${GATK} \
    -T SplitNCigarReads \
    -R "${reference}" \
    -I "${file3}" \
    -o "${file4}" \
    -rf ReassignOneMappingQuality \
    -RMQF 255 \
    -RMQT 60 \
    -U ALLOW_N_CIGAR_READS
java -jar -Xmx32g ${GATK} \
    -T BaseRecalibrator \
    -R "${reference}" \
    -I "${file4}" \
    -knownSites "${dbSNP_1}" \
    -o "${file5}"
java -jar -Xmx32g ${GATK} \
    -T BaseRecalibrator \
    -R "${reference}" \
    -I "${file4}" \
    -knownSites "${dbSNP_1}" \
    -BQSR "${file5}" \
    -o "${file6}"
java -jar -Xmx32g ${GATK} \
    -T AnalyzeCovariates \
    -R "${reference}" \
    -before "${file5}" \
    -after "${file6}" \
    -plots "${file1%_tsta.bam}_BQSR.pdf"
java -jar -Xmx32g ${GATK} \
    -T PrintReads \
    -R "${reference}" \
    -I "${file4}" \
    -BQSR "${file5}" \
    -o "${file7}"
ADD REPLY •
1
Entering edit mode

Can you please post another screenshot?

ADD REPLY •
0
Entering edit mode

My posting limit has reached max for 6 hours. I'll have to wait... I hope this edit reaches you.

Note that it's only after I try recalibrating my bases (in GATK) this happens.

java -jar -Xmx32g ${GATK} \
    -T SplitNCigarReads \
    -R "${reference}" \
    -I "${file3}" \
    -o "${file4}" \
    -rf ReassignOneMappingQuality \
    -RMQF 255 \
    -RMQT 60 \
    -U ALLOW_N_CIGAR_READS
java -jar -Xmx32g ${GATK} \
    -T BaseRecalibrator \
    -R "${reference}" \
    -I "${file4}" \
    -knownSites "${dbSNP_1}" \
    -o "${file5}"
java -jar -Xmx32g ${GATK} \
    -T BaseRecalibrator \
    -R "${reference}" \
    -I "${file4}" \
    -knownSites "${dbSNP_1}" \
    -BQSR "${file5}" \
    -o "${file6}"
java -jar -Xmx32g ${GATK} \
    -T AnalyzeCovariates \
    -R "${reference}" \
    -before "${file5}" \
    -after "${file6}" \
    -plots "${file1%_tsta.bam}_BQSR.pdf"
java -jar -Xmx32g ${GATK} \
    -T PrintReads \
    -R "${reference}" \
    -I "${file4}" \
    -BQSR "${file5}" \
    -o "${file7}"
ADD REPLY •
0
Entering edit mode

are you sure there are some reads in the zoomed region ?

ADD REPLY •
0
Entering edit mode

There are no reads in the zoomed region.

ADD REPLY •
1
Entering edit mode

my question was how do you know there are some read here ?

What is the output of

samtools view -c /Documents/Projects/120529_0165_L007_tstaidsr.bam "chr17:7677941-7678011"
ADD REPLY •
0
Entering edit mode

My mistake. I get

0

when I type

samtools view -c /Documents/Projects/120529_0165_L007_tstaidsr.bam "chr17:7677941-7678011"

EDIT: Upon closer inspection, the bam file after alignment was done correctly. It was during recalibration with GATK everything disappeared.

ADD REPLY •

Login before adding your answer.

Similar Posts
Loading Similar Posts
Traffic: 1604 users visited in the last hour
Content Search
Users
Tags
Badges
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6