Porechop output?

Question

[giulia.trauzzi]

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=[e-mail]

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

[Michael]

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

[giulia.trauzzi]

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

[Michael]

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.

[GenoMax]

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.

[GenoMax]

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

[giulia.trauzzi]

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

[GenoMax]

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

Answer 1

[shelkmike]

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.

[giulia.trauzzi]

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