<p>After dumping out the reads using <code>fasterq-dump</code>, you'll need to align them first and then extract those that map to your region of interest. I think <a href="https://github.com/lh3/minimap2" rel="nofollow noreferrer">minimap2</a> is an excellent choice. If you have HiFi reads, you can just use something like:</p>
<pre><code>minimap2 -ax map-hifi ref.fa pacbio-ccs.fq.gz &gt; aln.sam
</code></pre>
<p>Then, using <a href="http://www.htslib.org/doc/samtools.html" rel="nofollow noreferrer">samtools</a>, coordinate-<a href="http://www.htslib.org/doc/samtools-sort.html" rel="nofollow noreferrer">sort</a> and <a href="http://www.htslib.org/doc/samtools-index.html" rel="nofollow noreferrer">index</a> the sorted alignments:</p>
<pre><code>samtools sort -o aln.sorted.bam aln.sam
samtools index aln.sorted.bam
</code></pre>
<p>With an indexed BAM, you can then <a href="http://www.htslib.org/doc/samtools-view.html" rel="nofollow noreferrer">view</a> the alignments that overlap your region of interest:</p>
<pre><code>samtools view aln.sorted.bam chrZ:1-10000
</code></pre>
<p>If you instead write these to a file, you can then convert them using samtools  <a href="http://www.htslib.org/doc/samtools-fasta.html" rel="nofollow noreferrer">fasta</a>. I think something like the following should suffice:</p>
<pre><code>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
</code></pre>
<p>You could then upload the FASTA file to the NCBI <a href="https://blast.ncbi.nlm.nih.gov/Blast.cgi" rel="nofollow noreferrer">Web BLAST</a> server.</p>

CC BY-SA 4.0
2023-01-13T11:51:59.970Z
"After dumping out the reads using prefetch". Do you mean after fasterq-dump ends? And what will be looking at once that ends, a set of fastq files?
venkatesh war
CC BY-SA 4.0
2023-01-13T11:54:16.827Z
@venkateshwar Typo, sorry! Yes, or likely just a single FASTQ file since these are PacBio long reads.
Steve
CC BY-SA 4.0
2023-01-13T10:39:44.550Z
2023-01-13T11:54:22.083Z
2023-01-13T11:54:22.083Z
Steve
<p>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?</p>
<p>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.</p>
<p>Or should I subset after aligning the entire SRA reads data with the reference genome, and in that case what tools do I use?</p>

<assembly><reference-genome><pacbio><sratoolkit>
How to subset an SRA file for a single chromosome?
# How to subset an SRA file for a single chromosome?

(Asked by: username_0 on 2023-01-13)

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?

**Tags:** assembly, reference-genome, pacbio, sratoolkit

---

## 1 answers

### Answer by username_1 on 2023-01-13. Score: 3 | Accepted answer

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.

> **Comments:**
> **username_0**: "After dumping out the reads using prefetch". Do you mean after fasterq-dump ends? And what will be looking at once that ends, a set of fastq files?  
> **username_1**: @venkateshwar Typo, sorry! Yes, or likely just a single FASTQ file since these are PacBio long reads.  

---

