Question: csaw - Htslib import errors and missing reads/chromosomes
0
gravatar for reskejak
7 months ago by
reskejak10
Michigan State University
reskejak10 wrote:

I'm working through csaw for paired-end ATAC-seq data (see my other post for more background on my general workflow), and I don't think my BAM files are being properly read into R. I have used numerous files that are coordinate-sorted (also have tried name-sorted) and indexed (via either samtools or Rsamtools) and I keep getting these Htslib errors and warnings upon csaw::windowCounts or csaw::regionCounts:

[W::bam_hdr_read] bgzf_check_EOF: No such file or directory
[E::bgzf_read] Read block operation failed with error -1 after zd of zu bytes

As the command continues, the first warning appears many times. The files ultimately end up loading, but I believe they are only partially read, as none of my macs2 peaks on chr1 show up (0 counts across all samples) and chrX and chrY do not appear to be read at all.

Example for

> standard.chr <- paste0("chr", c(1:19, "X", "Y"))
> param <- readParam(max.frag=1000, pe="both", 
                     discard=blacklist, restrict=standard.chr)
> binned <- windowCounts(pe.bams, bin=TRUE, width=10000, param=param)
> binned@rowRanges@seqnames

factor-Rle of length 211698 with 19 runs
Lengths:     4 12711 11876 11694 11717 12151 10069  9504  9176  8749  5822 17839 15645 15223 14806 14637 14194 12570  3311
values :  chr1 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chr19 chr2  chr3  chr4 chr5  chr6  chr7  chr8  chr9
Levels(21): chr1 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 ... chr3 chr4 chr5 chr6 chr7 chr8 chr9 chrX chrY

As you can see, no 10kb bins were called for chrX or chrY, and only 4 were called for chr1.

I'm pretty certain the BAM file is not corrupted as the headers and EOF exist, and samtools quickcheck elicits no errors.

Example of header via samtools view -H file.sorted.bam

@HD     VN:1.5  SO:coordinate
@SQ     SN:chr1 LN:195471971
@SQ     SN:chr10        LN:130694993
@SQ     SN:chr11        LN:122082543
@SQ     SN:chr12        LN:120129022
@SQ     SN:chr13        LN:120421639
@SQ     SN:chr14        LN:124902244
@SQ     SN:chr15        LN:104043685
@SQ     SN:chr16        LN:98207768
@SQ     SN:chr17        LN:94987271
@SQ     SN:chr18        LN:90702639
@SQ     SN:chr19        LN:61431566
@SQ     SN:chr1_GL456210_random LN:169725
@SQ     SN:chr1_GL456211_random LN:241735
@SQ     SN:chr1_GL456212_random LN:153618
@SQ     SN:chr1_GL456213_random LN:39340
@SQ     SN:chr1_GL456221_random LN:206961
@SQ     SN:chr2 LN:182113224
@SQ     SN:chr3 LN:160039680
@SQ     SN:chr4 LN:156508116
@SQ     SN:chr4_GL456216_random LN:66673
@SQ     SN:chr4_JH584292_random LN:14945
@SQ     SN:chr4_GL456350_random LN:227966
@SQ     SN:chr4_JH584293_random LN:207968
@SQ     SN:chr4_JH584294_random LN:191905
@SQ     SN:chr4_JH584295_random LN:1976
@SQ     SN:chr5 LN:151834684
@SQ     SN:chr5_JH584296_random LN:199368
@SQ     SN:chr5_JH584297_random LN:205776
@SQ     SN:chr5_JH584298_random LN:184189
@SQ     SN:chr5_GL456354_random LN:195993
@SQ     SN:chr5_JH584299_random LN:953012
@SQ     SN:chr6 LN:149736546
@SQ     SN:chr7 LN:145441459
@SQ     SN:chr7_GL456219_random LN:175968
@SQ     SN:chr8 LN:129401213
@SQ     SN:chr9 LN:124595110
@SQ     SN:chrX LN:171031299
@SQ     SN:chrX_GL456233_random LN:336933
@SQ     SN:chrY LN:91744698
@SQ     SN:chrY_JH584300_random LN:182347
@SQ     SN:chrY_JH584301_random LN:259875
@SQ     SN:chrY_JH584302_random LN:155838
@SQ     SN:chrY_JH584303_random LN:158099
@SQ     SN:chrUn_GL456239       LN:40056
@SQ     SN:chrUn_GL456367       LN:42057
@SQ     SN:chrUn_GL456378       LN:31602
@SQ     SN:chrUn_GL456381       LN:25871
@SQ     SN:chrUn_GL456382       LN:23158
@SQ     SN:chrUn_GL456383       LN:38659
@SQ     SN:chrUn_GL456385       LN:35240
@SQ     SN:chrUn_GL456390       LN:24668
@SQ     SN:chrUn_GL456392       LN:23629
@SQ     SN:chrUn_GL456393       LN:55711
@SQ     SN:chrUn_GL456394       LN:24323
@SQ     SN:chrUn_GL456359       LN:22974
@SQ     SN:chrUn_GL456360       LN:31704
@SQ     SN:chrUn_GL456396       LN:21240
@SQ     SN:chrUn_GL456372       LN:28664
@SQ     SN:chrUn_GL456387       LN:24685
@SQ     SN:chrUn_GL456389       LN:28772
@SQ     SN:chrUn_GL456370       LN:26764
@SQ     SN:chrUn_GL456379       LN:72385
@SQ     SN:chrUn_GL456366       LN:47073
@SQ     SN:chrUn_GL456368       LN:20208
@SQ     SN:chrUn_JH584304       LN:114452
@PG     ID:bowtie2      PN:bowtie2      VN:2.2.6  ...
# rest of header...

And here's an example read via samtools view file.sorted.bam | head:

NS500653:255:H2YWNAFXY:3:11506:4441:16945       99      chr1    3000222 27      75M     =       3000302 15 GTCTAGGAAGTTGTCCATTTCATCCAGGTTTTCCTGGTTTTTTTTTAGTATAGCCTTTCATAGTAGAATCTGATG    AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEE/EEEEEAEAEEEEAE/E<A<AAA6EEEEA    MD:Z:65A9       PG:Z:MarkDuplicates     XG:i:0  NM:i:1  XM:i:1  XN:i:0  XO:i:0  AS:i:-4 XS:i:-41        YS:i:-5YT:Z:CP

Any suggestions for how to fix these issues when loading BAMs into csaw? I'm working on Windows, so maybe it would help to try this out on the Centos cluster.

library atac-seq chip-seq bam csaw • 411 views
ADD COMMENTlink modified 7 months ago • written 7 months ago by reskejak10
0
gravatar for reskejak
7 months ago by
reskejak10
Michigan State University
reskejak10 wrote:

UPDATE: ANSWER

I have not exactly identified the cause of this Rhtslib error, but I tried many alternative and repeated methods to generating these BAMs from scratch to no avail. This was ultimately remedied by using csaw on a Centos cluster with >100 gb RAM.

I have not tried submitting this command as a job with specific RAM designations, but my working theory is that 8 or 16 gb of RAM (on my local Windows machines) was not sufficient to fully load this experiment. This is further supported by the fact that another, much smaller data set could be successfully loaded on my local machine.

In conclusion: make sure your machine has sufficient RAM to load large experiments in csaw.

ADD COMMENTlink written 7 months ago by reskejak10
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 966 users visited in the last hour