Trying to extract .cram file read sequences into array

Question

[sacryt]

Hi I am working with .cram files and trying to use PySam to extract and view the read sequences.

This script works on a smaller CRAM file I found online:

cram = pysam.AlignmentFile(path, "rc")

for x in cram.fetch():
        print(x.query_sequence)

however, when I try and do this for a larger .cram file on the RAP on UK biobank, it gives me OSError: truncated file.

It doesn't work even if I define a region in cram.fetch() that lies within a chromosome which I can see contains several overlapping reads with IGV.

Could anyone please help with why this may be happening, or alternatively suggest an alternative way for me to extract all/certain reads from a .cram file into a new text file/array etc.? e.g. could I use samtools in the terminal, but I don't know how to extract just the query_sequence

Sorry in advance I'm very new to this kind of work

[Pierre Lindenbaum]

check the CRAM with samtools quickcheck the.cram && echo OK

[sacryt]

this returned OK, so does that mean there is no issue with the file itself?

[Pierre Lindenbaum]

at least the file is not truncated, now you should try:

samtools view -T /path/to/ref.fa in.cram > /dev/null && echo OK

if it's ok, then it could be a problem with RAP or with your code (shouldn't you provide a REFerence with the CRAM ?).

[sacryt]

I cannot find a reference genome within the RAP, but I can run samtools view x.cram without providing a reference directly, and this is currently working fine. I should note within the same directory there are also .cram.crai files, so is this why I don't need to provide a separate reference?

now I have no idea, does this just mean it's the pysam code that doesn't work specifically on this file for some reason?

[Pierre Lindenbaum]

and this is currently working fine

because the path is hard coded in the bam or samtools download the reference...

y there are also .cram.crai files

it's unrelated. crai files are only here if you need random-access.

[sacryt]

So shall I still try running the line you suggested and provide my own reference? e.g. using this https://ftp.ncbi.nlm.nih.gov/1000genomes/ftp/technical/reference/GRCh38_reference_genome/GRCh38_full_analysis_set_plus_decoy_hla.fa ?

[sacryt]

I have run the line you suggested and used the fasta reference I linked above, and this also ran successfully. Do you have any idea on what I can try to do to resolve this now? thanks for your help so far