How to subset an SRA file for a single chromosome?

Question

I used prefetch to get the Pacbio reads of chicken from the SRA database. I want to align these reads against a reference genome, but not all the reads. I am only interested in a particular region on Chromosome Z. Is there a way that I can subset the downloaded SRA file?

I am trying to use blast, but don't know where to start. How should I prepare the SRA file from prefetch for the blast? I have currently started fasterq-dump for the downloaded SRA file.

Or should I subset after aligning the entire SRA reads data with the reference genome, and in that case what tools do I use?

Answer 1

After dumping out the reads using fasterq-dump, you'll need to align them first and then extract those that map to your region of interest. I think minimap2 is an excellent choice. If you have HiFi reads, you can just use something like:

minimap2 -ax map-hifi ref.fa pacbio-ccs.fq.gz > aln.sam

Then, using samtools, coordinate-sort and index the sorted alignments:

samtools sort -o aln.sorted.bam aln.sam
samtools index aln.sorted.bam

With an indexed BAM, you can then view the alignments that overlap your region of interest:

samtools view aln.sorted.bam chrZ:1-10000

If you instead write these to a file, you can then convert them using samtools fasta. I think something like the following should suffice:

samtools view -bo chrZ_1_10000.bam aln.sorted.bam chrZ:1-10000
samtools fasta -0 chrZ_1_10000.fa -o /dev/null chrZ_1_10000.bam

You could then upload the FASTA file to the NCBI Web BLAST server.