RSEM error: RSEM currently does not support gapped alignments, sorry!
1
1
Entering edit mode
4.5 years ago
lfkhajavi ▴ 30

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

RSEM rsem-calculate-expression RNA-Seq • 3.5k views
ADD COMMENT
0
Entering edit mode

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.

ADD REPLY
2
Entering edit mode
4.5 years ago
Jake Warner ▴ 820

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
ADD COMMENT
0
Entering edit mode

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

ADD REPLY
1
Entering edit mode

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.

ADD REPLY
0
Entering edit mode

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?

ADD REPLY

Login before adding your answer.

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