realign bam files using bwa
Entering edit mode
7.6 years ago


I am trying to realign some of the bam files (paired end) using the following command. But I am getting error:

[fread] Unexpected end of file
bam.files <- list.files(pattern = ".bam$", full.names=T)
for( I in 1:length(bam.files)){
system(paste("bwa aln -t 22 ref.fa -b1 ", bam.files[i], ">", gsub(".bam", "_1.sai", bam.files[i])))
system(paste("bwa aln -t 22 ref.fa -b2 ", bam.files[i], ">", gsub(".bam", "_2.sai", bam.files[i])))

# For 2 .sai files (Paired-end)
sai.files <- list.files(pattern = ".sai$")
system(paste("bwa sampe", ref, sai.files[2 * I - 1], sai.files[2 * i], bam.files[i], bam.files[i], ">", gsub(".bam", ".sam", sai.files[2 * I - 1])))
alignment sequencing R • 3.1k views
Entering edit mode
7.6 years ago

You're doing an awful lot of things at once, embedding BWA in shell script in R in a loop.

Unexpected end of file means it tried to read something and had trouble. I don't know if that's BWA's error or R's. Try isolating those parts.

Does sai.files[2*i-1] exist?

Does system() even support >? That's not a system command, that's a bash command, and might be the problem, as R tries to read a file called > and can't find it. system() has its own method for piping results to files.

In general, R is not the right tool for this. BWA should be long-running and multithread and R's system() doesn't really like long running tasks. Get out of R and try it in BASH. Once you have a shell script that works, you can trigger it with R asynchronously.

Entering edit mode

Agreed! Pick your smallest file, run it in the bash terminal and if you are getting the same error post your command and the error here and we may be able to help you debug. One thing to check along the way is to ensure it is a BAM file and not a SAM file. Sometimes I get similar errors if I forget to convert SAM to BAM in a previous step but keep the BAM extension. You know it is SAM if you can read it in a text editor.


Login before adding your answer.

Traffic: 750 users visited in the last hour
Help About
Access RSS

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

Powered by the version 2.3.6