Issue : fastq data cat into output/tmp/reads having error.

[csittk]

I had tried to use both fastq.gz and fastq files on the origami-alignment.
The first step origami-alignment unzip and cat the fastq files to the output/tmp/left-reads.fq and output/tmp/right-reads.fq had always come with an error

This is cutadapt 1.15 with Python 2.7.6
Command line parameters: -n 3 --overlap 10 -e 0 --discard-untrimmed -m 20 -a CTGCTGTCCG -A CTGCTGTCCG -o output/tmp/l_same_aa.fq -p output/tmp/r_same_aa.fq output/tmp/left_reads.fq output/tmp/right_reads.fq
Running on 1 core
Trimming 2 adapters with at most 0.0% errors in paired-end mode ...
cutadapt: error: In read named 'SRR6010260.33 WIGTC-HISEQ:1:1212:1644:2156 length=50': length of quality sequence (48) and length of read (50) do not match
This is cutadapt 1.15 with Python 2.7.6
Command line parameters: -n 3 --overlap 10 -e 0 --discard-untrimmed -m 20 -a CTGCTGTCAT -A CTGCTGTCAT -o output/tmp/l_same_bb.fq -p output/tmp/r_same_bb.fq output/tmp/left_reads.fq output/tmp/right_reads.fq
Running on 1 core
Trimming 2 adapters with at most 0.0% errors in paired-end mode ...
cutadapt: error: In read named 'SRR6010260.33 WIGTC-HISEQ:1:1212:1644:2156 length=50': length of quality sequence (48) and length of read (50) do not match
This is cutadapt 1.15 with Python 2.7.6
Command line parameters: -n 3 --overlap 10 -e 0 --discard-untrimmed -m 20 -a CTGCTGTCCG -A CTGCTGTCAT -o output/tmp/l_diff_ab.fq -p output/tmp/r_diff_ab.fq output/tmp/left_reads.fq output/tmp/right_reads.fq
Running on 1 core
Trimming 2 adapters with at most 0.0% errors in paired-end mode ...
cutadapt: error: In read named 'SRR6010260.33 WIGTC-HISEQ:1:1212:1644:2156 length=50': length of quality sequence (48) and length of read (50) do not match
This is cutadapt 1.15 with Python 2.7.6
Command line parameters: -n 3 --overlap 10 -e 0 --discard-untrimmed -m 20 -a CTGCTGTCAT -A CTGCTGTCCG -o output/tmp/l_diff_ba.fq -p output/tmp/r_diff_ba.fq output/tmp/left_reads.fq output/tmp/right_reads.fq
Running on 1 core
Trimming 2 adapters with at most 0.0% errors in paired-end mode ...
cutadapt: error: In read named 'SRR6010260.33 WIGTC-HISEQ:1:1212:1644:2156 length=50': length of quality sequence (48) and length of read (50) do not match
This is cutadapt 1.15 with Python 2.7.6
Command line parameters: -n 3 --overlap 10 -e 0 --discard-trimmed -m 20 -o output/tmp/l_neither.fq -p output/tmp/r_neither.fq SRR6010260_R1.fastq SRR6010260_R2.fastq

I checked with the raw fastq files the number of quality sequence and length of read for 'SRR6010260.33 WIGTC-HISEQ:1:1212:1644:2156 length=50' is correct. However, in the output/tmp/right-reads.fq files, it removed 2 base of quality sequence and result in error.

fastq raw data

@SRR6010260.33 WIGTC-HISEQ:1:1212:1644:2156 length=50
GTGCTGGGGTCCATGTGGGCCAGATGCCCTGGGCCCTGGGCAGGGCCAGG
+SRR6010260.33 WIGTC-HISEQ:1:1212:1644:2156 length=50
@d0@0<<<C/<D<C1<1<<1<@ghhc??CGC1E<@/1CGEFCDHHH@E??

output/tmp/right-reads.fq

@SRR6010260.33 WIGTC-HISEQ:1:1212:1644:2156 length=50
GTGCTGGGGTCCATGTGGGCCAGATGCCCTGGGCCCTGGGCAGGGCCAGG
+SRR6010260.33 WIGTC-HISEQ:1:1212:1644:2156 length=50
@d0@0<<<C/<D<C1<1<<1<@ghhc??CGC1E<@CGEFCDHHH@E??

It appears that the cat process of fastq files had removed "\1" from the fastq data?
Please check on this issue. I am not sure how it affect downstream alignment if the cutadapt step does not works properly.

Thanks.

[danielsday]

Hi csittk,

Thank you for reporting the error. This was a piece of code that tried to automatically remove /1 and /2 from the paired-end read names (since usually processed FASTQ files that had one of these two strings as a postfix in the read name) because leaving them in caused downstream problems when re-pairing the reads after aligning, quantification, etc. I had written code into the alignment steps to automatically trim these characters, but it didn't generalize well to other FASTQ files (like in this case). I had removed the code before releasing it, but it appears there was one place where this trimming code still existed in the code. I just uploaded a branch where this should be removed now.

Could you try downloading, installing and running the code from this branch to see if this fixes your error? https://github.com/younglab/origami/tree/read-hotfix

[csittk]

I am having trouble again with the BamHeader.h files during configuration and installation. I am able to configure using CPPFLAGS=' -I/usr/local/include/bamtools', but during make process it reported BamHeader.h not found.

Skipping the installation process, I try to use the binary origami-alignment and change the command line options KEEPPOSTFIX=no -> yes and I am able to run the alignment without having the previous error.

[danielsday]

Could you do me a favor and post your ./configure and make output? And can you post the output of "Makefile" here?

[csittk]

linux_bash:~/tools/origami$ ./configure CPPFLAGS="$CPPFLAGS -I/usr/local/include/bamtools/"

configure: Checking for software tools needed by origami
checking for cutadapt... cutadapt
checking for macs2... macs2
checking for bowtie... bowtie
checking for bedtools... bedtools
checking for R... R
checking for Rscript... Rscript
checking for samtools... samtools
checking for perl... perl
checking for g++... g++
checking whether the C++ compiler works... yes
checking for C++ compiler default output file name... a.out
checking for suffix of executables...
checking whether we are cross compiling... no
checking for suffix of object files... o
checking whether we are using the GNU C++ compiler... yes
checking whether g++ accepts -g... yes
configure: Checking for bedtools version
configure: Checking for C++ libraries needed by origami
checking how to run the C++ preprocessor... g++ -E
checking for grep that handles long lines and -e... /usr/local/bin/grep
checking for egrep... /usr/local/bin/grep -E
checking for ANSI C header files... yes
checking for sys/types.h... yes
checking for sys/stat.h... yes
checking for stdlib.h... yes
checking for string.h... yes
checking for memory.h... yes
checking for strings.h... yes
checking for inttypes.h... yes
checking for stdint.h... yes
checking for unistd.h... yes
checking BamReader.h usability... yes
checking BamReader.h presence... no
configure: WARNING: BamReader.h: accepted by the compiler, rejected by the preprocessor!
configure: WARNING: BamReader.h: proceeding with the compiler's result
checking for BamReader.h... yes
checking BamWriter.h usability... yes
checking BamWriter.h presence... no
configure: WARNING: BamWriter.h: accepted by the compiler, rejected by the preprocessor!
configure: WARNING: BamWriter.h: proceeding with the compiler's result
checking for BamWriter.h... yes
configure: Checking for R libraries needed by origami
checking if R library GenomicRanges is installed... yes
checking if R library matrixStats is installed... yes
configure: creating ./config.status
config.status: creating Makefile
config.status: creating bin/origami-alignment
config.status: creating bin/origami-analysis

linux_bash:~/tools/origami$ make

cd src && g++ -I/usr/local/include/bamtools/api/ -c mapped-reads-merge.cpp
mapped-reads-merge.cpp:1:27: fatal error: api/BamReader.h: No such file or directory
#include <api/BamReader.h>
^
compilation terminated.
make: *** [bamexes] Error 1

This is the error.

[danielsday]

OK, thank you. This works for me:

dsday@tak4 ~/dsday/temp/origami$ ./configure CPPFLAGS="$CPPFLAGS -I/usr/local/include/bamtools"
configure: Checking for software tools needed by origami
checking for cutadapt... cutadapt
checking for macs2... macs2
checking for bowtie... bowtie
checking for bedtools... bedtools
checking for R... R
checking for Rscript... Rscript
checking for samtools... samtools
checking for perl... perl
checking for g++... g++
checking whether the C++ compiler works... yes
checking for C++ compiler default output file name... a.out
checking for suffix of executables...
checking whether we are cross compiling... no
checking for suffix of object files... o
checking whether we are using the GNU C++ compiler... yes
checking whether g++ accepts -g... yes
configure: Checking for bedtools version
configure: Checking for C++ libraries needed by origami
checking how to run the C++ preprocessor... g++ -E
checking for grep that handles long lines and -e... /bin/grep
checking for egrep... /bin/grep -E
checking for ANSI C header files... yes
checking for sys/types.h... yes
checking for sys/stat.h... yes
checking for stdlib.h... yes
checking for string.h... yes
checking for memory.h... yes
checking for strings.h... yes
checking for inttypes.h... yes
checking for stdint.h... yes
checking for unistd.h... yes
checking api/BamReader.h usability... yes
checking api/BamReader.h presence... yes
checking for api/BamReader.h... yes
checking api/BamWriter.h usability... yes
checking api/BamWriter.h presence... yes
checking for api/BamWriter.h... yes
configure: Checking for R libraries needed by origami
checking if R library GenomicRanges is installed... yes
checking if R library matrixStats is installed... yes
configure: creating ./config.status
config.status: creating Makefile
config.status: creating bin/origami-alignment
config.status: creating bin/origami-analysis

dsday@tak4 ~/dsday/temp/origami$ make
cd src && g++ -I/usr/local/include/bamtools -c mapped-reads-merge.cpp
cd src && g++ -o mapped-reads-merge mapped-reads-merge.o -lbamtools

What I don't understand on your version is this part (bolded):

cd src && g++ -I/usr/local/include/bamtools/api/ -c mapped-reads-merge.cpp

Did you happen to edit the generated Makefile to change this path? When I run the above ./configure, I get -I/usr/local/include/bamtools in that same position of the Makefile, so I am thinking either you are running ./configure differently or you edited the Makefile. Right now, it seems that you are telling g++ to look in the /usr/local/include/bamtools/api/ folder for the api/BamReader.h (which is the path in the #include preprocessing statement), and g++ is trying to find /usr/local/include/bamtools/api/api/BamReader.h which doesn't exist. If you point it to look in /usr/local/include/bamtools/, I think then it should work.

[csittk]

Thanks. The installation issue resolve. Still the issue persist with the new update origami-alignment using the ablinker-test file.

origami-alignment /home/tzeking/data/bowtie1_chr_hg19/bowtie1_chr_hg19 left-test.fq right-test.fq

Launching origami...
Running in long-linker trimming mode
This is cutadapt 1.15 with Python 2.7.6
Command line parameters: -n 3 -m 20 --overlap 10 -a forward=ACGCGATATCTTATCTGACT -a reverse=AGTCAGATAAGATATCGCGT -o output/tmp/l_t1.fq --untrimmed-output output/tmp/l_nt1.fq -p output/tmp/r_t1.fq --untrimmed-paired-output output/tmp/r_nt1.fq output/tmp/left_reads.fq output/tmp/right_reads.fq
Running on 1 core
Trimming 2 adapters with at most 10.0% errors in paired-end legacy mode ...
cutadapt: error: In read named 'mES-cohesin-rep1.28 HWI-ST193:495:C16HUACXX:6:1101:1752:2215 length=50': length of quality sequence (48) and length of read (50) do not match
This is cutadapt 1.15 with Python 2.7.6
Command line parameters: -n 3 -m 20 --overlap 10 -a forward=ACGCGATATCTTATCTGACT -a reverse=AGTCAGATAAGATATCGCGT -o output/tmp/r_t2.fq --untrimmed-output output/tmp/r_nt2.fq -p output/tmp/l_t2.fq --untrimmed-paired-output output/tmp/l_nt2.fq output/tmp/r_nt1.fq output/tmp/l_nt1.fq
Running on 1 core
Trimming 2 adapters with at most 10.0% errors in paired-end legacy mode ...
Finished in 0.00 s (152 us/read; 0.39 M reads/minute).

=== Summary ===

Total read pairs processed: 27
Read 1 with adapter: 0 (0.0%)
Read 2 with adapter: 0 (0.0%)
Pairs that were too short: 0 (0.0%)
Pairs written (passing filters): 27 (100.0%)

Total basepairs processed: 2,700 bp
Read 1: 1,350 bp
Read 2: 1,350 bp
Total written (filtered): 2,700 bp (100.0%)
Read 1: 1,350 bp
Read 2: 1,350 bp

=== First read: Adapter forward ===

Sequence: ACGCGATATCTTATCTGACT; Type: regular 3'; Length: 20; Trimmed: 0 times.

=== First read: Adapter reverse ===

Sequence: AGTCAGATAAGATATCGCGT; Type: regular 3'; Length: 20; Trimmed: 0 times.

This is cutadapt 1.15 with Python 2.7.6
Command line parameters: -n 3 -m 20 --overlap 10 -a forward=ACGCGATATCTTATCTGACT -a reverse=AGTCAGATAAGATATCGCGT -o output/tmp/r_t3.fq --untrimmed-output output/tmp/r_nt3.fq -p output/tmp/l_t3.fq --untrimmed-paired-output output/tmp/l_nt3.fq output/tmp/r_t1.fq output/tmp/l_t1.fq
Running on 1 core
Trimming 2 adapters with at most 10.0% errors in paired-end legacy mode ...
No reads processed! Either your input file is empty or you used the wrong -f/--format parameter.
Aligning reads

/usr/local/bin/origami-alignment: line 69: 7365 Aborted (core dumped) bowtie -n 1 -m 1 -p 6 --sam /home/tzeking/data/bowtie1_chr_hg19/bowtie1_chr_hg19 output/tmp/left_kept.fq > output/tmp/left_kept.sam 2> output/logs/first_read_alignment.txt
/usr/local/bin/origami-alignment: line 69: 7369 Aborted (core dumped) bowtie -n 1 -m 1 -p 6 --sam /home/tzeking/data/bowtie1_chr_hg19/bowtie1_chr_hg19 output/tmp/right_kept.fq > output/tmp/right_kept.sam 2> output/logs/second_read_alignment.txt
/usr/local/bin/origami-alignment: line 69: /usr/local/bin/../bin/mapped-reads-merge: No such file or directory
Filtering out duplicated PETs
[E::hts_open_format] fail to open file 'output/mapped_reads.bam'
samtools view: failed to open "output/mapped_reads.bam" for reading: No such file or directory
Processed 0 PETs, excluded 0 PETs because of duplication
Calling peaks
/usr/local/bin/origami-alignment: line 713: macs: command not found
INFO @ Tue, 06 Feb 2018 04:09:56:
Command line: callpeak -t output/mapped_reads.bam -n macs2 -q 0.05 -g hs --nomodel --extsize 100 --outdir output
ARGUMENTS LIST:
name = macs2
format = AUTO
ChIP-seq file = ['output/mapped_reads.bam']
control file = None
effective genome size = 2.70e+09
band width = 300
model fold = [5, 50]
qvalue cutoff = 5.00e-02
Larger dataset will be scaled towards smaller dataset.
Range for calculating regional lambda is: 10000 bps
Broad region calling is off
Paired-End mode is off

INFO @ Tue, 06 Feb 2018 04:09:56: #1 read tag files...
INFO @ Tue, 06 Feb 2018 04:09:56: #1 read treatment tags...
Traceback (most recent call last):
File "/usr/local/bin/macs2", line 5, in
pkg_resources.run_script('MACS2==2.1.1.20160309', 'macs2')
File "/usr/local/lib/python2.7/dist-packages/distribute-0.6.10-py2.7.egg/pkg_resources.py", line 461, in run_script
self.require(requires)[0].run_script(script_name, ns)
File "/usr/local/lib/python2.7/dist-packages/distribute-0.6.10-py2.7.egg/pkg_resources.py", line 1201, in run_script
exec script_code in namespace, namespace
File "/usr/local/lib/python2.7/dist-packages/MACS2-2.1.1.20160309-py2.7-linux-x86_64.egg/EGG-INFO/scripts/macs2", line 617, in

File "/usr/local/lib/python2.7/dist-packages/MACS2-2.1.1.20160309-py2.7-linux-x86_64.egg/EGG-INFO/scripts/macs2", line 57, in main

File "build/bdist.linux-x86_64/egg/MACS2/callpeak_cmd.py", line 73, in run
File "build/bdist.linux-x86_64/egg/MACS2/callpeak_cmd.py", line 395, in load_tag_files_options
File "MACS2/IO/Parser.pyx", line 58, in MACS2.IO.Parser.guess_parser (MACS2/IO/Parser.c:2786)
File "MACS2/IO/Parser.pyx", line 81, in MACS2.IO.Parser.guess_parser (MACS2/IO/Parser.c:2413)
File "MACS2/IO/Parser.pyx", line 774, in MACS2.IO.Parser.BAMParser.init (MACS2/IO/Parser.c:11146)
File "/usr/lib/python2.7/gzip.py", line 34, in open
return GzipFile(filename, mode, compresslevel)
File "/usr/lib/python2.7/gzip.py", line 94, in init
fileobj = self.myfileobj = builtin.open(filename, mode or 'rb')
IOError: [Errno 2] No such file or directory: 'output/mapped_reads.bam'
Failed to call MACS1 peaks correctly, see error messages

I think turn on the KEEPPOSTFIX can skip this issue, but the logic error still there.