Rsubread doesn't output a .bai file
1
0
Entering edit mode
@fa745d19
Last seen 11 weeks ago
United Kingdom

I'm trying to align some Illumina amplicon sequencing, generated on a MiSeq, using Rsubread. It outputs .bam, .indel.vcf, and .summary files, but no .bai file, so I can't display the results in IGV. Here's how I got to this point:

  1. Download the chr5.fa.gz reference sequence from UCSC (https://hgdownload.soe.ucsc.edu/goldenPath/hg38/chromosomes/).

  2. Create an index using: buildindex(basename="reference.ch5", reference="chr5.fa")

  3. Align my paired end sequencing files:

reads1.ch5 <- file.path(dirPath, "MSH3repeat_S1_L001_R1_001.fastq.gz")

reads2.ch5 <- file.path(dirPath, "MSH3repeat_S1_L001_R2_001.fastq.gz")

align(index="reference.ch5", readfile1=reads1.ch5, readfile2=reads2.ch5, output_file="alignResults.ch5.BAM")

How do I get it to make a .bai file, and how do I then display the results in IGV? Thanks!

Rsubread align Alignment • 145 views
ADD COMMENT
2
Entering edit mode
@gordon-smyth
Last seen 1 hour ago
WEHI, Melbourne, Australia

bai files are only produced when reads are sorted by genomic position, which Rsubread::align will do if you set sortReadsByCoordinates = TRUE.

Rather than running align again, you could use samtools to sort the bam files. That does not require re-alignment.

ADD COMMENT
0
Entering edit mode

Worked perfectly, thanks!

ADD REPLY

Login before adding your answer.

Traffic: 267 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6