2
$\begingroup$

I am trying to do adapter trimming, alignment and sorting for a range of large scale paired end fastq files. The code I am using is given below:

#!/bin/bash
# Do Trimming, Alignment according to Config file

#SBATCH --job-name=T_Align_single
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --cpus-per-task=32
#SBATCH --partition=batch
#SBATCH --account=markay
#SBATCH --time=72:00:00
#SBATCH --mem=16G
#SBATCH --mail-type=TIME_LIMIT_80,BEGIN,END
#SBATCH --mail-user=<email-redacted>

module load bowtie2/2.4.1 
module load samtools/1.9

mkdir -p Mapping/

config=${INPUT_FILE}

file1=$(cat $config | awk '$1 ~ /fastq/ {print $1}')
file2=$(cat $config | awk '$1 ~ /fastq/ {print $2}')
file3=$(cat $config | awk '$1 ~ /index=/ {print $1}'| sed 's/ref_index=//')


#adapter trimming
trimfile1=Mapping/"${file1%.fastq}".trim.fastq
trimfile2=Mapping/"${file2%.fastq}".trim.fastq

 java -jar /labs/markay/Aranyak/Adriana_intro_CnT/Cnt_Test/Test-analysis/fastq_processing/Trimmomatic-0.39/trimmomatic-0.39.jar  PE -phred33 $file1 $file2 \
   $trimfile1 Mapping/"${file1%.fastq}".unpaired.fastq \
   $trimfile2 Mapping/"${file2%.fastq}".unpaired.fastq \
   LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15 \
   MINLEN:20 ILLUMINACLIP: /labs/markay/Aranyak/Adriana_intro_CnT/adapter.fa:2:30:10 -threads ${SLURM_CPUS_PER_TASK}

# mapping using bowtie2
sortbam=Mapping/"${file1%_R1_001.fastq}".sort.bam
sample_ID=$(cat $config | awk '$1 ~ /fastq/ {print $4}')
bowtie2 -p ${SLURM_CPUS_PER_TASK} --local --very-sensitive-local --no-unal \
   --no-mixed --no-discordant --phred33 -I 10 -X 700 \
   -x $file3 \
   -1 $trimfile1 -2 $trimfile2 --rg-id $sample_ID | samtools view -Sb - | samtools sort -@${SLURM_CPUS_PER_TASK} -o $sortbam

#Index the sorting by coordinate bam file
samtools index ${sortbam}

My original fastq files before adapter trimming looked like this

R3-LiD-M-LK_S16_L002_R2_001.fastq.gz
R3-LiD-M-LK_S16_L002_R1_001.fastq.gz

However after adapter trimming I removed the lane names and my file looks like this

R3-LiD-M-LK_S16_R2_001.fastq.gz
R3-LiD-M-LK_S16_R1_001.fastq.gz

The code that worked for the original file is not generating a proper bam file at the end though I am getting a list of trimmed fastq files. The error message I am getting is:

[E::hts_open_format] Failed to open file R1-LiD-H-DJ_S3.sort.rmdup.bam
samtools index: failed to open "R1-LiD-H-DJ_S3.sort.rmdup.bam": No such file or directory
[E::hts_open_format] Failed to open file R1-LiD-H-DJ_S3.sort.bam
The file 'R1-LiD-H-DJ_S3.sort.bam' does not exist
[E::hts_open_format] Failed to open file R1-LiD-H-DJ_S3.sort.rmdup.bam
The file 'R1-LiD-H-DJ_S3.sort.rmdup.bam' does not exist
R1-LiD-H-DJ_S3 DONE.

How should I change the code to generate proper alignment and bam files?

$\endgroup$
4
  • 3
    $\begingroup$ This should be a simple learning exercise for you - all that's changed is that the lane number part is gone, but your error doesn't seem related to that at all, and no variables use the lane number part anywhere. Add a set -x after the #SBATCH lines and look at the job's STDERR file to see the exact commands being executed. $\endgroup$
    – Ram RS
    Sep 19, 2022 at 22:00
  • 1
    $\begingroup$ Or just print out your variables so you can see what is happening. You are trying to open a file named R1-LiD-H-DJ_S3.sort.bam and a file named R1-LiD-H-DJ_S3.sort.rmdup.bam and neither of them exist. You also don't mention rmdub anywhere in your code, is this really the code you are running? $\endgroup$
    – terdon
    Sep 21, 2022 at 16:40
  • $\begingroup$ Thanks for your help, I have been able to trouble shoot the error with the naming of the files. $\endgroup$ Sep 21, 2022 at 23:54
  • $\begingroup$ Hi @AranyakGoswami, can you please write up your solution as an answer to this question (and make it an accepted answer)? That will help the StackExchange bots do the right thing with regards to not hunting for more answers to an already-solved question. $\endgroup$
    – gringer
    Oct 13, 2022 at 23:38

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Browse other questions tagged or ask your own question.