realign bam files using bwa
1
1
Entering edit mode
5.8 years ago

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])))
}
sequencing alignment R • 2.6k views
ADD COMMENT
2
Entering edit mode
5.8 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.

ADD COMMENT
0
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.

ADD REPLY

Login before adding your answer.

Traffic: 2180 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