Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

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

Open
gcabebe opened this issue Feb 15, 2024 · 2 comments
Open

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

gcabebe opened this issue Feb 15, 2024 · 2 comments

Comments

@gcabebe
Copy link

gcabebe commented Feb 15, 2024

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
Copy link
Owner

arq5x commented Feb 15, 2024

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
Copy link
Author

gcabebe commented Mar 14, 2024

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants