arq5x / bedtools2

bedtools - the swiss army knife for genome arithmetic
MIT License
939 stars 287 forks source link

genomecov -d -pc -ibam showing no coverage for every position #1080

Open gcabebe opened 9 months ago

gcabebe commented 9 months ago

I have bam files that I've obtained through STAR alignment. Here is the first 10 lines of one of them:

GWNJ-1013:671:GW2304233104th:2:1126:15501:17879 0       chromosome      1       255     11S139M *       0       0       ATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTCACGCTACAACAGG  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFF  NH:i:1  HI:i:1  AS:i:137        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1126:15130:18427 0       chromosome      1       255     11S139M *       0       0       ATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTCACGCTACAACAGG  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF::FFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF  NH:i:1  HI:i:1  AS:i:137        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1206:9272:29857  0       chromosome      1       255     28S108M14S      *       0       0       GGCGGTTTGGCGTGCGTATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGAGATCGGAAGAGCG  FFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF  NH:i:1  HI:i:1  AS:i:106        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1266:30734:10692 0       chromosome      1       255     24S126M *       0       0       GTTTGGCGTGCGTATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTC  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFF  NH:i:1  HI:i:1  AS:i:124        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1266:31412:33411 0       chromosome      1       255     24S126M *       0       0       GTTTGGCGTGCGTATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTC  FFFFFFFFFFFFFFFF:FFFFFFFFFFFF:FFFF::FFFFFFF:FFFFFFF:FFFFFFFFFFFFFFFF:FFFF,FFFFFFF:FFFFFFFFF:FFFFFFFFFFFFFFF:FFF::FFFFFFFF,FFFFFFFF:FFFFFFFFFFF:F,FFFFF  NH:i:1  HI:i:1  AS:i:124        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1330:21938:10708 0       chromosome      1       255     16S134M *       0       0       TGCGTATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTCACGCTACA  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFF,FFFFFFFFFFFFFFFFF:F  NH:i:1  HI:i:1  AS:i:132        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1367:17074:4726  0       chromosome      1       255     11S82M  *       0       0       ATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGG   FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF  NH:i:1   HI:i:1  AS:i:80 nM:i:0
GWNJ-1013:671:GW2304233104th:2:1367:11939:11584 0       chromosome      1       255     11S82M  *       0       0       ATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGG   FFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFF:F:FFFFFFFFFFFFFFFFFFF  NH:i:1   HI:i:1  AS:i:80 nM:i:0
GWNJ-1013:671:GW2304233104th:2:1372:24740:30843 0       chromosome      1       255     16S48M  *       0       0       TGCGTATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAA        FFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF        NH:i:1  HI:i:1  AS:i:47 nM:i:0
GWNJ-1013:671:GW2304233104th:2:1401:11107:19319 0       chromosome      1       255     5S144M1S        *       0       0       ATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTCACGCTACAACAGGAATTGA  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFF  NH:i:1  HI:i:1  AS:i:142        nM:i:0

This is the script I'm using for the genomecov command:

BAM_DIR='./sorted_bam_files'

# Loop over each BAM file in the directory
for file in "$BAM_DIR"/*_high_mapq_reads.sorted.bam
do
    echo "Calculating coverage for $file..."
    # Get the basename of the file
    base=$(basename "$file")
    # Calculate coverage
    bedtools genomecov -d -pc -ibam $file > "${base%.sorted.bam}.coverage"
done

This is the output I keep getting. I've tried this with at least 30 other bam files but the output is always the same.

chromosome      1       0
chromosome      2       0
chromosome      3       0
chromosome      4       0
chromosome      5       0
chromosome      6       0
chromosome      7       0
chromosome      8       0
chromosome      9       0
chromosome      10      0

Any ideas on why this keeps happening? I don't think there's something wrong with my bam files because I was able to use them for DGE analyses and checked to see if most reads were getting unmapped but they weren't. Any help is appreciated.

arq5x commented 9 months ago

The FLAG column in the examples you show is 255, which means that the read and its mate are unmapped (plug in 255 here: https://broadinstitute.github.io/picard/explain-flags.html). If all alignments are "unmapped" (I don't know why), then there is no coverage, technically. genomecov requires the reads to be mapped.

gcabebe commented 8 months ago

Someone on StackExchanged answered my question here:

It was because I included the argument -pc, which is only for paired-end reads. The second column on the bam file I showed here is all 0, which indicates single-end reads. So I removed that argument and was able to get the tables I needed, which look like this:

chromosome      1       169
chromosome      2       169
chromosome      3       171
chromosome      4       171
chromosome      5       173
chromosome      6       178
chromosome      7       178
chromosome      8       180
chromosome      9       212
chromosome      10      212