realign bam files using bwa

Question

[ifudontmind_plzz]

Hi,

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])))
}

Answer 1

[karl.stamm]

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.

[alolex]

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.