RSEM error: RSEM currently does not support gapped alignments, sorry!

Question

[lfkhajavi]

hello,

i'm using STAR v2.5.2b and RSEM v1.2.31

STAR command:

for i in RNAseqFastq/*_R1_fastq.gz; do STAR --runThreadN 16 --genomeDir STAR/STAR_indices --sjdbGTFfile hg38_annotation.gtf --readFilesIn $i ${i%R1_fastq.gz}R2_fastq.gz --outFileNamePrefix
STAR/BAMfiles/$(basename $i fastq.gz) --readFilesCommand zcat --sjdbOverhang 50 --outFilterType BySJout --outFilterMultimapNmax 20 --outMultimapperOrder Random --alignSJoverhangMin 8 --alignSJDBoverhangMin 1
--outFilterMismatchNmax 999 --outFilterMismatchNoverLmax 0.04 --alignIntronMin 20 --alignIntronMax 1000000 --alignMatesGapMax 1000000 --quantMode TranscriptomeSAM --alignSoftClipAtReferenceEnds No
--outSAMstrandField intronMotif --outSAMmultNmax 1 --clip5pNbases 3 --outFilterIntronMotifs RemoveNoncanonical --outSAMtype BAM SortedByCoordinate; done

RSEM prepare reference:

rsem-prepare-reference -p 4 --gtf hg38_annotation.gtf hg38_genome.fa RSEM/hg38

RSEM calculate expression:

for i in STAR/*.toTranscriptome.out.bam; do rsem-calculate-expression --append-names --calc-ci --output-genome-bam --alignments --paired-end $i RSEM/hg38 RSEM/parALL/RNAseq_quals/$(basename $i
_Aligned.toTranscriptome.out.bam); done

ERROR: RSEM currently does not support gapped alignments, sorry!

I can't figure out where my problem lies... any help will be greatly appreciated. Thank you in advance, Leila

[GenoMax]

STAR is producing spliced alignment which RSEM does not like. @Alex Dobin (author of STAR) recommends this to turn off spliced alignments.

From this page.

However, please note that RSEM does * not * support gapped alignments. So make sure that your aligner does not produce alignments with intersions/deletions.

Answer 1

[Jake Warner]

The error is pretty verbose. RSEM does not support gapped alignment. So when you map to a genome (which has introns) you need to use an end to end aligner like bowtie or set up STAR appropriately. Try these options for STAR to prohibit gaps:

--alignEndsType EndToEnd --alignIntronMax 1 --alignIntronMin 2 --scoreDelOpen -10000 --scoreInsOpen -10000

[lfkhajavi]

Thank you... i'll give it a try.

[GenoMax]

Before you do this consider @Alex's answer (which I had included above).

if you generated the genome indexes with annotations, you will splicing to annotated junctions, as --alignIntronMax only controls the annotated junctions. You can either (i) use genome with annotations and --alignSJDBoverhangMin 999 (any number > read length) while mapping, or (ii) re-generate the genome without annotations.

[lfkhajavi]

hello,

changing option --alignSJDBoverhangMin 1 to --alignSJDBoverhangMin 999 did not fix the problem... unfortunately. re-mapping using option --alignEndsType EndToEnd didn't remedy the problem

the commands listed were working well with RSEM until I re-ran star with the addition of --outMultimapperOrder Random & --clip5pNbases 3 It's with this new alignment that I'm having issues with RSEM...

Could the option --clip5pNbases 3 be introducing soft- or hard- clippings?