Good evening, everyone.

I've been working on a pipeline to convert FASTQ files obtained from RNA-seq into VCF files, but I'm having trouble getting it to work properly. I successfully trimmed the FASTQ files using “trimmomatic” and generated the BAM files. However, I'm encountering errors when trying to generate the VCF file from that BAM file.

The steps are as follows:

1. Create the reference index (.fai and .dict files)
2. Sort the BAM file by read name (queryname)
3. Identify and mark PCR-duplicated reads (create .markdup.bam and .metrics.txt)
4. Remove unmapped reads
5. Re-sort the data in chromosome order (seat-standardized)
6. Create a CSI index file for the BAM file (.sort.csi file)
7. Finally, detect individual sample variants (SNPs) and create a VCF file

The .markdup.bam file, .sort.bam file, and its CSI index file appear to have been created correctly. However, when attempting to execute the final step 6, an error occurs where GATK fails to detect the created CSI index file. Below are the respective commands.

1~6

cmd = (
f"gatk SortSam --java-options '-Xmx32G' "
f"-I {input_bam} "
f"-O {sorted_bam} "
f"-SO queryname "
f"--VALIDATION_STRINGENCY SILENT")

gatk_call = (
f"gatk MarkDuplicates --java-options '-Xmx32G' "
f"-I {input_sorted} "
f"-O {output_bam} "
f"-M {metrics_file} "
f"--ASSUME_SORT_ORDER queryname "
f"--VALIDATION_STRINGENCY SILENT")

cmd_filter = f"samtools view -b -F 4 -o {filtered_bam} {markdup_bam}"

cmd_sort = (
f"gatk SortSam --java-options '-Xmx32G' "
f"-I {filtered_bam} "
f"-O {final_bam} "
f"-SO coordinate "
f"--CREATE_INDEX false " 
f"--VALIDATION_STRINGENCY SILENT")

cmd_index = f"samtools index -c -m 14 {final_bam} -@ 8"

gatk_call = (
f"gatk HaplotypeCaller --java-options '-Xmx32G' "
f"--reference {reference} "
f"--emit-ref-confidence GVCF "
f"--input {input_sort_bam} "
f"--read-index {input_csi} "
f"--output {output_vcf} "
f"--smith-waterman FASTEST_AVAILABLE "
f"--native-pair-hmm-threads 22").

Apologies for the length, but I'd appreciate any help.