1
$\begingroup$

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?

$\endgroup$
0

1 Answer 1

3
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ "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? $\endgroup$ Jan 13, 2023 at 11:51
  • $\begingroup$ @venkateshwar Typo, sorry! Yes, or likely just a single FASTQ file since these are PacBio long reads. $\endgroup$
    – Steve
    Jan 13, 2023 at 11:54

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.