This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Rsubread: Aligning multiple paired-end files

Hello brain trust, Thanks for reading, please help. I have 18 paired-end files (9 forward and 9 reverse, file extension is .fq, each has a unique name) and I am trying to align them using Rsubread. I am working through a hpc cluster, not in R directly. I slurm/sbatch to submit the R script.

If I align a set of paired-end files individually, everything works well. However, when I try to align the entire 18 files, I get an error message.

Here is my R script for the individual PE files: (this works fine)

align(index="mm10",readfile1="Control-1_R1_001_paired.fq",readfile2="Control-1_R2_001_paired.fq",input_format="FASTQ",output_file="rsubread.bam",output_format="BAM",nthreads=16)

Here is my R script for all of the files: (ps: I edited the path to the directory to "path to file" for de-identification reasons)

align(index="mm10",readfile1="/path to file/forward_reads_list",readfile2="/path to file/reverse_reads_list",input_format="FASTQ",output_file="rsubread.bam",output_format="BAM",nthreads=16)

The output_file "rsubread.bam" does not get listed in the working directory. For the slurm/sbatch, there is no error message as the .err file is empty. However, my output, .out is below: (Again, I edited the path to the directory to "path to file" for de-identification reasons).

 Function      : Read alignment (RNA-Seq)                                   ||
|| Input file 1  : /path to file/forward_ ... ||
|| Input file 2  : /path to file/reverse_ ... ||
|| Output file   : rsubread.bam (BAM)                                         ||
|| Index name    : mm10                                                       ||
||                                                                            ||
||                    ------------------------------------                    ||
||                                                                            ||
||                       Threads : 16                                         ||
||                  Phred offset : 33                                         ||
||       # of extracted subreads : 10                                         ||
||                Min read1 vote : 3                                          ||
||                Min read2 vote : 1                                          ||
||             Max fragment size : 600                                        ||
||             Min fragment size : 50                                         ||
||    Maximum allowed mismatches : 3                                          ||
||   Maximum allowed indel bases : 5                                          ||
|| # of best alignments reported : 1                                          ||
||                Unique mapping : no                                         ||
||                                                                            ||
\\===================== http://subread.sourceforge.net/ ======================//

//================ Running (22-Aug-2018 06:13:12, pid=36855) =================\\
||                                                                            ||
|| WARNING format issue in file '/path to file/ ... ||
||         The required format is : FASTQ or FASTA                            ||
||         The file format is unknown.                                        ||
|| A wrong format may result in wrong results or crash the program.           ||
|| Please refer to the manual for file format options.                        ||
|| If the file is in the correct format, please ignore this message.          ||
||                                                                            ||
|| WARNING format issue in file '/path to file/ ... ||
||         The required format is : FASTQ or FASTA                            ||
||         The file format is unknown.                                        ||
|| A wrong format may result in wrong results or crash the program.           ||
|| Please refer to the manual for file format options.                        ||
|| If the file is in the correct format, please ignore this message.

If you are still reading at this point, thank you :). Any ideas on how to solve this?

rna-seq alignment rsubread
• 11,015 views • link
written by imesi • imesi

Looks like subread does not accept a list of files (as input). Why no use a loop to go through the files with your single file command that works.

• 0 views • link
written by GenoMax GenoMax

Apologies for the delay. I was barred from responding because apparently their is a 5-post in 6 hours limit for newcomers. Had no idea. Here is my script. I slurm submit the script in unix, so it doesn't run on the command line directly.

#!/usr/bin/env Rscript
library(Rsubread)

reads1 <- list.files( path = "/path/to/file/forward_reads_list", pattern = "*_R1_001.fq$" )
reads2 <- list.files( path = "/path/to/file/reverse_reads_list", pattern = "*_R2_001.fq$" )
align(index="mm10",readfile1=reads1,readfile2=reads2,input_format="FASTQ",output_format="BAM",nthreads=16)

Thanks.

• 1 views • link
updated by GenoMax GenoMax • written by imesi • imesi

Hello, when I run this script I get the ERROR: unable to open file 'NA'. File name might be incorrect, or you do not have the permission to read the file.

ibrary(Rsubread)
reads1 <- list.files( path="/home/weiyan/Desktop/RNA_SEQ/files/fastq/trimmed/rif1", pattern = "*_1.fq.gz$" )
reads2 <- list.files( path="/home/weiyan/Desktop/RNA_SEQ/files/fastq/trimmed/rif1", pattern = "*_2.fq.gz$" )
all.equal(length(reads1),length(reads2))
align(index="mm10",readfile1=reads1, readfile2=reads2, input_format="gzFASTQ", 
output_file=paste(as.character(reads1)), output_format="BAM", nthreads=20)

Any suggestions? Thanks!

• 0 views • link
written by weiyanjia2008 • weiyanjia2008

2 answers

According to the documentation https://www.rdocumentation.org/packages/Rsubread/versions/1.22.2/topics/align

readfile1: a character vector including names of files that include sequence reads to be aligned. For paired-end reads, this gives the list of files including first reads in each library. File format is FASTQ/FASTA by default. See input_format option for more supported formats.

readfile2: a character vector giving names of files that include second reads in paired-end read data. Files included in readfile2 should be in the same order as their mate files included in readfile1 .NULL by default.

So you need to pass the names of files explicitly, not a list (or directory name).

• 1 views • link
written by Santosh Anand Santosh Anand

Thanks Santosh. Each of my files has a unique name, I guess I can write them in as *_R1_001.fq and *_R2_001.fq. I have tried that, but no luck.

• 1 views • link
written by imesi • imesi

No, you can not write *_R1_001.fq and *_R2_001.fq with the align() function. R is not a Unix shell, and glob expansions do not work generally.

• 1 views • link
written by h.mon h.mon

@hmon is very right.. R will not expand the globs (*, ? etc.) automatically. Instead you need to use list.files

list.files(pattern="*_R1_001.fq")
• 1 views • link
written by Santosh Anand Santosh Anand

I find using Rsubread cumbersome and it adds R overhead to an already heavy task. In addition, I had experienced Rsubread crashes, where subread-align is rock-solid.

You have to pass a vector with all file names, and either leave out output_file=, or pass a vector of output file names:

reads1 <- list.files( path = "/path/to/files", pattern = "*_1.fastq$" )
reads2 <- list.files( path = "/path/to/files", pattern = "*_2.fastq$" )
align(index="mm10",readfile1=reads1,readfile2=reads2,input_format="FASTQ",output_format="BAM",nthreads=16)
• 1 views • link
written by h.mon h.mon

Thank you h.mon I just ran your suggested script and this is the error:

Error: unexpected string constant in: " reads1 <- list.files( path = "/path/to/files", pattern = ""

• 1 views • link
written by imesi • imesi

You have to adjust the script to match file names you have. _R1_001.fq$ instead of _1.fastq$ as in example. Did you do that?

• 1 views • link
written by GenoMax GenoMax

Eeek! I did not. Thanks!

• 1 views • link
written by imesi • imesi

You anonymize too much your commands and error messages, and becomes difficult to help. You probably need pattern = "*_R1_001_paired.fq" and pattern = "*_R2_001_paired.fq", or something similar.

• 1 views • link
written by h.mon h.mon

True, h.mon, I am paranoid like that, probably because of healthcare privacy laws. I realize how that can be a hindrance in bioinformatics. I will do better. Thanks.

So I reran the script, and the error message is:

Error in align(index = "mm10", readfile1 = reads1, readfile2 = reads2, :

The number of input file names is different from the number of output file names.

Execution halted

This is the order of my files in the list forward_reads_list

Control-1_R1_001_paired.fq
Control-2_R1_001_paired.fq
Control-3_R1_001_paired.fq
IFNg-1_R1_001_paired.fq
IFNg-2_R1_001_paired.fq
IFNg-3_R1_001_paired.fq
M10_R1_001_paired.fq
M6_R1_001_paired.fq
M7_R1_001_paired.fq

reverse_reads_list

Control-1_R2_001_paired.fq
Control-2_R2_001_paired.fq
Control-3_R2_001_paired.fq
IFNg-1_R2_001_paired.fq
IFNg-2_R2_001_paired.fq
IFNg-3_R2_001_paired.fq
M10_R2_001_paired.fq
M6_R2_001_paired.fq
M7_R2_001_paired.fq
• 1 views • link
updated by h.mon h.mon • written by imesi • imesi

Can you copy / paste the command you run? Did you use the output_file= parameter?

• 1 views • link
written by h.mon h.mon

check if lengths of your lists are same and you are using latest Rsubread package. Try following (change the path appropriately):

reads1 <- list.files(path=".",pattern = "_R1_.*.fq$" )
reads2 <- list.files(path=".",pattern = "_R2_.*.fq$" )
all.equal(length(reads1),length(reads2))

library(Rsubread)
align(index="mm10",readfile1=reads1,readfile2=reads2,input_format="FASTQ",output_format="BAM",nthreads=6).

From the error "The file format is unknown." , it seems to me there is an error in one of your files. Validate your fastq files and/or make sure that no other files exist in the raw data directories other than fq/fastq files.

• 1 views • link
written by cpad0112 cpad0112

Thanks a gazillion cpad0112, it is running!!!! I am doing my happy dance! This is the script:

reads1 <- list.files(path=".",pattern = "_R1_.*.fq$" )
reads2 <- list.files(path=".",pattern = "_R2_.*.fq$" )
all.equal(length(reads1),length(reads2))
align(index="mm10",readfile1=reads1,readfile2=reads2,input_format="FASTQ",output_format="BAM",nthreads=16)

Thank you hmon, genomax and Santosh Anand, you guys are awesome!

• 1 views • link
updated by GenoMax GenoMax • written by imesi • imesi

good to know that things are working. If script is running in current directory, then you can as well remove, path="." to shorten further. By default, R searches in current directory. imesi

• 1 views • link
written by cpad0112 cpad0112

It is running in the current directory, but I won't touch it: if ain't broke..... I tried t modify the script to specify an output_file per the manual but it failed to run.

• 1 views • link
written by imesi • imesi

Okay. I agree with you. If it is not broken, don't fix it. Did the script run successfully till the end?

• 1 views • link
written by cpad0112 cpad0112

It is still running and will for a little while. I have 9 files and each is about 170 million reads. The first 2 files were completed in about 280 minutes each. I guess I am looking at a total run time of over a day.

• 1 views • link
written by imesi • imesi

Log in to answer this question.

More posts like this