search

deweylab / RSEM

RSEM: accurate quantification of gene and isoform expression from RNA-Seq data
http://deweylab.biostat.wisc.edu/rsem/
GNU General Public License v3.0
408 stars 118 forks source link

rsem-parse-alignments error #30

Closed jwisrael closed 7 years ago

jwisrael commented 7 years ago

I am using RSEM v.1.2.31 and I’ve run the following:

rsem/rsem-prepare-reference --num-threads 8 --star --gtf /ref/gencode.v19.annotation.hs37d5_chr.gtf /ref/hs37d5.fa /result/rsem_ref

rsem/rsem-calculate-expression --paired-end --star --num-threads 5 /fastq/20150507_HBR_2_mRNA_L1.LB2_1.clipped.fastq /fastq/20150507_HBR_2_mRNA_L1.LB2_2.clipped.fastq /ref/rsem_ref /output/20150507_HBR_2_mRNA_L1.LB2

The STAR mapping works fine, but the pipeline fails on rsem-parse-alignments (see below). Do you have any suggestions for getting past this step?

rsem-parse-alignments /ref/rsem_ref /output/20150507_HBR_2_mRNA_L1.LB2.temp/20150507_HBR_2_mRNA_L1.LB2 /output/20150507_HBR_2_mRNA_L1.LB2.stat/20150507_HBR_2_mRNA_L1.LB2 /output/20150507_HBR_2_mRNA_L1.LB2.temp/20150507_HBR_2_mRNA_L1.LB2.bam 3 -tag XM Parsed 1000000 entries Parsed 2000000 entries Parsed 3000000 entries Parsed 4000000 entries Parsed 5000000 entries Parsed 6000000 entries Parsed 7000000 entries Parsed 8000000 entries Parsed 9000000 entries Parsed 10000000 entries Parsed 11000000 entries Parsed 12000000 entries Parsed 13000000 entries Parsed 14000000 entries Parsed 15000000 entries Parsed 16000000 entries Parsed 17000000 entries Warning: Detected a read pair whose two mates have different names--HWI:1:X:1:1216:4127:50166 and HWI:1:X:1:1216:4180:50180! Read HWI:1:X:1:1216:4127:50166: RSEM currently does not support gapped alignments, sorry! "rsem-parse-alignments /ref/rsem_ref /output/20150507_HBR_2_mRNA_L1.LB2.temp/20150507_HBR_2_mRNA_L1.LB2 /output/20150507_HBR_2_mRNA_L1.LB2.stat/20150507_HBR_2_mRNA_L1.LB2 /output/20150507_HBR_2_mRNA_L1.LB2.temp/20150507_HBR_2_mRNA_L1.LB2.bam 3 -tag XM" failed! Plase check if you provide correct parameters/options for the pipeline!

bli25wisc commented 7 years ago

Hi @jwisrael, can you use samtools view to extract the BAM record of "HWI:1:X:1:1216:4127:50166" and "HWI:1:X:1:1216:4180:50180" ?

jwisrael commented 7 years ago

Hi @bli25wisc here are the records:

HWI:1:X:1:1216:4127:50166 345 ENST00000262493.6 3092 255 77M268435404S * 0 0 AACTTAAAAAAAAAAGGAAAAAAGAAAAAAAAGCCCACGGGTCTTTGTAAG DA>CDBDDFHJJJIHFGIJJJIHBIJJJJIHBHGHIJIGHHHHFDEFFCCC NH:i:1 HI:i:1

HWI:1:X:1:1216:4180:50180 419 ENST00000467733.1 268 0 50M = 326 109 GAAATTTGTTCTCAAGCACAAAGTAGCACAGAAGAGGGAAGATGCTGTTT B@BDFFFFHHHFHHIJIGIIII?FFDHCHECHG3C>;DDCBG@@GD<?E? NH:i:6 HI:i:1 HWI:1:X:1:1216:4180:50180 339 ENST00000467733.1 326 0 51M = 268 -109 GTGACTCGAAAACTTTCTGAAGCTGATAATAGAAAGATGTCTCGGAAGGAG FFHIHIJIJJJIJJJJJJIIHGIEIIJJJIJJIJJJIGHHFHDDDDFFCCC NH:i:6 HI:i:1 HWI:1:X:1:1216:4180:50180 419 ENST00000496050.1 240 0 50M = 298 109 GAAATTTGTTCTCAAGCACAAAGTAGCACAGAAGAGGGAAGATGCTGTTT B@BDFFFFHHHFHHIJIGIIII?FFDHCHECHG3C>;DDCBG@@GD<?E? NH:i:6 HI:i:2 HWI:1:X:1:1216:4180:50180 339 ENST00000496050.1 298 0 51M = 240 -109 GTGACTCGAAAACTTTCTGAAGCTGATAATAGAAAGATGTCTCGGAAGGAG FFHIHIJIJJJIJJJJJJIIHGIEIIJJJIJJIJJJIGHHFHDDDDFFCCC NH:i:6 HI:i:2 HWI:1:X:1:1216:4180:50180 419 ENST00000463503.1 216 0 50M = 274 109 GAAATTTGTTCTCAAGCACAAAGTAGCACAGAAGAGGGAAGATGCTGTTT B@BDFFFFHHHFHHIJIGIIII?FFDHCHECHG3C>;DDCBG@@GD<?E? NH:i:6 HI:i:3 HWI:1:X:1:1216:4180:50180 339 ENST00000463503.1 274 0 51M = 216 -109 GTGACTCGAAAACTTTCTGAAGCTGATAATAGAAAGATGTCTCGGAAGGAG FFHIHIJIJJJIJJJJJJIIHGIEIIJJJIJJIJJJIGHHFHDDDDFFCCC NH:i:6 HI:i:3 HWI:1:X:1:1216:4180:50180 419 ENST00000476217.1 276 0 50M = 334 109 GAAATTTGTTCTCAAGCACAAAGTAGCACAGAAGAGGGAAGATGCTGTTT B@BDFFFFHHHFHHIJIGIIII?FFDHCHECHG3C>;DDCBG@@GD<?E? NH:i:6 HI:i:4 HWI:1:X:1:1216:4180:50180 339 ENST00000476217.1 334 0 51M = 276 -109 GTGACTCGAAAACTTTCTGAAGCTGATAATAGAAAGATGTCTCGGAAGGAG FFHIHIJIJJJIJJJJJJIIHGIEIIJJJIJJIJJJIGHHFHDDDDFFCCC NH:i:6 HI:i:4 HWI:1:X:1:1216:4180:50180 419 ENST00000467789.1 300 0 50M = 358 109 GAAATTTGTTCTCAAGCACAAAGTAGCACAGAAGAGGGAAGATGCTGTTT B@BDFFFFHHHFHHIJIGIIII?FFDHCHECHG3C>;DDCBG@@GD<?E? NH:i:6 HI:i:5 HWI:1:X:1:1216:4180:50180 339 ENST00000467789.1 358 0 51M = 300 -109 GTGACTCGAAAACTTTCTGAAGCTGATAATAGAAAGATGTCTCGGAAGGAG FFHIHIJIJJJIJJJJJJIIHGIEIIJJJIJJIJJJIGHHFHDDDDFFCCC NH:i:6 HI:i:5 HWI:1:X:1:1216:4180:50180 419 ENST00000265044.2 335 0 50M = 393 109 GAAATTTGTTCTCAAGCACAAAGTAGCACAGAAGAGGGAAGATGCTGTTT B@BDFFFFHHHFHHIJIGIIII?FFDHCHECHG3C>;DDCBG@@GD<?E? NH:i:6 HI:i:6 HWI:1:X:1:1216:4180:50180 339 ENST00000265044.2 393 0 51M = 335 -109 GTGACTCGAAAACTTTCTGAAGCTGATAATAGAAAGATGTCTCGGAAGGAG FFHIHIJIJJJIJJJJJJIIHGIEIIJJJIJJIJJJIGHHFHDDDDFFCCC NH:i:6 HI:i:6

bli25wisc commented 7 years ago

Hi @jwisrael , I think this is a STAR aligner bug. The first alignment has CIGAR string "77M268435404S", which should not appear for transcript-based alignments. If you are not using the latest version of STAR, maybe you can consider update your program.

Hope it helps, Bo

jwisrael commented 7 years ago

Hi @bli25wisc yes I think that fixed the problem. Thanks!