Hi,
After running Kneaddata with Bowtie2 on paired-end data, the output I’m getting from the final output seems to be unpaired (the first read has over 9x the amount of reads as the second). I’m curious to know if there’s a way to force paired-end reads in the analysis and throw out any reads which are unpaired.
The command I ran was the following:
Thanks!
Hi, Thanks for the post. Kneaddata should by default track read pairs if a pair of input files are provided. You should see pair output files (with the same number of reads) and orphan files. I think in your case kneaddata is possibly having an issue tracking the pairs due to sequence identifiers of an unexpected format. Can you check to see if there are spaces in the sequence ids or possibly they are missing the read number?
On our end we will work on updating kneaddata to catch the case where sequence ids of an unexpected format are provided and throw an informative error message. Sorry for the confusion.
Thank you,
Hi Lauren,
Thank you so much for you reply. I changed the sequence identifiers and sure enough, that solved the issue. Thanks again!
Best,
@lauren.j.mciver I have a follow up question to this thread. Kneaddata v.0.7.4 is not identifying the paired ends when run in the biobakery_workflow v.3.0.0-alpha.7. My file names are “ABCD_S70_R1.fastq.gz”, so I added the command line options --pair-identifier “_R1” since it does not follow the “.R1” format in the user guide on github.
My sequence identifiers in the fastq look like this: @EAS139:136:FC706VJ:2:2104:15343:197393 1:N:18:ATCACG. Would the space between the first bit and the 1:N... cause the workflow to miss the pair?
The other oddity is that when I run kneaddata on the same data files outside of the workflow and without the --pair-identifier flag, it runs successfully and correctly merges R1/R2 into a single fastq.
Do you have any guidance? I’m not sure how to troubleshoot given that kneaddata independently runs correctly, but fails to merge the pairs when run as part of the workflow.
Hi @ewissel , Thank you for the detailed post. I just checked the Kneaddata code for the latest version and it should catch your sequence identifier format with the space plus “:1”. The --pair-identifier flag is just used for the bioBakery workflows so the tool can pick up the paired files to pass on to Kneaddata. If you don’t use that flag and the identifier does not match the default then the workflow will process the reads as single end. Have you tried updating to the latest version of Kneaddata v0.10.0?
Thank you,
Thanks, Lauren! Maybe the kneaddata version is the issue. I download the workflow via conda, and the current version of kneaddata on conda is 0.7.4. Is it possible to update the conda install to v0.10.0?
Hi @ewissel , I pushed the latest version of kneaddata to conda. Please try it out and let me know if it resolves the issue you are seeing.
Thank you,
Thanks, Lauren!
When I run conda update -c biobakery kneaddata, I still get v0.7.4. I know that bioconda and biobakery channels have kneaddata, but I thought that directing to the biobakery channel would ensure conda looks there to update kneaddata.
Do you have any advice on this?
Hi Emily, I think conda will pick the best version based on your current environment which is not always the latest version. Can you try adding the specific version into your command kneaddata=0.10.0 and see if that will get the latest one?
Thanks!
Hey Lauren,
I was successfully able to update kneaddata with conda install -c biobakery kneaddata=0.10.0 (looks like conda upgrade doesn’t like version info). However, this did not resolve kneaddata not matching paired end files properly.
My files are name “sampleID_SX_R1_001.fastq.gz”, and for the workflow I have tried the following pair identifier arguments:
“_R1”
“_R1_001”
“_R1_”
“R1_”
“R1_001”
Is my identifier argument wrong? Should I be doing something differently for anadama2/the workflow?
Also, I am using biobakery_workflows v3.0.0-alpha.7, which is the latest on conda. Should I be using a different version?
Hi Emily, Thank you for the follow up. I am glad you were able to update to the latest kneaddata version. I think any of those pair identifiers should work with the workflow. If you could send me (feel free to send it directly to me) your log file I can dig in a bit more detail to see what might be going on.
Thank you,
CK_zhu
July 17, 2021, 1:59am
12
Hi lauren I will continue this problem
First I use kneaddata 0.7.2 it give me unpaired reads
singularity --debug exec
/software/kneaddata_0.7.2.sif kneaddata \
--remove-intermediate-output --threads 4 --bypass-trim \
--input ../raw_data/SRR527911_1.fastq.gz --input ../raw_data/SRR527911_2.fastq.gz \
--output ../temp/01_kneaddata --reference-db ../databases/kneaddata/human_genome \
--bowtie2-options "--very-sensitive --dovetail"
paired 1 737M and paired 2 82M are not same
-rw-r--r-- 1 ckzhu sample_lib 0 Jul 17 09:07 SRR527911_1_kneaddata_unmatched_2.fastq
-rw-r--r-- 1 ckzhu sample_lib 0 Jul 17 09:07 SRR527911_1_kneaddata_Homo_sapiens_bowtie2_unmatched_1_contam.fastq
-rw-r--r-- 1 ckzhu sample_lib 13K Jul 17 09:07 SRR527911_1_kneaddata_unmatched_1.fastq
-rw-r--r-- 1 ckzhu sample_lib 82M Jul 17 09:07 SRR527911_1_kneaddata_paired_2.fastq
-rw-r--r-- 1 ckzhu sample_lib 737M Jul 17 09:07 SRR527911_1_kneaddata_paired_1.fastq
-rw-r--r-- 1 ckzhu sample_lib 1.4K Jul 17 09:07 SRR527911_1_kneaddata_Homo_sapiens_bowtie2_unmatched_2_contam.fastq
-rw-r--r-- 1 ckzhu sample_lib 448 Jul 17 09:07 SRR527911_1_kneaddata_Homo_sapiens_bowtie2_paired_contam_2.fastq
-rw-r--r-- 1 ckzhu sample_lib 4.0K Jul 17 09:07 SRR527911_1_kneaddata_Homo_sapiens_bowtie2_paired_contam_1.fastq
-rw-r--r-- 1 ckzhu sample_lib 11K Jul 17 09:07 SRR527911_1_kneaddata.log
Then I use kneaddata 0.10.0 from docker images
singularity --debug exec
/software/kneaddata_0.10.0.sif kneaddata \
--remove-intermediate-output --threads 4 --bypass-trim \
--input ../raw_data/SRR527911_1.fastq.gz --input ../raw_data/SRR527911_2.fastq.gz \
--output ../temp/01_kneaddata --reference-db ../databases/kneaddata/human_genome \
--bowtie2-options "--very-sensitive --dovetail"
Decompressing gzipped file ...
Decompressing gzipped file ...
Reformatting file sequence identifiers ...
Reformatting file sequence identifiers ...
Initial number of reads ( /public/home/sample_lib/ckzhu/software/Snakemake_singularity/test/pipeline/temp/01_kneaddata/reformatted_identifiers8t6864hl_decompressed_zex5i232_SRR527911_1 ): 1921490.0
Initial number of reads ( /public/home/sample_lib/ckzhu/software/Snakemake_singularity/test/pipeline/temp/01_kneaddata/reformatted_identifiersrphkmryn_decompressed_gf1ann04_SRR527911_2 ): 1921490.0
Bypass trimming
Total reads after trimming ( /public/home/sample_lib/ckzhu/software/Snakemake_singularity/test/pipeline/temp/01_kneaddata/reformatted_identifiers8t6864hl_decompressed_zex5i232_SRR527911_1 ): 1921490.0
Total reads after trimming ( /public/home/sample_lib/ckzhu/software/Snakemake_singularity/test/pipeline/temp/01_kneaddata/reformatted_identifiersrphkmryn_decompressed_gf1ann04_SRR527911_2 ): 1921490.0
ERROR: Unable to write file: /public/home/sample_lib/ckzhu/software/Snakemake_singularity/test/pipeline/temp/01_kneaddata/reformatted_identifiers8t6864hl_decompressed_zex5i232_SRR527911_1
DEBUG [U=5169,P=28402] Master() Child exited with exit status 1
why I can’t write tmp file in that dir? should root to run singularity?
conda create -n kneaddata kneaddata=0.10.0
source activate kneaddata
kneaddata \
--remove-intermediate-output --threads 4 --bypass-trim \
--input ../raw_data/SRR527911_1.fastq.gz --input ../raw_data/SRR527911_2.fastq.gz \
--output ../temp/01_kneaddata --reference-db ../databases/kneaddata/human_genome \
--bowtie2-options "--very-sensitive --dovetail"
Decompressing gzipped file ...
Decompressing gzipped file ...
Reformatting file sequence identifiers ...
Reformatting file sequence identifiers ...
Initial number of reads ( /public/home/sample_lib/ckzhu/software/Snakemake_singularity/test/pipeline/temp/01_kneaddata/reformatted_identifiersjsysc7lb_decompressed_sff5u2hr_SRR527911_1 ): 1921490.0
Initial number of reads ( /public/home/sample_lib/ckzhu/software/Snakemake_singularity/test/pipeline/temp/01_kneaddata/reformatted_identifiersf3j5gvzj_decompressed_hab_axlq_SRR527911_2 ): 1921490.0
Bypass trimming
Total reads after trimming ( /public/home/sample_lib/ckzhu/software/Snakemake_singularity/test/pipeline/temp/01_kneaddata/reformatted_identifiersjsysc7lb_decompressed_sff5u2hr_SRR527911_1 ): 1921490.0
Total reads after trimming ( /public/home/sample_lib/ckzhu/software/Snakemake_singularity/test/pipeline/temp/01_kneaddata/reformatted_identifiersf3j5gvzj_decompressed_hab_axlq_SRR527911_2 ): 1921490.0
ERROR: Unable to write file: /public/home/sample_lib/ckzhu/software/Snakemake_singularity/test/pipeline/temp/01_kneaddata/reformatted_identifiersjsysc7lb_decompressed_sff5u2hr_SRR527911_1
I answered this elsewhere on the forum , but tldr: The sequence names apparently must have something like .R1. or .R2. as the input file name.
Hello, Thank you for the detailed post and sorry for the slow response. I don’t think you need to run as root with singularity. It looks like it writes a couple of files before it fails. Is it possible you are running out of disk space? Kneaddata can use a bit of disk space (up to 4x the original input size) if if needs to decompress input files and reformat the sequence identifiers. I think in the first case with the older kneaddata version it is likely having an issue tracking the pairs which should be fixed in the newer version.
Thank you,
Hi Lauren,
And the outputs are
I am not sure if this is a sequence identifiers issue or something else.
The kneadata version is kneaddata v0.12.0
I can also send you my log file if that helps with answering my question.
Thank you!!!
This is the log file.
02/22/2023 11:23:51 PM - kneaddata.utilities - INFO: Decompressing gzipped file …
02/22/2023 11:26:09 PM - kneaddata.utilities - DEBUG: Checking output file from Trimmomatic : /panfs/jay/groups/29/gallaher/jiang329/practice/kneaddata_out/10co_S53.R1_kneaddata.trimmed.1.fastq
02/22/2023 11:28:35 PM - kneaddata.utilities - INFO: Total memory = 503.452335358 GB
02/22/2023 11:28:35 PM - kneaddata.utilities - INFO: Total memory = 503.452335358 GB
Thanks for the detailed error post. Kneadata only expects periods for file extensions so it is getting a bit confused in naming the output file with the “.R1” and “.R2”. If you would replace the “.” with a “_” in the file names that should hopefully fix the issue with the output file name.
Thank you,
Hello,
kneaddata --input1 S1_R1.fq.gz --input2 S1_R2.fq.gz --threads 10 --trimmomatic ~/software/Trimmomatic-0.33/dist/jar/ -db /home/yask/reference_data/knead_database/mouse --output knead_output
here are the examples of sequence headers:@V350094545L2C001R0020000295 :0:0:0:0 1:N:0:GCGATCTA_TCGCCTTA@V350094545L2C001R0020000295 :0:0:0:0 2:N:0:GCGATCTA_TCGCCTTA
Here is the log:
05/19/2023 02:09:55 PM - kneaddata.utilities - INFO: Decompressing gzipped file …clean.fastq --al-single /home/yask/raw_data/C3_ko_microbiome/all_files/knead_output/S1_R1_kneaddata_mouse_C57BL_6NJ_bowtie2_unmatched %_contam.fastq -S /dev/null
It would be great if you could help me with that!
Regards
Artem
1 Like
gmark
July 3, 2023, 1:21am
19
Hi Artem,
I’m running into the same issue you are regarding the unpaired reads. Did you ever find a solution?
Best,
Hi Artem,
I also have the exact same problem and am stumped as to a solution – did you end up finding one?
Thanks!
1 Like
