Porechop output?
1
0
Entering edit mode
9 months ago

Hello everyone,

I have been trying to run porechop on my nanopore data from the Promethion. I am trying to remove the adapters from the ligation sequencing kit LSK-114 prior to mapping the reads against the reference genome using STAR (although I may use bowtie2 or bwa-mem as STAR keeps running out of time or memory).

This is my code for Porechop

#!/bin/bash #SBATCH --mem=80G #SBATCH --ntasks=20 #SBATCH --tasks-per-node=20 #SBATCH -t 12:00:00 #SBATCH -o OUT/01porechop.%J #SBATCH -e ERR/01porechop.%J #SBATCH --job-name=porechop #SBATCH --partition=c_compute_cg1 #SBATCH --account=scw1179 #SBATCH --mail-type=ALL #SBATCH --mail-user=xxx@cardiff.ac.uk

module load porechop/0.2.4 for f in ../resources/fastq/sub/*.fastq.gz do porechop -i $f --format fastq.gz -t 20 -o ../resources/fastq/sub/trimmed/${f}_trimmed_pore.fastq.gz done

I am running this on 19 GB of data and after 12 hours it ran out of time. Should it usually take that long? I am currently running it with increased time, but I am facing another issue, I do not have any output in the specified output folder. I have a temporary file in my bin folder for each fastq file but not the trimmed reads. Am I doing something wrong? Or does it mean that porechop is not finding any adapter in my reads?

Thanks,

Giulia

porechop Nanopore • 1.4k views
ADD COMMENT
2
Entering edit mode

Is it RNA-seq? Otherwise, I would use minimap2 for aligning long reads to the genome.

ADD REPLY
0
Entering edit mode

Hi, no it isn't. It is WGS from Nanopore. I did some more reading and saw that Minimap2 was better for this. The problem with minimap2 is that it is outputting a truncated sam file that I cannot convert into bam, sort nor index. I am really not sure why. I did the same with another dataset and all worked well. I am low-key wondering if the issue may be the quality of the data?

This is the error message I am getting when trying to convert

[E::sam_parse1] missing SAM header
[W::sam_read1] Parse error at line 3
[main_samview] truncated file.

I read others had this issue, but not sure how to solve it? Should I attempt with bwa-mem? or bowtie2?

Thanks,

Giulia

ADD REPLY
0
Entering edit mode

I have never experienced this problem, but I am using it mostly inside other pipelines, e.g. Quast. Can you make a new topic where you explain the error in detail? Possibly you are just a missing switch in minimap2 or samtools. Also, I have seen a lot of these errors happen when people run pipelines through screen in the wrong way.

ADD REPLY
0
Entering edit mode

I read others had this issue, but not sure how to solve it?

Did you use the -a option with minimap2. By default minimap2 outputs PAF format alignments.

ADD REPLY
1
Entering edit mode

If you can, use dorado instead. It will trim the adapters and will allow you to specify the kit.

ADD REPLY
0
Entering edit mode

Thanks, I will try it. :) Actually I did try, but my job keeps failing and creating empty files for some reason.

Would the --trim subcommand under basecaller? I just cannot find it in the version of dorado I am using (0.3.4).

Thanks

ADD REPLY
1
Entering edit mode

dorado is now in v.0.5.2. So definitely upgrade. Many new features.

ADD REPLY
1
Entering edit mode
9 months ago
shelkmike ★ 1.4k

Porechop is poorly parallelizable. Instead of running it with "-t 20" for files consecutively, it's better to run it for many files in parallel using "-t 1" for each.

ADD COMMENT
0
Entering edit mode

Thanks for the tip :) Will give it a try!

ADD REPLY

Login before adding your answer.

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